C++ ArrayRCP类代码示例

int main(int argc, char* argv[])
{

  using namespace std;
  using namespace Teuchos;

  const int num_samples = 1;
  const int num_loops = 500000;
  const int size = 10;
  const int num_vectors = 3;

  TEST_FOR_EXCEPTION(num_loops * size != 5000000, std::logic_error,
		     "Work amount is not constant!");

  // Make all vectors in a contiguous block

  MyVector<double>* vector_array = new MyVector<double>[num_vectors * size];
  ArrayRCP< MyVector<double> > a = 
    arcp< MyVector<double> >(vector_array, 0, size, false);
  ArrayRCP< MyVector<double> > b = 
    arcp< MyVector<double> >(&vector_array[size], 0, size, false);
  ArrayRCP< MyVector<double> > c = 
    arcp< MyVector<double> >(&vector_array[2*size], 0, size, false);

#ifdef HAVE_PHALANX_TVMET
  tvmet::Vector<double, 3>* tvmet_array = 
    new tvmet::Vector<double, 3>[num_vectors * size];
  ArrayRCP< tvmet::Vector<double, 3> > d = 
    arcp< tvmet::Vector<double, 3> >(tvmet_array, 0, size, false);
  ArrayRCP< tvmet::Vector<double, 3> > e = 
    arcp< tvmet::Vector<double, 3> >(&tvmet_array[size], 0, size, false);
  ArrayRCP< tvmet::Vector<double, 3> > f = 
    arcp< tvmet::Vector<double, 3> >(&tvmet_array[2*size], 0, size, false);
#endif

  double* raw_array = new double[num_vectors * size * 3];

  double* raw_a = raw_array;
  double* raw_b = &raw_array[size];
  double* raw_c = &raw_array[2*size];

  for (int i=0; i < a.size(); ++i)
    a[i] = 1.0;
  for (int i=0; i < b.size(); ++i)
    b[i] = 2.0;
  for (int i=0; i < c.size(); ++i)
    c[i] = 3.0;
#ifdef HAVE_PHALANX_TVMET
  for (int i=0; i < d.size(); ++i)
    d[i] = 1.0;
  for (int i=0; i < e.size(); ++i)
    e[i] = 2.0;
  for (int i=0; i < f.size(); ++i)
    f[i] = 3.0;
#endif
  for (int i=0; i < size; ++i) {
    int offset = i * 3;
    for (int j=0; j < 3; ++j) {
      raw_a[offset + j] = 1.0;
      raw_b[offset + j] = 2.0;
      raw_c[offset + j] = 3.0;
    }
  }

  RCP<Time> vector_time = TimeMonitor::getNewTimer("Vector Time");
  RCP<Time> update_time = TimeMonitor::getNewTimer("Update Time");
#ifdef HAVE_PHALANX_TVMET
  RCP<Time> tvmet_time = TimeMonitor::getNewTimer("TVMET Time");
#endif
  RCP<Time> raw_time = TimeMonitor::getNewTimer("Raw Time");
  RCP<Time> raw2_time = TimeMonitor::getNewTimer("Raw2 Time");

  for (int sample = 0; sample < num_samples; ++sample) {
    
    cout << "Vector" << endl;
    {
      TimeMonitor t(*vector_time);
      for (int i=0; i < num_loops; ++i)
	for (int j=0; j < c.size(); ++j)
	  c[j] = a[j] * b[j];
    } 
    
    cout << "Update" << endl;
    {
      TimeMonitor t(*update_time);
      for (int i=0; i < num_loops; ++i)
	for (int j=0; j < c.size(); ++j)
	  c[j].update_multiply(a[j], b[j]);
    }
    
#ifdef HAVE_PHALANX_TVMET
    cout << "TVMET" << endl;
    {
      TimeMonitor t(*tvmet_time);
      for (int i=0; i < num_loops; ++i)
	for (int j=0; j < d.size(); ++j)
	  f[j] = d[j] * e[j];
    }
#endif
    
//.........这里部分代码省略.........

  void RepartitionFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  DeterminePartitionPlacement(const Matrix& A, GOVector& decomposition, GO numPartitions) const {
    RCP<const Map> rowMap = A.getRowMap();

    RCP<const Teuchos::Comm<int> > comm = rowMap->getComm()->duplicate();
    int numProcs = comm->getSize();

    RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
    TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null, Exceptions::RuntimeError, "Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
    RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();

    const Teuchos::ParameterList& pL = GetParameterList();

    // maxLocal is a constant which determins the number of largest edges which are being exchanged
    // The idea is that we do not want to construct the full bipartite graph, but simply a subset of
    // it, which requires less communication. By selecting largest local edges we hope to achieve
    // similar results but at a lower cost.
    const int maxLocal = pL.get<int>("repartition: remap num values");
    const int dataSize = 2*maxLocal;

    ArrayRCP<GO> decompEntries;
    if (decomposition.getLocalLength() > 0)
      decompEntries = decomposition.getDataNonConst(0);

    // Step 1: Sort local edges by weight
    // Each edge of a bipartite graph corresponds to a triplet (i, j, v) where
    //   i: processor id that has some piece of part with part_id = j
    //   j: part id
    //   v: weight of the edge
    // We set edge weights to be the total number of nonzeros in rows on this processor which
    // correspond to this part_id. The idea is that when we redistribute matrix, this weight
    // is a good approximation of the amount of data to move.
    // We use two maps, original which maps a partition id of an edge to the corresponding weight,
    // and a reverse one, which is necessary to sort by edges.
    std::map<GO,GO> lEdges;
    for (LO i = 0; i < decompEntries.size(); i++)
      lEdges[decompEntries[i]] += A.getNumEntriesInLocalRow(i);

    // Reverse map, so that edges are sorted by weight.
    // This results in multimap, as we may have edges with the same weight
    std::multimap<GO,GO> revlEdges;
    for (typename std::map<GO,GO>::const_iterator it = lEdges.begin(); it != lEdges.end(); it++)
      revlEdges.insert(std::make_pair(it->second, it->first));

    // Both lData and gData are arrays of data which we communicate. The data is stored
    // in pairs, so that data[2*i+0] is the part index, and data[2*i+1] is the corresponding edge weight.
    // We do not store processor id in data, as we can compute that by looking on the offset in the gData.
    Array<GO> lData(dataSize, -1), gData(numProcs * dataSize);
    int numEdges = 0;
    for (typename std::multimap<GO,GO>::reverse_iterator rit = revlEdges.rbegin(); rit != revlEdges.rend() && numEdges < maxLocal; rit++) {
      lData[2*numEdges+0] = rit->second; // part id
      lData[2*numEdges+1] = rit->first;  // edge weight
      numEdges++;
    }

    // Step 2: Gather most edges
    // Each processors contributes maxLocal edges by providing maxLocal pairs <part id, weight>, which is of size dataSize
    MPI_Datatype MpiType = MpiTypeTraits<GO>::getType();
    MPI_Allgather(static_cast<void*>(lData.getRawPtr()), dataSize, MpiType, static_cast<void*>(gData.getRawPtr()), dataSize, MpiType, *rawMpiComm);

    // Step 3: Construct mapping

    // Construct the set of triplets
    std::vector<Triplet<int,int> > gEdges(numProcs * maxLocal);
    size_t k = 0;
    for (LO i = 0; i < gData.size(); i += 2) {
      GO part   = gData[i+0];
      GO weight = gData[i+1];
      if (part != -1) {                     // skip nonexistent edges
        gEdges[k].i = i/dataSize;           // determine the processor by its offset (since every processor sends the same amount)
        gEdges[k].j = part;
        gEdges[k].v = weight;
        k++;
      }
    }
    gEdges.resize(k);

    // Sort edges by weight
    // NOTE: compareTriplets is actually a reverse sort, so the edges weight is in decreasing order
    std::sort(gEdges.begin(), gEdges.end(), compareTriplets<int,int>);

    // Do matching
    std::map<int,int> match;
    std::vector<char> matchedRanks(numProcs, 0), matchedParts(numProcs, 0);
    int numMatched = 0;
    for (typename std::vector<Triplet<int,int> >::const_iterator it = gEdges.begin(); it != gEdges.end(); it++) {
      GO rank = it->i;
      GO part = it->j;
      if (matchedRanks[rank] == 0 && matchedParts[part] == 0) {
        matchedRanks[rank] = 1;
        matchedParts[part] = 1;
        match[part] = rank;
        numMatched++;
      }
    }
    GetOStream(Statistics0) << "Number of unassigned paritions before cleanup stage: " << (numPartitions - numMatched) << " / " << numPartitions << std::endl;

    // Step 4: Assign unassigned partitions
    // We do that through random matching for remaining partitions. Not all part numbers are valid, but valid parts are a subset of [0, numProcs).
    // The reason it is done this way is that we don't need any extra communication, as we don't need to know which parts are valid.
//.........这里部分代码省略.........

 /*! \brief Return the graph metric values.
  *  \param values on return is the array of values.
  */
 ArrayRCP<const GraphMetricValues<scalar_t> > getGraphMetrics() const{
     if(graphMetricsConst_.is_null()) return graphMetrics_;
     return graphMetricsConst_;
 }

void globalWeightedCutsMessagesHopsByPart(
    const RCP<const Environment> &env,
    const RCP<const Comm<int> > &comm,
    const RCP<const GraphModel<typename Adapter::base_adapter_t> > &graph,
    const ArrayView<const typename Adapter::part_t> &parts,
    typename Adapter::part_t &numParts,
    ArrayRCP<RCP<BaseClassMetrics<typename Adapter::scalar_t> > > &metrics,
    ArrayRCP<typename Adapter::scalar_t> &globalSums,
    const RCP <const MachineRep> machine)
{
  env->debug(DETAILED_STATUS, "Entering globalWeightedCutsMessagesHopsByPart");
  //////////////////////////////////////////////////////////
  // Initialize return values

  typedef typename Adapter::lno_t t_lno_t;
  typedef typename Adapter::gno_t t_gno_t;
  typedef typename Adapter::scalar_t t_scalar_t;
  typedef typename Adapter::part_t part_t;
  typedef typename Adapter::node_t t_node_t;


  typedef typename Zoltan2::GraphModel<typename Adapter::base_adapter_t>::input_t t_input_t;

  t_lno_t localNumVertices = graph->getLocalNumVertices();
  t_gno_t globalNumVertices = graph->getGlobalNumVertices();
  t_lno_t localNumEdges = graph->getLocalNumEdges();

  ArrayView<const t_gno_t> Ids;
  ArrayView<t_input_t> v_wghts;
  graph->getVertexList(Ids, v_wghts);

  typedef GraphMetrics<t_scalar_t> mv_t;

  //get the edge ids, and weights
  ArrayView<const t_gno_t> edgeIds;
  ArrayView<const t_lno_t> offsets;
  ArrayView<t_input_t> e_wgts;
  graph->getEdgeList(edgeIds, offsets, e_wgts);


  std::vector <t_scalar_t> edge_weights;
  int numWeightPerEdge = graph->getNumWeightsPerEdge();

  int numMetrics = 4;                   // "edge cuts", messages, hops, weighted hops
  if (numWeightPerEdge) numMetrics += numWeightPerEdge * 2;   // "weight n", weighted hops per weight n

  // add some more metrics to the array
  typedef typename ArrayRCP<RCP<BaseClassMetrics<typename Adapter::scalar_t> > >::size_type array_size_type;
  metrics.resize( metrics.size() + numMetrics );

  for( array_size_type n = metrics.size() - numMetrics; n < metrics.size(); ++n ){
    mv_t * newMetric = new mv_t;                  // allocate the new memory
    env->localMemoryAssertion(__FILE__,__LINE__,1,newMetric);   // check errors
    metrics[n] = rcp( newMetric);         // create the new members
  }
  array_size_type next = metrics.size() - numMetrics; // MDM - this is most likely temporary to preserve the format here - we are now filling a larger array so we may not have started at 0

  std::vector <part_t> e_parts (localNumEdges);
#ifdef HAVE_ZOLTAN2_MPI
  if (comm->getSize() > 1)
  {
    Zoltan_DD_Struct *dd = NULL;

    MPI_Comm mpicomm = Teuchos::getRawMpiComm(*comm);
    int size_gnot = Zoltan2::TPL_Traits<ZOLTAN_ID_PTR, t_gno_t>::NUM_ID;

    int debug_level = 0;
    Zoltan_DD_Create(&dd, mpicomm,
        size_gnot, 0,
        sizeof(part_t), localNumVertices, debug_level);

    ZOLTAN_ID_PTR ddnotneeded = NULL;  // Local IDs not needed
    Zoltan_DD_Update(
        dd,
        (ZOLTAN_ID_PTR) Ids.getRawPtr(),
        ddnotneeded,
        (char *) &(parts[0]),
        NULL,
        int(localNumVertices));

    Zoltan_DD_Find(
        dd,
        (ZOLTAN_ID_PTR) edgeIds.getRawPtr(),
        ddnotneeded,
        (char *)&(e_parts[0]),
        NULL,
        localNumEdges,
        NULL
        );
    Zoltan_DD_Destroy(&dd);
  } else
#endif
  {

    std::map<t_gno_t,t_lno_t> global_id_to_local_index;

    //else everything is local.
    //we need a globalid to local index conversion.
    //this does not exists till this point, so we need to create one.
    for (t_lno_t i = 0; i < localNumVertices; ++i){
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession session(&argc, &argv);
  RCP<const Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  int rank = comm->getRank();

  Teuchos::RCP<Teuchos::FancyOStream> outStream = 
    Teuchos::VerboseObjectBase::getDefaultOStream();
  Teuchos::EVerbosityLevel v=Teuchos::VERB_EXTREME;

  typedef Tpetra::CrsMatrix<zscalar_t,zlno_t,zgno_t,znode_t> tmatrix_t;
  typedef Tpetra::CrsGraph<zlno_t,zgno_t,znode_t> tgraph_t;
  typedef Tpetra::Vector<zscalar_t,zlno_t,zgno_t,znode_t> tvector_t;
  typedef Tpetra::MultiVector<zscalar_t,zlno_t,zgno_t,znode_t> tmvector_t;
  typedef Xpetra::CrsMatrix<zscalar_t,zlno_t,zgno_t,znode_t> xmatrix_t;
  typedef Xpetra::CrsGraph<zlno_t,zgno_t,znode_t> xgraph_t;
  typedef Xpetra::Vector<zscalar_t,zlno_t,zgno_t,znode_t> xvector_t;
  typedef Xpetra::MultiVector<zscalar_t,zlno_t,zgno_t,znode_t> xmvector_t;
  typedef Xpetra::TpetraMap<zlno_t,zgno_t,znode_t> xtmap_t;

  // Create object that can give us test Tpetra and Xpetra input.

  RCP<UserInputForTests> uinput;

  try{
    uinput = 
      rcp(new UserInputForTests(testDataFilePath,std::string("simple"), comm, true));
  }
  catch(std::exception &e){
    TEST_FAIL_AND_EXIT(*comm, 0, string("input ")+e.what(), 1);
  }

  /////////////////////////////////////////////////////////////////
  //   Tpetra::CrsMatrix
  //   Tpetra::CrsGraph
  //   Tpetra::Vector
  //   Tpetra::MultiVector
  /////////////////////////////////////////////////////////////////

  // XpetraTraits<Tpetra::CrsMatrix<zscalar_t, zlno_t, zgno_t, znode_t> > 
  {
    RCP<tmatrix_t> M;
  
    try{
      M = uinput->getUITpetraCrsMatrix();
    }
    catch(std::exception &e){
      TEST_FAIL_AND_EXIT(*comm, 0, 
        string("getTpetraCrsMatrix ")+e.what(), 1);
    }
  
    if (rank== 0)
      std::cout << "Original Tpetra matrix " << M->getGlobalNumRows()
        << " x " << M->getGlobalNumCols() << std::endl;

    M->describe(*outStream,v);

    RCP<const xtmap_t> xmap(new xtmap_t(M->getRowMap()));

    ArrayRCP<zgno_t> newRowIds = roundRobinMap(xmap);
  
    zgno_t localNumRows = newRowIds.size();
  
    RCP<const tmatrix_t> newM;
    try{
      newM = Zoltan2::XpetraTraits<tmatrix_t>::doMigration(*M,
        localNumRows, newRowIds.getRawPtr());
    }
    catch(std::exception &e){
      TEST_FAIL_AND_EXIT(*comm, 0, 
        string(" Zoltan2::XpetraTraits<tmatrix_t>::doMigration ")+e.what(), 1);
    }

    if (rank== 0)
      std::cout << "Migrated Tpetra matrix" << std::endl;
  
    newM->describe(*outStream,v);
  }

  // XpetraTraits<Tpetra::CrsGraph<zscalar_t, zlno_t, zgno_t, znode_t> > 
  {
    RCP<tgraph_t> G;
  
    try{
      G = uinput->getUITpetraCrsGraph();
    }
    catch(std::exception &e){
      TEST_FAIL_AND_EXIT(*comm, 0, 
        string("getTpetraCrsGraph ")+e.what(), 1);
    }
  
    if (rank== 0)
      std::cout << "Original Tpetra graph" << std::endl;
  
    G->describe(*outStream,v);
  
    RCP<const xtmap_t> xmap(new xtmap_t(G->getRowMap()));
    ArrayRCP<zgno_t> newRowIds = roundRobinMap(xmap);
  
    zgno_t localNumRows = newRowIds.size();
//.........这里部分代码省略.........

  void CoordinatesTransferFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Build(Level & fineLevel, Level &coarseLevel) const {
    FactoryMonitor m(*this, "Build", coarseLevel);

    GetOStream(Runtime0, 0) << "Transferring coordinates" << std::endl;

    const ParameterList  & pL = GetParameterList();
    int                 writeStart = pL.get< int >("write start");
    int                 writeEnd   = pL.get< int >("write end");

    RCP<Aggregates>     aggregates = Get< RCP<Aggregates> > (fineLevel, "Aggregates");
    RCP<MultiVector>    fineCoords = Get< RCP<MultiVector> >(fineLevel, "Coordinates");
    RCP<const Map>      coarseMap  = Get< RCP<const Map> >  (fineLevel, "CoarseMap");

    // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
    // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
    // logical blocks in the matrix

    ArrayView<const GO> elementAList = coarseMap->getNodeElementList();
    LO                  blkSize      = 1;
    if (rcp_dynamic_cast<const StridedMap>(coarseMap) != Teuchos::null)
      blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap)->getFixedBlockSize();

    GO                  indexBase    = coarseMap->getIndexBase();
    size_t              numElements  = elementAList.size() / blkSize;
    Array<GO>           elementList(numElements);

    // Amalgamate the map
    for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
      elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;

    RCP<const Map> coarseCoordMap = MapFactory        ::Build(coarseMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, coarseMap->getComm());
    RCP<MultiVector> coarseCoords = MultiVectorFactory::Build(coarseCoordMap, fineCoords->getNumVectors());

    // Maps
    RCP<const Map> uniqueMap    = fineCoords->getMap();
    RCP<const Map> nonUniqueMap = aggregates->GetMap();

    // Create overlapped fine coordinates to reduce global communication
    RCP<const Import>     importer = ImportFactory     ::Build(uniqueMap, nonUniqueMap);
    RCP<MultiVector> ghostedCoords = MultiVectorFactory::Build(nonUniqueMap, fineCoords->getNumVectors());
    ghostedCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);

    // Get some info about aggregates
    int                         myPID        = uniqueMap->getComm()->getRank();
    LO                          numAggs      = aggregates->GetNumAggregates();
    ArrayRCP<LO>                aggSizes     = aggregates->ComputeAggregateSizes();
    const ArrayRCP<const LO>    vertex2AggID = aggregates->GetVertex2AggId()->getData(0);
    const ArrayRCP<const LO>    procWinner   = aggregates->GetProcWinner()->getData(0);

    // Fill in coarse coordinates
    for (size_t j = 0; j < fineCoords->getNumVectors(); j++) {
      ArrayRCP<const Scalar> fineCoordsData = ghostedCoords->getData(j);
      ArrayRCP<Scalar>     coarseCoordsData = coarseCoords->getDataNonConst(j);

      for (LO lnode = 0; lnode < vertex2AggID.size(); lnode++)
        if (procWinner[lnode] == myPID)
          coarseCoordsData[vertex2AggID[lnode]] += fineCoordsData[lnode];

      for (LO agg = 0; agg < numAggs; agg++)
        coarseCoordsData[agg] /= aggSizes[agg];
    }

    Set<RCP<MultiVector> >(coarseLevel, "Coordinates", coarseCoords);
    if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd) {
      std::ostringstream buf;
      buf << fineLevel.GetLevelID();
      std::string fileName = "coordinates_before_rebalance_level_" + buf.str() + ".m";
      Utils::Write(fileName,*fineCoords);
    }
    if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd) {
      std::ostringstream buf;
      buf << coarseLevel.GetLevelID();
      std::string fileName = "coordinates_before_rebalance_level_" + buf.str() + ".m";
      Utils::Write(fileName,*coarseCoords);
    }

  } // Build

  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();
  }

size_t computeLocalEdgeList(
  const RCP<const Environment> &env, const RCP<const Comm<int> > &comm,
  size_t numLocalEdges,           // local edges
  size_t numLocalGraphEdges,      // edges in "local" graph
  RCP<const IdentifierMap<User> > &idMap,
  ArrayRCP<const typename InputTraits<User>::zgid_t> &allEdgeIds, // in
  ArrayRCP<const typename InputTraits<User>::gno_t> &allEdgeGnos, // in
  ArrayRCP<int> &allProcs,                                 // in
  ArrayRCP<const typename InputTraits<User>::lno_t> &allOffs,    // in
  ArrayRCP<StridedData<typename InputTraits<User>::lno_t,
                       typename InputTraits<User>::scalar_t> > &allWeights,// in
  ArrayRCP<const typename InputTraits<User>::lno_t> &edgeLocalIds, //
  ArrayRCP<const typename InputTraits<User>::lno_t> &offsets,      // out
  ArrayRCP<StridedData<typename InputTraits<User>::lno_t,
    typename InputTraits<User>::scalar_t> > &eWeights)             // out
{
  typedef typename InputTraits<User>::zgid_t zgid_t;
  typedef typename InputTraits<User>::gno_t gno_t;
  typedef typename InputTraits<User>::scalar_t scalar_t;
  typedef typename InputTraits<User>::lno_t lno_t;
  typedef StridedData<lno_t, scalar_t> input_t;
  int rank = comm->getRank();

  bool gnosAreGids = idMap->gnosAreGids();

  edgeLocalIds = ArrayRCP<const lno_t>(Teuchos::null);
  eWeights = ArrayRCP<input_t>(Teuchos::null);
  offsets = ArrayRCP<const lno_t>(Teuchos::null);

  if (numLocalGraphEdges == 0) {
    // Set the offsets array and return
    size_t allOffsSize = allOffs.size();
    lno_t *offs = new lno_t [allOffsSize];
    env->localMemoryAssertion(__FILE__, __LINE__, allOffsSize, offs);
    for (size_t i = 0; i < allOffsSize; i++) offs[i] = 0;
    offsets = arcp(offs, 0, allOffsSize, true);
    return 0;
  }

  if (numLocalGraphEdges == numLocalEdges){

    // Entire graph is local.

    lno_t *lnos = new lno_t [numLocalGraphEdges];
    env->localMemoryAssertion(__FILE__, __LINE__, numLocalGraphEdges, lnos);
    if (comm->getSize() == 1) {
      // With one rank, Can use gnos as local index.
      if (gnosAreGids)
        for (size_t i=0; i < numLocalEdges; i++) lnos[i] = allEdgeIds[i];
      else
        for (size_t i=0; i < numLocalEdges; i++) lnos[i] = allEdgeGnos[i];
    }
    else {
      ArrayRCP<gno_t> gnoArray;

      if (gnosAreGids){
        ArrayRCP<const gno_t> gnosConst =
                 arcp_reinterpret_cast<const gno_t>(allEdgeIds);
        gnoArray = arcp_const_cast<gno_t>(gnosConst);
      }
      else {
        gnoArray = arcp_const_cast<gno_t>(allEdgeGnos);
      }

      // Need to translate to gnos to local indexing
      ArrayView<lno_t> lnoView(lnos, numLocalGraphEdges);
      try {
        idMap->lnoTranslate(lnoView,
                            gnoArray.view(0,numLocalGraphEdges),
                            TRANSLATE_LIB_TO_APP);
      }
      Z2_FORWARD_EXCEPTIONS;
    }
    edgeLocalIds = arcp(lnos, 0, numLocalGraphEdges, true);
    offsets = allOffs;
    eWeights = allWeights;

  }

  void ZoltanInterface<LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::Build(Level& level) const {
    FactoryMonitor m(*this, "Build", level);

    RCP<Matrix>      A        = Get< RCP<Matrix> >     (level, "A");
    RCP<const Map>   rowMap   = A->getRowMap();

    RCP<MultiVector> Coords   = Get< RCP<MultiVector> >(level, "Coordinates");
    size_t           dim      = Coords->getNumVectors();

    GO               numParts = level.Get<GO>("number of partitions");

    if (numParts == 1) {
      // Running on one processor, so decomposition is the trivial one, all zeros.
      RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, true);
      Set(level, "Partition", decomposition);
      return;
    }

    float zoltanVersion_;
    Zoltan_Initialize(0, NULL, &zoltanVersion_);

    RCP<const Teuchos::MpiComm<int> >            dupMpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(rowMap->getComm()->duplicate());
    RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > zoltanComm = dupMpiComm->getRawMpiComm();

    RCP<Zoltan> zoltanObj_ = rcp(new Zoltan((*zoltanComm)()));  //extract the underlying MPI_Comm handle and create a Zoltan object
    if (zoltanObj_ == Teuchos::null)
      throw Exceptions::RuntimeError("MueLu::Zoltan : Unable to create Zoltan data structure");

    // Tell Zoltan what kind of local/global IDs we will use.
    // In our case, each GID is two ints and there are no local ids.
    // One can skip this step if the IDs are just single ints.
    int rv;
    if ((rv = zoltanObj_->Set_Param("num_gid_entries", "1")) != ZOLTAN_OK)
      throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'num_gid_entries' returned error code " + Teuchos::toString(rv));
    if ((rv = zoltanObj_->Set_Param("num_lid_entries", "0") ) != ZOLTAN_OK)
      throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'num_lid_entries' returned error code " + Teuchos::toString(rv));
    if ((rv = zoltanObj_->Set_Param("obj_weight_dim", "1") ) != ZOLTAN_OK)
      throw Exceptions::RuntimeError("MueLu::Zoltan::Setup : setting parameter 'obj_weight_dim' returned error code "  + Teuchos::toString(rv));

    if (GetVerbLevel() & Statistics1) zoltanObj_->Set_Param("debug_level", "1");
    else                              zoltanObj_->Set_Param("debug_level", "0");

    zoltanObj_->Set_Param("num_global_partitions", toString(numParts));

    zoltanObj_->Set_Num_Obj_Fn(GetLocalNumberOfRows,      (void *) &*A);
    zoltanObj_->Set_Obj_List_Fn(GetLocalNumberOfNonzeros, (void *) &*A);
    zoltanObj_->Set_Num_Geom_Fn(GetProblemDimension,      (void *) &dim);
    zoltanObj_->Set_Geom_Multi_Fn(GetProblemGeometry,     (void *) Coords.get());

    // Data pointers that Zoltan requires.
    ZOLTAN_ID_PTR import_gids = NULL;  // Global nums of objs to be imported
    ZOLTAN_ID_PTR import_lids = NULL;  // Local indices to objs to be imported
    int   *import_procs       = NULL;  // Proc IDs of procs owning objs to be imported.
    int   *import_to_part     = NULL;  // Partition #s to which imported objs should be assigned.
    ZOLTAN_ID_PTR export_gids = NULL;  // Global nums of objs to be exported
    ZOLTAN_ID_PTR export_lids = NULL;  // local indices to objs to be exported
    int   *export_procs       = NULL;  // Proc IDs of destination procs for objs to be exported.
    int   *export_to_part     = NULL;  // Partition #s for objs to be exported.
    int   num_imported;                // Number of objs to be imported.
    int   num_exported;                // Number of objs to be exported.
    int   newDecomp;                   // Flag indicating whether the decomposition has changed
    int   num_gid_entries;             // Number of array entries in a global ID.
    int   num_lid_entries;

    {
      SubFactoryMonitor m1(*this, "Zoltan RCB", level);
      rv = zoltanObj_->LB_Partition(newDecomp, num_gid_entries, num_lid_entries,
                                    num_imported, import_gids, import_lids, import_procs, import_to_part,
                                    num_exported, export_gids, export_lids, export_procs, export_to_part);
      if (rv == ZOLTAN_FATAL)
        throw Exceptions::RuntimeError("Zoltan::LB_Partition() returned error code");
    }

    // TODO check that A's row map is 1-1.  Zoltan requires this.

    RCP<Xpetra::Vector<GO, LO, GO, NO> > decomposition;
    if (newDecomp) {
      decomposition = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(rowMap, false); // Don't initialize, will be overwritten
      ArrayRCP<GO> decompEntries = decomposition->getDataNonConst(0);

      int mypid = rowMap->getComm()->getRank();
      for (typename ArrayRCP<GO>::iterator i = decompEntries.begin(); i != decompEntries.end(); ++i)
        *i = mypid;

      LO blockSize = A->GetFixedBlockSize();
      for (int i = 0; i < num_exported; ++i) {
        // We have assigned Zoltan gids to first row GID in the block
        // NOTE: Zoltan GIDs are different from GIDs in the Coordinates vector
        LO  localEl = rowMap->getLocalElement(export_gids[i]);
        int partNum = export_to_part[i];
        for (LO j = 0; j < blockSize; ++j)
          decompEntries[localEl + j] = partNum;
      }
    }

    Set(level, "Partition", decomposition);

    zoltanObj_->LB_Free_Part(&import_gids, &import_lids, &import_procs, &import_to_part);
    zoltanObj_->LB_Free_Part(&export_gids, &export_lids, &export_procs, &export_to_part);

//.........这里部分代码省略.........

// ********************************************************
int main(int argc, char *argv[]) 
{
  using namespace std;
  using namespace Teuchos;
  using namespace PHX;
  
  GlobalMPISession mpi_session(&argc, &argv);

  try {
    
    RCP<Time> total_time = TimeMonitor::getNewTimer("Total Run Time");
    TimeMonitor tm(*total_time);

    // *********************************************************************
    // Start of MDField Testing
    // *********************************************************************
    {

      typedef MDField<double,Cell,Node>::size_type size_type;

      std::vector<size_type> dims(3);
      dims[0] = 10;
      dims[1] = 4;
      dims[2] = 3;

      RCP<DataLayout> quad_vector = 
	rcp(new MDALayout<Cell,Quadrature,Dim>(dims[0],dims[1],dims[2]));
      
      int size = quad_vector->size();

      TEUCHOS_TEST_FOR_EXCEPTION(size != dims[0]*dims[1]*dims[2], std::runtime_error, 
			 "Size mismatch on MDField!");

      ArrayRCP<double> a_mem = arcp<double>(size);
      ArrayRCP<double> b_mem = arcp<double>(size);

      for (int i=0; i < a_mem.size(); ++i)
	a_mem[i] = static_cast<double>(i);

      for (int i=0; i < b_mem.size(); ++i)
	b_mem[i] = static_cast<double>(i);

      MDField<double,Cell,Point,Dim> a("density",quad_vector);
      MDField<double> b("density",quad_vector);

      a.setFieldData(a_mem);
      b.setFieldData(b_mem);

      simulated_intrepid_integrate(a);     
      simulated_intrepid_integrate(b);     

      // ***********************
      // Shards tests
      // ***********************

      ArrayRCP<double> c_mem = arcp<double>(size);
      ArrayRCP<double> d_mem = arcp<double>(size);

      for (int i=0; i < c_mem.size(); ++i)
	c_mem[i] = static_cast<double>(i);

      for (int i=0; i < d_mem.size(); ++i)
	d_mem[i] = static_cast<double>(i);

      shards::Array<double,shards::NaturalOrder,Cell,Node,Dim> c(c_mem.get(),
								 dims[0], 
								 dims[1],
								 dims[2]);

      size_type rank = dims.size();

      const ArrayRCP<const shards::ArrayDimTag*> tags = 
	arcp<const shards::ArrayDimTag*>(rank);
      tags[0] = &Cell::tag();
      tags[1] = &Point::tag();
      tags[2] = &Dim::tag();
      
      shards::Array<double,shards::NaturalOrder> d(d_mem.get(),rank,
						   &dims[0],tags.get());
      
      simulated_intrepid_integrate(d); 
      simulated_intrepid_integrate((const shards::Array<double,shards::NaturalOrder>&)(c));    

    }

    // *********************************************************************
    // *********************************************************************
    std::cout << "\nTest passed!\n" << std::endl; 
    // *********************************************************************
    // *********************************************************************

  }
  catch (const std::exception& e) {
    std::cout << "************************************************" << endl;
    std::cout << "************************************************" << endl;
    std::cout << "Exception Caught!" << endl;
    std::cout << "Error message is below\n " << e.what() << endl;
    std::cout << "************************************************" << endl;
  }
//.........这里部分代码省略.........

  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
    FactoryMonitor m(*this, "Setup Smoother", currentLevel);

    if (SmootherPrototype::IsSetup() == true)
      this->GetOStream(Warnings0) << "MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";

    // Extract blocked operator A from current level
    A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
    RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
    TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
                               "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");

    // Store map extractors
    rangeMapExtractor_  = bA->getRangeMapExtractor();
    domainMapExtractor_ = bA->getDomainMapExtractor();

    // Store the blocks in local member variables
    A00_ = bA->getMatrix(0,0);
    A01_ = bA->getMatrix(0,1);
    A10_ = bA->getMatrix(1,0);
    A11_ = bA->getMatrix(1,1);

    const ParameterList& pL = Factory::GetParameterList();
    SC omega = pL.get<SC>("Damping factor");

#if 0 // old code
    // Create the inverse of the diagonal of F
    D_ = VectorFactory::Build(A00_->getRowMap());

    ArrayRCP<SC> diag;
    if (pL.get<bool>("lumping") == false)
      diag = Utilities::GetMatrixDiagonal      (*A00_);
    else
      diag = Utilities::GetLumpedMatrixDiagonal(*A00_);

    SC one = Teuchos::ScalarTraits<SC>::one();

    ArrayRCP<SC> Ddata = D_->getDataNonConst(0);
    for (GO row = 0; row < Ddata.size(); row++)
      Ddata[row] = one / (diag[row]*omega);*/
#else
    // Create the inverse of the diagonal of F
    // TODO add safety check for zeros on diagonal of F!
    RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
    if (pL.get<bool>("lumping") == false) {
      A00_->getLocalDiagCopy(*diagFVector);       // extract diagonal of F
    } else {
      diagFVector = Utilities::GetLumpedMatrixDiagonal(A00_);
    }
    diagFVector->scale(omega);
    D_ = Utilities::GetInverse(diagFVector);
#endif

    // Set the Smoother
    // carefully switch to the SubFactoryManagers (defined by the users)
    {
      SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
      smoo_ = currentLevel.Get<RCP<SmootherBase> >("PreSmoother", FactManager_->GetFactory("Smoother").get());
      S_    = currentLevel.Get<RCP<Matrix> >      ("A",           FactManager_->GetFactory("A").get());
    }

    SmootherPrototype::IsSetup(true);
  }

//.........这里部分代码省略.........
    std::vector<RCP<const Xpetra::Map<LO,GO,Node> > > xmaps;
    xmaps.push_back(xstridedvelmap);
    xmaps.push_back(xstridedpremap);

    RCP<const Xpetra::MapExtractor<Scalar,LO,GO,Node> > map_extractor = Xpetra::MapExtractorFactory<Scalar,LO,GO,Node>::Build(xstridedfullmap,xmaps);

    /////////////////////////////////////// build blocked transfer operator
    // using the map extractor
    RCP<Xpetra::BlockedCrsMatrix<Scalar,LO,GO,Node> > bOp = rcp(new Xpetra::BlockedCrsMatrix<Scalar,LO,GO,Node>(map_extractor,map_extractor,10));
    bOp->setMatrix(0,0,xA11);
    bOp->setMatrix(0,1,xA12);
    bOp->setMatrix(1,0,xA21);
    bOp->setMatrix(1,1,xA22);

    bOp->fillComplete();

    //////////////////////////////////////// prepare setup
    ParameterListInterpreter mueLuFactory(xmlFile, *comm);


    RCP<Hierarchy> H = mueLuFactory.CreateHierarchy();
    H->setDefaultVerbLevel(VERB_HIGH);
    H->SetMaxCoarseSize(maxCoarseSize);

    RCP<MueLu::Level> Finest = H->GetLevel(0);
    Finest->setDefaultVerbLevel(VERB_HIGH);
    Finest->Set("A",           rcp_dynamic_cast<Matrix>(bOp));


    ////////////////////////////////////////// prepare null space for A11
    RCP<MultiVector> nullspace11 = MultiVectorFactory::Build(xstridedvelmap, 2);  // this is a 2D standard null space

    for (int i=0; i<nDofsPerNode-1; ++i) {
      ArrayRCP<Scalar> nsValues = nullspace11->getDataNonConst(i);
      int numBlocks = nsValues.size() / (nDofsPerNode - 1);
      for (int j=0; j< numBlocks; ++j) {
        nsValues[j*(nDofsPerNode - 1) + i] = 1.0;
      }
    }

    Finest->Set("Nullspace1",nullspace11);

    ////////////////////////////////////////// prepare null space for A22
    RCP<MultiVector> nullspace22 = MultiVectorFactory::Build(xstridedpremap, 1);  // this is a 2D standard null space
    ArrayRCP<Scalar> nsValues22 = nullspace22->getDataNonConst(0);
    for (int j=0; j< nsValues22.size(); ++j) {
      nsValues22[j] = 1.0;
    }

    Finest->Set("Nullspace2",nullspace22);

    /////////////////////////////////// BEGIN setup

    mueLuFactory.SetupHierarchy(*H);

    ///////////////////////////////////// END setup

    *out << std::endl;

    RCP<MultiVector> xLsg = MultiVectorFactory::Build(xstridedfullmap,1);

    // Use AMG directly as an iterative method
    {
      xLsg->putScalar( (SC) 0.0);

      // Epetra_Vector -> Xpetra::Vector

 EpetraCrsMatrixT<EpetraGlobalOrdinal>::EpetraCrsMatrixT(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &colMap, const ArrayRCP< const size_t > &NumEntriesPerRowToAlloc, ProfileType pftype, const Teuchos::RCP< Teuchos::ParameterList > &plist)
   : isFillResumed_(false)
 {
   Teuchos::Array<int> numEntriesPerRowToAlloc(NumEntriesPerRowToAlloc.begin(), NumEntriesPerRowToAlloc.end()); // convert array of "size_t" to array of "int"
   mtx_ = Teuchos::rcp(new Epetra_CrsMatrix(Copy, toEpetra(rowMap), toEpetra(colMap), numEntriesPerRowToAlloc.getRawPtr(), toEpetra(pftype)));
 }

//.........这里部分代码省略.........
        for (int i = 0; i < num_neighbors; i++)
          send_maps[i] = &(sendDatas[i][0]);

        // Debugging
        printMaps(comm, sendDatas, amgx2muelu, neighbors, *importer->getTargetMap(), "send_map_vector");

        // Construct recv arrays
        std::vector<std::vector<int> > recvDatas (num_neighbors);
        std::vector<int>               recv_sizes(num_neighbors, 0);
        for (int i = 0; i < importPIDs.size(); i++)
          if (importPIDs[i] != -1) {
            int index = hashTable.get(importPIDs[i]);
            recvDatas [index].push_back(muelu2amgx_[i]);
            recv_sizes[index]++;
        }
        // FIXME: recvDatas must be sorted (based on GIDs)

        std::vector<const int*> recv_maps(num_neighbors);
        for (int i = 0; i < num_neighbors; i++)
          recv_maps[i] = &(recvDatas[i][0]);

        // Debugging
        printMaps(comm, recvDatas, amgx2muelu, neighbors, *importer->getTargetMap(), "recv_map_vector");

        AMGX_SAFE_CALL(AMGX_matrix_comm_from_maps_one_ring(A_, 1, num_neighbors, neighbors, &send_sizes[0], &send_maps[0], &recv_sizes[0], &recv_maps[0]));

        AMGX_vector_bind(X_, A_);
        AMGX_vector_bind(Y_, A_);
      }

      RCP<Teuchos::Time> matrixTransformTimer = Teuchos::TimeMonitor::getNewTimer("MueLu: AMGX: transform matrix");
      matrixTransformTimer->start();

      ArrayRCP<const size_t> ia_s;
      ArrayRCP<const int>    ja;
      ArrayRCP<const double> a;
      inA->getAllValues(ia_s, ja, a);

      ArrayRCP<int> ia(ia_s.size());
      for (int i = 0; i < ia.size(); i++)
        ia[i] = Teuchos::as<int>(ia_s[i]);

      N_      = inA->getNodeNumRows();
      int nnz = inA->getNodeNumEntries();

      matrixTransformTimer->stop();
      matrixTransformTimer->incrementNumCalls();


      // Upload matrix
      // TODO Do we need to pin memory here through AMGX_pin_memory?
      RCP<Teuchos::Time> matrixTimer = Teuchos::TimeMonitor::getNewTimer("MueLu: AMGX: transfer matrix  CPU->GPU");
      matrixTimer->start();
      if (numProcs == 1) {
        AMGX_matrix_upload_all(A_, N_, nnz, 1, 1, &ia[0], &ja[0], &a[0], NULL);

      } else {
        // Transform the matrix
        std::vector<int>    ia_new(ia.size());
        std::vector<int>    ja_new(ja.size());
        std::vector<double> a_new (a.size());

        ia_new[0] = 0;
        for (int i = 0; i < N_; i++) {
          int oldRow = amgx2muelu[i];

  void globalWeightedCutsByPart(
    const RCP<const Environment> &env,
    const RCP<const Comm<int> > &comm,
    const RCP<const GraphModel<typename Adapter::base_adapter_t> > &graph,
    const ArrayView<const typename Adapter::part_t> &part,
    typename Adapter::part_t &numParts,
    ArrayRCP<RCP<BaseClassMetrics<typename Adapter::scalar_t> > > &metrics,
    ArrayRCP<typename Adapter::scalar_t> &globalSums)
{
  env->debug(DETAILED_STATUS, "Entering globalWeightedCutsByPart");
  //////////////////////////////////////////////////////////
  // Initialize return values

  numParts = 0;

  int ewgtDim = graph->getNumWeightsPerEdge();

  int numMetrics = 1;                   // "edge cuts"
  if (ewgtDim) numMetrics += ewgtDim;   // "weight n"

  typedef typename Adapter::scalar_t scalar_t;
  typedef typename Adapter::gno_t gno_t;
  typedef typename Adapter::lno_t lno_t;
  typedef typename Adapter::node_t node_t;
  typedef typename Adapter::part_t part_t;
  typedef StridedData<lno_t, scalar_t> input_t;

  typedef GraphMetrics<scalar_t> mv_t;
  typedef Tpetra::CrsMatrix<part_t,lno_t,gno_t,node_t>  sparse_matrix_type;
  typedef Tpetra::Vector<part_t,lno_t,gno_t,node_t>     vector_t;
  typedef Tpetra::Map<lno_t, gno_t, node_t>                map_type;
  typedef Tpetra::global_size_t GST;
  const GST INVALID = Teuchos::OrdinalTraits<GST>::invalid ();

  using Teuchos::as;

  // add some more metrics to the array
  typedef typename ArrayRCP<RCP<BaseClassMetrics<typename Adapter::scalar_t> > >::size_type array_size_type;
  metrics.resize( metrics.size() + numMetrics );

  for( array_size_type n = metrics.size() - numMetrics; n < metrics.size(); ++n ) {
    mv_t * newMetric = new mv_t;									// allocate the new memory
    env->localMemoryAssertion(__FILE__,__LINE__,1,newMetric);		// check errors
    metrics[n] = rcp( newMetric); 				// create the new members
    }
  array_size_type next = metrics.size() - numMetrics; // MDM - this is most likely temporary to preserve the format here - we are now filling a larger array so we may not have started at 0


  //////////////////////////////////////////////////////////
  // Figure out the global number of parts in use.
  // Verify number of vertex weights is the same everywhere.

  lno_t localNumObj = part.size();
  part_t localNum[2], globalNum[2];
  localNum[0] = static_cast<part_t>(ewgtDim);
  localNum[1] = 0;

  for (lno_t i=0; i < localNumObj; i++)
    if (part[i] > localNum[1]) localNum[1] = part[i];

  try{
    reduceAll<int, part_t>(*comm, Teuchos::REDUCE_MAX, 2,
      localNum, globalNum);
  }
  Z2_THROW_OUTSIDE_ERROR(*env)

  env->globalBugAssertion(__FILE__,__LINE__,
    "inconsistent number of edge weights",
    globalNum[0] == localNum[0], DEBUG_MODE_ASSERTION, comm);

  part_t nparts = globalNum[1] + 1;

  part_t globalSumSize = nparts * numMetrics;
  scalar_t * sumBuf = new scalar_t [globalSumSize];
  env->localMemoryAssertion(__FILE__, __LINE__, globalSumSize, sumBuf);
  globalSums = arcp(sumBuf, 0, globalSumSize);

  //////////////////////////////////////////////////////////
  // Calculate the local totals by part.

  scalar_t *localBuf = new scalar_t [globalSumSize];
  env->localMemoryAssertion(__FILE__,__LINE__,globalSumSize,localBuf);
  memset(localBuf, 0, sizeof(scalar_t) * globalSumSize);

  scalar_t *cut = localBuf;              // # of cuts

  ArrayView<const gno_t> Ids;
  ArrayView<input_t> vwgts;
  //size_t nv =
  graph->getVertexList(Ids, vwgts);

  ArrayView<const gno_t> edgeIds;
  ArrayView<const lno_t> offsets;
  ArrayView<input_t> wgts;
  //size_t numLocalEdges =
  graph->getEdgeList(edgeIds, offsets, wgts);
  // **************************************************************************
  // *************************** BUILD MAP FOR ADJS ***************************
  // **************************************************************************

//.........这里部分代码省略.........