本文整理汇总了C++中rcp函数的典型用法代码示例。如果您正苦于以下问题：C++ rcp函数的具体用法？C++ rcp怎么用？C++ rcp使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了rcp函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc, char* argv[]) {
#ifdef HAVE_MPI
MPI::Init(argc, argv);
#endif
RCP<MxComm> myComm = rcp(new MxComm());
#if 0
#ifdef HAVE_MPI
MPI::Init(argc, argv);
//MPI_Init(argc, argv);
Epetra_MpiComm myComm(MPI_COMM_WORLD);
#else
Epetra_SerialComm myComm;
#endif
#endif
// input file method
#if 1
std::string inFile;
Teuchos::CommandLineProcessor cmdp(false, true);
cmdp.setOption("infile", &inFile, "XML format input file.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
if (inFile == "") {
std::cout << "Please specify an input file using --infile=your_file.mx\n";
exit(0);
}
// now read the input file with trilinos XML reader
Teuchos::XMLObject xmlObj(Teuchos::FileInputSource(inFile).getObject());
// get simulation dimension
int dim = atoi(MxUtil::XML::getAttr("dim", xmlObj).c_str());
if (dim < 1 or dim > 3) {
std::cout << "Simulation dimension invalid or not given, using 3D.\n";
dim = 3;
}
// get simulation type
std::string domain = MxUtil::XML::getAttr("domain", xmlObj).c_str();
if (domain != "frequency" and domain != "time") {
std::cout << "Simulation domain invalid or not given, using frequency-domain.\n";
domain = "frequency";
}
// create problem
MxProblem<1> * prob1d;
MxProblem<2> * prob2d;
MxProblem<3> * prob3d;
switch (dim) {
case 1:
prob1d = new MxProblem<1>(xmlObj, myComm);
prob1d->solve();
delete prob1d;
break;
case 2:
prob2d = new MxProblem<2>(xmlObj, myComm);
prob2d->solve();
delete prob2d;
break;
case 3:
prob3d = new MxProblem<3>(xmlObj, myComm);
prob3d->solve();
delete prob3d;
break;
}
#endif
#if 0
// epetra stuff test
MxMap map(10, 0, myComm);
Epetra_CrsMatrix mat(Copy, map, 0);
int ind = 2;
double val = 0;
mat.InsertGlobalValues(1, 1, &val, &ind);
ind = 3;
val = 4;
mat.InsertGlobalValues(1, 1, &val, &ind);
mat.FillComplete(map, map);
Epetra_Vector myvec(map);
myvec.Random();
std::cout << myvec;
mat.Apply(myvec, myvec);
std::cout << myvec;
Epetra_CrsMatrix copy(mat);
std::cout << mat;
MxUtil::Epetra::stripZeros(mat);
std::cout << mat;
//.........这里部分代码省略.........
开发者ID:achellies，项目名称:maxwell，代码行数:101，代码来源:maxwell.cpp
示例2: CreateHierarchy
virtual RCP<Hierarchy> CreateHierarchy() const {
return rcp(new Hierarchy());
}
开发者ID:rorypeck，项目名称:trilinos，代码行数:3，代码来源:MueLu_HierarchyManager.hpp
示例3: rcp
void CrsMatrixWrapper<ST>::solve(const Teuchos::ArrayView<ST>& x,
const Teuchos::ArrayView<const ST>& b,
escript::SolverBuddy& sb) const
{
typedef VectorType<ST> Vector;
RCP<Vector> X = rcp(new Vector(mat.getDomainMap(), 1));
RCP<Vector> B = rcp(new Vector(mat.getRangeMap(), b, b.size(), 1));
RCP<const Matrix> A = rcpFromRef(mat);
if (escript::isDirectSolver(sb.getSolverMethod())) {
RCP<DirectSolverType<Matrix,Vector> > solver(m_direct);
if (solver.is_null()) {
solver = createDirectSolver<Matrix,Vector>(sb, A, X, B);
m_direct = solver;
if (sb.isVerbose()) {
std::cout << "Using " << solver->description() << std::endl;
std::cout << "Performing symbolic factorization..." << std::flush;
}
solver->symbolicFactorization();
if (sb.isVerbose()) {
std::cout << "done\nPerforming numeric factorization..." << std::flush;
}
solver->numericFactorization();
if (sb.isVerbose()) {
std::cout << "done\n" << std::flush;
}
} else {
if (sb.isVerbose()) {
std::cout << "Using " << solver->description() << std::endl;
}
if (m_resetCalled) {
// matrix structure never changes
solver->setA(A, Amesos2::SYMBFACT);
m_resetCalled = false;
}
solver->setX(X);
solver->setB(B);
}
if (sb.isVerbose()) {
std::cout << "Solving system..." << std::flush;
}
solver->solve();
if (sb.isVerbose()) {
std::cout << "done" << std::endl;
RCP<Teuchos::FancyOStream> fos(Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)));
solver->printTiming(*fos, Teuchos::VERB_HIGH);
}
} else { // iterative solver
double t0 = Teuchos::Time::wallTime();
RCP<ProblemType<ST> > problem(m_solver);
if (problem.is_null()) {
problem = rcp(new ProblemType<ST>(A, X, B));
m_solver = problem;
RCP<OpType<ST> > prec = createPreconditioner<ST>(A, sb);
m_preconditioner = prec;
if (!prec.is_null()) {
// Trilinos BiCGStab does not support left preconditioners
if (sb.getSolverMethod() == escript::SO_METHOD_BICGSTAB)
problem->setRightPrec(prec);
else
problem->setLeftPrec(prec);
}
problem->setHermitian(sb.isSymmetric());
problem->setProblem();
} else {
for (auto t: problem->getTimers()) {
t->reset();
}
if (m_resetCalled) {
// special case for MueLu preconditioner - call Reuse...
// which honours the "reuse: type" parameter.
RCP<MueLu::TpetraOperator<ST,LO,GO,NT> > mlOp =
Teuchos::rcp_dynamic_cast<MueLu::TpetraOperator<ST,LO,GO,NT> >(m_preconditioner);
if (mlOp.get()) {
RCP<Matrix> A_(Teuchos::rcp_const_cast<Matrix>(A));
MueLu::ReuseTpetraPreconditioner(A_, *mlOp);
}
}
problem->setProblem(X, B);
}
double t1 = Teuchos::Time::wallTime();
RCP<SolverType<ST> > solver = createSolver<ST>(sb);
if (sb.isVerbose()) {
std::cout << "Using " << solver->description() << std::endl;
}
solver->setProblem(problem);
Belos::ReturnType result = solver->solve();
double t2 = Teuchos::Time::wallTime();
const int numIters = solver->getNumIters();
double tol = sb.getTolerance();
try {
tol = solver->achievedTol();
} catch (...) {
}
if (sb.isVerbose()) {
if (result == Belos::Converged) {
sb.updateDiagnostics("converged", true);
//.........这里部分代码省略.........
开发者ID:svn2github，项目名称:Escript，代码行数:101，代码来源:CrsMatrixWrapper.cpp
示例4: rcp
RCP<SimpleFunctionObject<OperandType> >
DivisionFunctionXMLConverter<OperandType>::getSpecificSimpleFunction(
OperandType operand) const
{
return rcp(new DivisionFunction<OperandType>(operand));
}
开发者ID:00liujj，项目名称:trilinos，代码行数:6，代码来源:Teuchos_StandardFunctionObjectXMLConverters.hpp
示例5: TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(point_values, md_field_evaluate)
{
typedef panzer::Traits::FadType ScalarType;
typedef PHX::MDField<ScalarType> ArrayType;
typedef PHX::KokkosViewFactory<ScalarType,PHX::Device> ViewFactory;
typedef PHX::MDField<double>::size_type size_type;
Teuchos::RCP<shards::CellTopology> topo =
Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
const int num_cells = 4;
const int base_cell_dimension = 2;
const panzer::CellData cell_data(num_cells,topo);
int num_points = 3;
RCP<PointRule> point_rule = rcp(new PointRule("RandomPoints",num_points, cell_data));
TEST_EQUALITY(point_rule->num_points,num_points);
panzer::PointValues<ScalarType,PHX::MDField<ScalarType> > point_values;
panzer::MDFieldArrayFactory af("prefix_");
point_values.setupArrays(point_rule,af);
// Set up node coordinates. Here we assume the following
// ordering. This needs to be consistent with shards topology,
// otherwise we will get negative determinates
// 3(0,1)---2(1,1)
// | 0 |
// | |
// 0(0,0)---1(1,0)
const size_type derivative_dim = 4;
const std::vector<PHX::index_size_type> ddims(1,derivative_dim);
const int num_vertices = point_rule->topology->getNodeCount();
ArrayType node_coordinates = af.buildArray<ScalarType,Cell,NODE,Dim>("node_coordinates",num_cells, num_vertices, base_cell_dimension);
node_coordinates.setFieldData(ViewFactory::buildView(node_coordinates.fieldTag(),ddims));
const size_type x = 0;
const size_type y = 1;
for (size_type cell = 0; cell < node_coordinates.dimension(0); ++cell) {
int xleft = cell % 2;
int yleft = int(cell/2);
node_coordinates(cell,0,x) = xleft*0.5;
node_coordinates(cell,0,y) = yleft*0.5;
node_coordinates(cell,1,x) = (xleft+1)*0.5;
node_coordinates(cell,1,y) = yleft*0.5;
node_coordinates(cell,2,x) = (xleft+1)*0.5;
node_coordinates(cell,2,y) = (yleft+1)*0.5;
node_coordinates(cell,3,x) = xleft*0.5;
node_coordinates(cell,3,y) = (yleft+1)*0.5;
out << "Cell " << cell << " = ";
for(int i=0;i<4;i++)
out << "(" << node_coordinates(cell,i,x) << ", "
<< node_coordinates(cell,i,y) << ") ";
out << std::endl;
}
// Build the evaluation points
ArrayType point_coordinates = af.buildArray<ScalarType,IP,Dim>("points",num_points, base_cell_dimension);
point_coordinates.setFieldData(ViewFactory::buildView(point_coordinates.fieldTag(),ddims));
point_coordinates(0,0) = 0.0; point_coordinates(0,1) = 0.0; // mid point
point_coordinates(1,0) = 0.5; point_coordinates(1,1) = 0.5; // mid point of upper left quadrant
point_coordinates(2,0) = -0.5; point_coordinates(2,1) = 0.0; // mid point of line from center to left side
point_values.coords_ref.setFieldData(ViewFactory::buildView(point_values.coords_ref.fieldTag(),ddims));
point_values.node_coordinates.setFieldData(ViewFactory::buildView(point_values.node_coordinates.fieldTag(),ddims));
point_values.point_coords.setFieldData(ViewFactory::buildView(point_values.point_coords.fieldTag(),ddims));
point_values.jac.setFieldData(ViewFactory::buildView(point_values.jac.fieldTag(),ddims));
point_values.jac_inv.setFieldData(ViewFactory::buildView(point_values.jac_inv.fieldTag(),ddims));
point_values.jac_det.setFieldData(ViewFactory::buildView(point_values.jac_det.fieldTag(),ddims));
point_values.evaluateValues(node_coordinates,point_coordinates);
// check the reference values (ensure copying)
for(int p=0;p<num_points;p++)
for(size_type d=0;d<base_cell_dimension;d++)
TEST_EQUALITY(point_values.coords_ref(p,d).val(),point_coordinates(p,d).val());
// check the shifted values (ensure physical mapping)
for(int c=0;c<num_cells;c++) {
double dx = 0.5;
double dy = 0.5;
for(int p=0;p<num_points;p++) {
double x = dx*(point_coordinates(p,0).val()+1.0)/2.0 + node_coordinates(c,0,0).val();
double y = dy*(point_coordinates(p,1).val()+1.0)/2.0 + node_coordinates(c,0,1).val();
TEST_FLOATING_EQUALITY(point_values.point_coords(c,p,0).val(),x,1e-10);
TEST_FLOATING_EQUALITY(point_values.point_coords(c,p,1).val(),y,1e-10);
}
}
//.........这里部分代码省略.........
开发者ID:uppatispr，项目名称:trilinos-official，代码行数:101，代码来源:point_values.cpp
示例6: m
void UncoupledAggregationFactory<LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Build(Level ¤tLevel) const {
FactoryMonitor m(*this, "Build", currentLevel);
const ParameterList& pL = GetParameterList();
bDefinitionPhase_ = false; // definition phase is finished, now all aggregation algorithm information is fixed
bool bUseOnePtAggregationAlgorithm = pL.get<bool>("UseOnePtAggregationAlgorithm");
bool bUseSmallAggregationAlgorithm = pL.get<bool>("UseSmallAggregatesAggregationAlgorithm");
bool bUsePreserveDirichletAggregationAlgorithm = pL.get<bool>("UsePreserveDirichletAggregationAlgorithm");
bool bUseUncoupledAggregationAglorithm = pL.get<bool>("UseUncoupledAggregationAlgorithm");
bool bUseMaxLinkAggregationAlgorithm = pL.get<bool>("UseMaxLinkAggregationAlgorithm");
bool bUseIsolatedNodeAggregationAglorithm = pL.get<bool>("UseIsolatedNodeAggregationAlgorithm");
bool bUseEmergencyAggregationAlgorithm = pL.get<bool>("UseEmergencyAggregationAlgorithm");
// define aggregation algorithms
RCP<const FactoryBase> graphFact = GetFactory("Graph");
// TODO Can we keep different aggregation algorithms over more Build calls?
algos_.clear();
if (bUseOnePtAggregationAlgorithm) algos_.push_back(rcp(new OnePtAggregationAlgorithm (graphFact)));
if (bUseSmallAggregationAlgorithm) algos_.push_back(rcp(new SmallAggregationAlgorithm (graphFact)));
if (bUseUncoupledAggregationAglorithm) algos_.push_back(rcp(new UncoupledAggregationAlgorithm (graphFact)));
if (bUseMaxLinkAggregationAlgorithm) algos_.push_back(rcp(new MaxLinkAggregationAlgorithm (graphFact)));
if (bUsePreserveDirichletAggregationAlgorithm) algos_.push_back(rcp(new PreserveDirichletAggregationAlgorithm (graphFact)));
if (bUseIsolatedNodeAggregationAglorithm) algos_.push_back(rcp(new IsolatedNodeAggregationAlgorithm (graphFact)));
if (bUseEmergencyAggregationAlgorithm) algos_.push_back(rcp(new EmergencyAggregationAlgorithm (graphFact)));
std::string mapOnePtName = pL.get<std::string>("OnePt aggregate map name"), mapSmallAggName = pL.get<std::string>("SmallAgg aggregate map name");
RCP<const Map> OnePtMap, SmallAggMap;
if (mapOnePtName.length()) {
RCP<const FactoryBase> mapOnePtFact = GetFactory("OnePt aggregate map factory");
OnePtMap = currentLevel.Get<RCP<const Map> >(mapOnePtName, mapOnePtFact.get());
}
if (mapSmallAggName.length()) {
RCP<const FactoryBase> mapSmallAggFact = GetFactory("SmallAgg aggregate map factory");
SmallAggMap = currentLevel.Get<RCP<const Map> >(mapSmallAggName, mapSmallAggFact.get());
}
RCP<const GraphBase> graph = Get< RCP<GraphBase> >(currentLevel, "Graph");
// Build
RCP<Aggregates> aggregates = rcp(new Aggregates(*graph));
aggregates->setObjectLabel("UC");
const LO nRows = graph->GetNodeNumVertices();
// construct aggStat information
std::vector<unsigned> aggStat(nRows, NodeStats::READY);
ArrayRCP<const bool> dirichletBoundaryMap = graph->GetBoundaryNodeMap();
if (dirichletBoundaryMap != Teuchos::null) {
for (LO i = 0; i < nRows; i++)
if (dirichletBoundaryMap[i] == true)
aggStat[i] = NodeStats::BOUNDARY;
}
LO nDofsPerNode = Get<LO>(currentLevel, "DofsPerNode");
GO indexBase = graph->GetDomainMap()->getIndexBase();
if (SmallAggMap != Teuchos::null || OnePtMap != Teuchos::null) {
for (LO i = 0; i < nRows; i++) {
// reconstruct global row id (FIXME only works for contiguous maps)
GO grid = (graph->GetDomainMap()->getGlobalElement(i)-indexBase) * nDofsPerNode + indexBase;
if (SmallAggMap != null) {
for (LO kr = 0; kr < nDofsPerNode; kr++) {
if (SmallAggMap->isNodeGlobalElement(grid + kr))
aggStat[i] = MueLu::NodeStats::SMALLAGG;
}
}
if (OnePtMap != null) {
for (LO kr = 0; kr < nDofsPerNode; kr++) {
if (OnePtMap->isNodeGlobalElement(grid + kr))
aggStat[i] = MueLu::NodeStats::ONEPT;
}
}
}
}
const RCP<const Teuchos::Comm<int> > comm = graph->GetComm();
GO numGlobalRows = 0;
if (IsPrint(Statistics1))
sumAll(comm, as<GO>(nRows), numGlobalRows);
LO numNonAggregatedNodes = nRows;
GO numGlobalAggregatedPrev = 0, numGlobalAggsPrev = 0;
for (size_t a = 0; a < algos_.size(); a++) {
std::string phase = algos_[a]->description();
SubFactoryMonitor sfm(*this, "Algo \"" + phase + "\"", currentLevel);
algos_[a]->BuildAggregates(pL, *graph, *aggregates, aggStat, numNonAggregatedNodes);
if (IsPrint(Statistics1)) {
GO numLocalAggregated = nRows - numNonAggregatedNodes, numGlobalAggregated = 0;
GO numLocalAggs = aggregates->GetNumAggregates(), numGlobalAggs = 0;
sumAll(comm, numLocalAggregated, numGlobalAggregated);
sumAll(comm, numLocalAggs, numGlobalAggs);
//.........这里部分代码省略.........
开发者ID:jgoldfar，项目名称:trilinos，代码行数:101，代码来源:MueLu_UncoupledAggregationFactory_def.hpp
示例7: TEUCHOS_UNIT_TEST
// quad tests
TEUCHOS_UNIT_TEST(tSquareQuadMeshDOFManager_edgetests, buildTest_quad_edge_orientations_fail)
{
// build global (or serial communicator)
#ifdef HAVE_MPI
stk_classic::ParallelMachine Comm = MPI_COMM_WORLD;
#else
stk_classic::ParallelMachine Comm = WHAT_TO_DO_COMM;
#endif
int numProcs = stk_classic::parallel_machine_size(Comm);
TEUCHOS_ASSERT(numProcs==1);
// build a geometric pattern from a single basis
RCP<const panzer::FieldPattern> patternI1
= buildFieldPattern<Intrepid2::Basis_HCURL_QUAD_I1_FEM<double,FieldContainer> >();
out << *patternI1 << std::endl;
RCP<panzer::ConnManager<int,int> > connManager = buildQuadMesh(Comm,2,2,1,1);
RCP<panzer::DOFManagerFEI<int,int> > dofManager = rcp(new panzer::DOFManagerFEI<int,int>());
dofManager->setOrientationsRequired(true);
TEST_EQUALITY(dofManager->getOrientationsRequired(),true);
TEST_EQUALITY(dofManager->getConnManager(),Teuchos::null);
dofManager->setConnManager(connManager,MPI_COMM_WORLD);
TEST_EQUALITY(dofManager->getConnManager(),connManager);
dofManager->addField("b",patternI1);
dofManager->buildGlobalUnknowns();
for(int i=0;i<4;i++) {
const int * indices = connManager->getConnectivity(i);
TEST_EQUALITY(connManager->getConnectivitySize(i),8);
out << "cell = " << i << ": ";
for(int j=0;j<4;j++)
out << indices[j+4] << " ";
out << std::endl;
}
out << "GIDS" << std::endl;
for(int i=0;i<4;i++) {
std::vector<int> gids;
dofManager->getElementGIDs(i,gids);
TEST_EQUALITY(gids.size(),4);
out << "cell = " << i << ": ";
for(int j=0;j<4;j++)
out << gids[j] << " ";
out << std::endl;
}
std::vector<int> total;
dofManager->getOwnedIndices(total);
TEST_EQUALITY(total.size(),12);
dofManager->printFieldInformation(out);
}
开发者ID:cihanuq，项目名称:Trilinos，代码行数:64，代码来源:tSquareQuadMeshDOFManager_edgetests_FEI.cpp
示例8: Setup
void Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetupSchwarz(Level& currentLevel) {
if (this->IsSetup() == true)
this->GetOStream(Warnings0) << "MueLu::Ifpack2Smoother::Setup(): Setup() has already been called" << std::endl;
// If we are doing "user" partitioning, we assume that what the user
// really wants to do is make tiny little subdomains with one row
// asssigned to each subdomain. The rows used for these little
// subdomains correspond to those in the 2nd block row. Then,
// if we overlap these mini-subdomains, we will do something that
// looks like Vanka (grabbing all velocities associated with each
// each pressure unknown). In addition, we put all Dirichlet points
// as a little mini-domain.
ParameterList& paramList = const_cast<ParameterList&>(this->GetParameterList());
bool isBlockedMatrix = false;
RCP<Matrix> merged2Mat;
std::string sublistName = "subdomain solver parameters";
if (paramList.isSublist(sublistName)) {
ParameterList& subList = paramList.sublist(sublistName);
std::string partName = "partitioner: type";
if (subList.isParameter(partName) && subList.get<std::string>(partName) == "user") {
isBlockedMatrix = true;
RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
"Matrix A must be of type BlockedCrsMatrix.");
size_t numVels = bA->getMatrix(0,0)->getNodeNumRows();
size_t numPres = bA->getMatrix(1,0)->getNodeNumRows();
size_t numRows = A_->getNodeNumRows();
ArrayRCP<LocalOrdinal> blockSeeds(numRows, Teuchos::OrdinalTraits<LocalOrdinal>::invalid());
size_t numBlocks = 0;
for (size_t rowOfB = numVels; rowOfB < numVels+numPres; ++rowOfB)
blockSeeds[rowOfB] = numBlocks++;
RCP<BlockedCrsMatrix> bA2 = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
TEUCHOS_TEST_FOR_EXCEPTION(bA2.is_null(), Exceptions::BadCast,
"Matrix A must be of type BlockedCrsMatrix.");
RCP<CrsMatrix> mergedMat = bA2->Merge();
merged2Mat = rcp(new CrsMatrixWrap(mergedMat));
// Add Dirichlet rows to the list of seeds
ArrayRCP<const bool> boundaryNodes;
boundaryNodes = Utilities::DetectDirichletRows(*merged2Mat, 0.0);
bool haveBoundary = false;
for (LO i = 0; i < boundaryNodes.size(); i++)
if (boundaryNodes[i]) {
// FIXME:
// 1. would not this [] overlap with some in the previos blockSeed loop?
// 2. do we need to distinguish between pressure and velocity Dirichlet b.c.
blockSeeds[i] = numBlocks;
haveBoundary = true;
}
if (haveBoundary)
numBlocks++;
subList.set("partitioner: map", blockSeeds);
subList.set("partitioner: local parts", as<int>(numBlocks));
}
}
RCP<const Tpetra::RowMatrix<SC, LO, GO, NO> > tpA;
if (isBlockedMatrix == true) tpA = Utilities::Op2NonConstTpetraRow(merged2Mat);
else tpA = Utilities::Op2NonConstTpetraRow(A_);
prec_ = Ifpack2::Factory::create(type_, tpA, overlap_);
SetPrecParameters();
prec_->initialize();
prec_->compute();
}
开发者ID:crtrott，项目名称:Trilinos，代码行数:75，代码来源:MueLu_Ifpack2Smoother_def.hpp
示例9: rcp
RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Ifpack2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
RCP<Ifpack2Smoother> smoother = rcp(new Ifpack2Smoother(*this) );
smoother->SetParameterList(this->GetParameterList());
return smoother;
}
开发者ID:crtrott，项目名称:Trilinos，代码行数:5，代码来源:MueLu_Ifpack2Smoother_def.hpp
示例10: insertRow
void
dft_PolyA22_Tpetra_Operator<Scalar,MatrixType>::
finalizeProblemValues
()
{
if (isLinearProblemSet_)
{
return; // nothing to do
}
insertRow(); // Dump any remaining entries
if (firstTime_) {
RCP<ParameterList> pl = rcp(new ParameterList(parameterList_->sublist("fillCompleteList")));
pl->set( "Preserve Local Graph", true );
cmsOnCmsMatrix_->fillComplete(pl);
ArrayRCP<size_t> numEntriesPerRow(cmsMap_->getNodeNumElements());
for (LocalOrdinal i = OTLO::zero(); i < cmsMap_->getNodeNumElements(); ++i) {
numEntriesPerRow[i] = cmsOnCmsMatrix_->getNumEntriesInLocalRow( i );
}
cmsOnCmsGraph_ = rcp(new GRAPH(cmsMap_, cmsOnCmsMatrix_->getColMap(), numEntriesPerRow, Tpetra::StaticProfile));
for (LocalOrdinal i = OTLO::zero(); i < cmsMap_->getNodeNumElements(); ++i) {
ArrayView<const GlobalOrdinal> indices;
ArrayView<const MatScalar> values;
cmsOnCmsMatrix_->getLocalRowView( i, indices, values );
cmsOnCmsGraph_->insertLocalIndices( i, indices );
}
cmsOnCmsGraph_->fillComplete();
cmsOnCmsMatrixStatic_ = rcp(new MAT(cmsOnCmsGraph_));
cmsOnCmsMatrixStatic_->setAllToScalar(STMS::zero());
for (LocalOrdinal i = OTLO::zero(); i < cmsMap_->getNodeNumElements(); ++i) {
ArrayView<const GlobalOrdinal> indices;
ArrayView<const MatScalar> values;
cmsOnCmsMatrix_->getLocalRowView( i, indices, values );
cmsOnCmsMatrixStatic_->sumIntoLocalValues( i, indices(), values() );
}
cmsOnCmsMatrixStatic_->fillComplete();
cmsOnCmsMatrixOp_ = rcp(new MMOP(cmsOnCmsMatrixStatic_));
}
if (!cmsOnCmsMatrixStatic_->isFillComplete()) {
RCP<ParameterList> pl = rcp(new ParameterList(parameterList_->sublist("fillCompleteList")));
cmsOnCmsMatrixStatic_->fillComplete(pl);
}
RCP<ParameterList> pl = rcp(new ParameterList(parameterList_->sublist("fillCompleteList")));
if (!isFLinear_) {
insertRow(); // Dump any remaining entries
cmsOnDensityMatrix_->fillComplete(densityMap_, cmsMap_, pl);
}
if (!hasDensityOnCms_) // Confirm that densityOnCmsMatrix is zero
{
// Scalar normvalue = densityOnCmsMatrix_->normInf();
// TEUCHOS_TEST_FOR_EXCEPT(normvalue!=0.0);
} else {
insertRow(); // Dump any remaining entries
densityOnCmsMatrix_->fillComplete(cmsMap_, densityMap_, pl);
}
// Form the inverse of the densityOnDensityMatrix
densityOnDensityInverse_->reciprocal(*densityOnDensityMatrix_);
// Use a diagonal preconditioner for the cmsOnCmsMatrix
if (firstTime_) {
RCP<const MAT> const_matrix = Teuchos::rcp_implicit_cast<const MAT>(cmsOnCmsMatrixStatic_);
cmsOnCmsInverse_ = rcp(new DIAGONAL(const_matrix));
cmsOnCmsInverseOp_ = rcp(new DIAGONAL_OP(cmsOnCmsInverse_));
#ifdef KDEBUG
TEUCHOS_TEST_FOR_EXCEPT(cmsOnCmsInverse_==Teuchos::null);
#endif
cmsOnCmsInverse_->initialize();
}
cmsOnCmsInverse_->compute();
// Compute the total number of entries in the A22 block
if (firstTime_) {
nnz_ = cmsOnCmsMatrixStatic_->getGlobalNumEntries() + \
cmsOnDensityMatrix_->getGlobalNumEntries() + \
densityOnCmsMatrix_->getGlobalNumEntries() + \
densityOnDensityMatrix_->getGlobalLength();
}
isLinearProblemSet_ = true;
firstTime_ = false;
} //end finalizeProblemValues
开发者ID:symashayak，项目名称:tramonto，代码行数:89，代码来源:dft_PolyA22_Tpetra_Operator.cpp
示例11: TEUCHOS_UNIT_TEST
// triangle tests
TEUCHOS_UNIT_TEST(tFieldPattern, test_equals)
{
out << note << std::endl;
RCP<Intrepid::Basis<double,FieldContainer> > basisA;
RCP<Intrepid::Basis<double,FieldContainer> > basisB;
basisA = rcp(new Intrepid::Basis_HGRAD_HEX_C1_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_HEX_C1_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_HEX_C1_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_HEX_C2_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(not intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_HEX_C2_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_HEX_C2_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_QUAD_C1_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_HEX_C1_FEM<double,FieldContainer>);
TEST_ASSERT(not intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(not intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_QUAD_C1_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_QUAD_C1_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_QUAD_C2_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_QUAD_C1_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(not intrepid_equals(basisA,basisB,__FILE__,__LINE__));
basisA = rcp(new Intrepid::Basis_HGRAD_QUAD_C2_FEM<double,FieldContainer>);
basisB = rcp(new Intrepid::Basis_HGRAD_QUAD_C2_FEM<double,FieldContainer>);
TEST_ASSERT(intrepid_same_geom(basisA,basisB,__FILE__,__LINE__));
TEST_ASSERT(intrepid_equals(basisA,basisB,__FILE__,__LINE__));
}
开发者ID:00liujj，项目名称:trilinos，代码行数:43，代码来源:tFieldPattern.cpp
示例12: getRcp
virtual RCP<ExprBase> getRcp() {return rcp(this);}
开发者ID:00liujj，项目名称:trilinos，代码行数:1，代码来源:SundanceDiscreteFunctionStub.hpp
示例13: pow_expand
RCP<const Basic> pow_expand(const RCP<const Pow> &self)
{
if (is_a<Integer>(*self->exp_)) {
if (is_a<Add>(*self->base_)) {
map_vec_mpz r;
int n = rcp_static_cast<const Integer>(self->exp_)->as_int();
RCP<const Add> base = rcp_static_cast<const Add>(self->base_);
umap_basic_int base_dict = base->dict_;
if (! (base->coef_->is_zero())) {
// Add the numerical coefficient into the dictionary. This
// allows a little bit easier treatment below.
insert(base_dict, base->coef_, one);
}
int m = base_dict.size();
multinomial_coefficients_mpz(m, n, r);
umap_basic_int rd;
// This speeds up overall expansion. For example for the benchmark
// (y + x + z + w)^60 it improves the timing from 135ms to 124ms.
rd.reserve(2*r.size());
RCP<const Number> add_overall_coeff=zero;
for (auto &p: r) {
auto power = p.first.begin();
auto i2 = base_dict.begin();
map_basic_basic d;
RCP<const Number> overall_coeff=one;
for (; power != p.first.end(); ++power, ++i2) {
if (*power > 0) {
RCP<const Integer> exp = rcp(new Integer(*power));
RCP<const Basic> base = i2->first;
if (is_a<Integer>(*base)) {
imulnum(outArg(overall_coeff),
rcp_static_cast<const Number>(
rcp_static_cast<const Integer>(base)->powint(*exp)));
} else if (is_a<Symbol>(*base)) {
Mul::dict_add_term(d, exp, base);
} else {
RCP<const Basic> exp2, t, tmp;
tmp = pow(base, exp);
Mul::as_base_exp(tmp, outArg(exp2), outArg(t));
Mul::dict_add_term(d, exp2, t);
}
if (!(i2->second->is_one())) {
imulnum(outArg(overall_coeff),
pownum(i2->second,
rcp_static_cast<const Number>(exp)));
}
}
}
RCP<const Basic> term = Mul::from_dict(overall_coeff, d);
RCP<const Number> coef2 = rcp(new Integer(p.second));
if (is_a_Number(*term)) {
iaddnum(outArg(add_overall_coeff),
mulnum(rcp_static_cast<const Number>(term), coef2));
} else {
if (is_a<Mul>(*term) &&
!(rcp_static_cast<const Mul>(term)->coef_->is_one())) {
// Tidy up things like {2x: 3} -> {x: 6}
imulnum(outArg(coef2),
rcp_static_cast<const Mul>(term)->coef_);
term = Mul::from_dict(one,
rcp_static_cast<const Mul>(term)->dict_);
}
Add::dict_add_term(rd, coef2, term);
}
}
RCP<const Basic> result = Add::from_dict(add_overall_coeff, rd);
return result;
}
}
return self;
}
开发者ID:curiousguy13，项目名称:csympy，代码行数:72，代码来源:pow.cpp
示例14: TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(point_values, intrepid_container_dfad)
{
PHX::KokkosDeviceSession session;
Teuchos::RCP<shards::CellTopology> topo =
Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));
const int num_cells = 4;
const int base_cell_dimension = 2;
const panzer::CellData cell_data(num_cells,topo);
int num_points = 3;
RCP<PointRule> point_rule = rcp(new PointRule("RandomPoints",num_points, cell_data));
TEST_EQUALITY(point_rule->num_points,num_points);
typedef panzer::Traits::FadType ScalarType;
panzer::PointValues<ScalarType,Intrepid2::FieldContainer<ScalarType> > point_values;
panzer::Intrepid2FieldContainerFactory af;
point_values.setupArrays(point_rule,af);
// Set up node coordinates. Here we assume the following
// ordering. This needs to be consistent with shards topology,
// otherwise we will get negative determinates
// 3(0,1)---2(1,1)
// | 0 |
// | |
// 0(0,0)---1(1,0)
const int num_vertices = point_rule->topology->getNodeCount();
Intrepid2::FieldContainer<ScalarType> node_coordinates(num_cells, num_vertices,
base_cell_dimension);
typedef panzer::ArrayTraits<ScalarType,FieldContainer<ScalarType> >::size_type size_type;
const size_type x = 0;
const size_type y = 1;
for (size_type cell = 0; cell < node_coordinates.dimension(0); ++cell) {
int xleft = cell % 2;
int yleft = int(cell/2);
node_coordinates(cell,0,x) = xleft*0.5;
node_coordinates(cell,0,y) = yleft*0.5;
node_coordinates(cell,1,x) = (xleft+1)*0.5;
node_coordinates(cell,1,y) = yleft*0.5;
node_coordinates(cell,2,x) = (xleft+1)*0.5;
node_coordinates(cell,2,y) = (yleft+1)*0.5;
node_coordinates(cell,3,x) = xleft*0.5;
node_coordinates(cell,3,y) = (yleft+1)*0.5;
out << "Cell " << cell << " = ";
for(int i=0; i<4; i++)
out << "(" << node_coordinates(cell,i,x) << ", "
<< node_coordinates(cell,i,y) << ") ";
out << std::endl;
}
// Build the evaluation points
Intrepid2::FieldContainer<ScalarType> point_coordinates(num_points, base_cell_dimension);
point_coordinates(0,0) = 0.0;
point_coordinates(0,1) = 0.0; // mid point
point_coordinates(1,0) = 0.5;
point_coordinates(1,1) = 0.5; // mid point of upper left quadrant
point_coordinates(2,0) = -0.5;
point_coordinates(2,1) = 0.0; // mid point of line from center to left side
point_values.evaluateValues(node_coordinates,point_coordinates);
}
开发者ID:crtrott，项目名称:Trilinos，代码行数:73，代码来源:point_values.cpp
示例15: getRcp
/** Return a ref count pointer to self */
virtual RCP<CurveBase> getRcp()
{
return rcp(this);
}
开发者ID:cakeisalie，项目名称:oomphlib_003，代码行数:5，代码来源:SundanceCircle.hpp
示例16: toXpetra
RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > toXpetra(const RCP<const Tpetra::Export<LocalOrdinal,GlobalOrdinal,Node> >& exp) {
if (!exp.is_null())
return rcp(new TpetraExport<LocalOrdinal,GlobalOrdinal,Node>(exp));
return Teuchos::null;
}
开发者ID:uppatispr，项目名称:trilinos-official，代码行数:6，代码来源:Xpetra_TpetraExport.hpp
示例17: TEUCHOS_UNIT_TEST
/*
* This test was created at the request of Chris Siefert to verify
* that some inexplicable behaviour in MueLu was not due to a faulty
* assumption in the Matrix Matrix Multiply Kernel.
* KLN 15/06/2011
*/
TEUCHOS_UNIT_TEST(Tpetra_MatMat, range_row_test){
RCP<const Comm<int> > comm = DefaultPlatform::getDefaultPlatform().getComm();
ParameterList defaultParameters;
RCP<node_type> node = rcp(new node_type(defaultParameters));
int numProcs = comm->getSize();
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//THIS NUMBER MUST BE EVEN SO THAT WHEN I CALCULATE THE NUMBER
//OF ROWS IN THE DOMAIN MAP I DON'T ENCOUNTER ANY
//WEIRD RESULTS DUE TO INTEGER DIVISION
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
int numRowsPerProc = 4;
int rank = comm->getRank();
global_size_t globalNumRows = numRowsPerProc*numProcs;
RCP<CrsMatrix<double,int,int,node_type> > identityMatrix =
getIdentityMatrix<double,int,int,int>(globalNumRows, comm, node);
//Create "B"
Array<int> myRows = tuple<int>(
rank*numRowsPerProc,
rank*numRowsPerProc+1,
rank*numRowsPerProc+2,
rank*numRowsPerProc+3);
Array<int> rangeElements;
if(rank == 0){
rangeElements = tuple<int>(
(numProcs-1)*numRowsPerProc+1,
(numProcs-1)*numRowsPerProc+2,
(numProcs-1)*numRowsPerProc,
(numProcs-1)*numRowsPerProc+3);
}
else{
rangeElements = tuple<int>(
(rank-1)*numRowsPerProc+1,
(rank-1)*numRowsPerProc+2,
(rank-1)*numRowsPerProc,
(rank-1)*numRowsPerProc+3);
}
RCP<const Map<int,int,node_type> > bRowMap =
Tpetra::createNonContigMapWithNode<int,int,node_type>(myRows, comm, node);
RCP<const Map<int,int,node_type> > bRangeMap =
Tpetra::createNonContigMapWithNode<int,int,node_type>(rangeElements, comm, node);
//We divide by 2 to make the matrix tall and "skinny"
RCP<const Map<int,int,node_type> > bDomainMap =
Tpetra::createUniformContigMapWithNode<int,int,node_type>(
globalNumRows/2, comm, node);
RCP<CrsMatrix<double,int,int,node_type> > bMatrix =
Tpetra::createCrsMatrix<double,int,int,node_type>(bRowMap, 1);
for(
ArrayView<const int>::iterator it =
bRowMap->getNodeElementList().begin();
it != bRowMap->getNodeElementList().end();
++it)
{
Array<int> col(1,(*it)/2);
Array<double> val(1,3.0);
bMatrix->insertGlobalValues(*it, col(), val());
}
bMatrix->fillComplete(bDomainMap, bRangeMap);
out << "Regular I*P" << std::endl;
mult_test_results results = multiply_test(
"Different Range and Row Maps",
identityMatrix,
bMatrix,
false,
false,
bMatrix,
comm,
out);
if(verbose){
out << "Results:" <<std::endl;
out << "\tEpsilon: " << results.epsilon << std::endl;
out << "\tcNorm: " << results.cNorm << std::endl;
out << "\tcompNorm: " << results.compNorm << std::endl;
}
TEST_COMPARE(results.epsilon, <, defaultEpsilon)
RCP<CrsMatrix<double,int,int,node_type> > identity2 =
getIdentityMatrix<double,int,int,int>(globalNumRows/2, comm, node);
RCP<const Map<int,int,node_type> > bTransRowMap =
Tpetra::createUniformContigMapWithNode<int,int,node_type>(globalNumRows/2,comm,node);
RCP<CrsMatrix<double,int,int,node_type> > bTrans =
Tpetra::createCrsMatrix<double,int,int,node_type>(bTransRowMap, 1);
Array<int> bTransRangeElements;
if(rank == 0){
bTransRangeElements = tuple<int>(
(numProcs-1)*(numRowsPerProc/2)+1,
(numProcs-1)*(numRowsPerProc/2));
}
//.........这里部分代码省略.........
开发者ID:crtrot |
客服电话
APP下载
官方微信
请发表评论