• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

C++ rcp函数代码示例

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

本文整理汇总了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 &currentLevel) 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

鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
C++ rcr0函数代码示例发布时间:2022-05-30
下一篇:
C++ rcc_wait_for_osc_ready函数代码示例发布时间:2022-05-30
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap