C++ RCP::Comm方法代码示例

 void ISVDUDV::Initialize( const Teuchos::RCP< Teuchos::ParameterList >& params,
                           const Teuchos::RCP< const Epetra_MultiVector >& ss,
                           const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > >& fileio) 
 {
   workU_ = Teuchos::rcp( new Epetra_MultiVector(ss->Map(),maxBasisSize_,false) );
   Epetra_LocalMap lclmap(ss->NumVectors(),0,ss->Comm());
   workV_ = Teuchos::rcp( new Epetra_MultiVector(lclmap,maxBasisSize_,false) );
 }

ZoltanLibClass::ZoltanLibClass(Teuchos::RCP<const Epetra_MultiVector> input_coords,
                               Teuchos::RCP<const Epetra_MultiVector> weights,
                               int inputType):
  Library(input_coords, weights, inputType)
{
  int weightDim = weights->NumVectors();

  if (weightDim > 1){
    if (input_coords->Comm().MyPID() == 0){
      std::cout << "WARNING: Zoltan will only use the first weight of the "<< weightDim << " supplied for each object" << std::endl;
    }
  }
}

//EpetraMap_To_TpetraMap: takes in Epetra_Map object, converts it to its equivalent Tpetra::Map object,
//and returns an RCP pointer to this Tpetra::Map
Teuchos::RCP<const Tpetra_Map> Petra::EpetraMap_To_TpetraMap(const Teuchos::RCP<const Epetra_Map>& epetraMap_,
                                                      const Teuchos::RCP<const Teuchos::Comm<int> >& commT_)
{
  const std::size_t numElements = Teuchos::as<std::size_t>(epetraMap_->NumMyElements());
  const auto indexBase = Teuchos::as<GO>(epetraMap_->IndexBase());
  if (epetraMap_->DistributedGlobal() || epetraMap_->Comm().NumProc() == Teuchos::OrdinalTraits<int>::one()) {
    Teuchos::Array<Tpetra_GO> indices(numElements);
    int *epetra_indices = epetraMap_->MyGlobalElements();
    for(LO i=0; i < numElements; i++)
       indices[i] = epetra_indices[i];
    const Tpetra::global_size_t computeGlobalElements = Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
    return Teuchos::rcp(new Tpetra_Map(computeGlobalElements, indices, indexBase, commT_));
  } else {
    return Teuchos::rcp(new Tpetra_Map(numElements, indexBase, commT_, Tpetra::LocallyReplicated));
  }
}

LOCA::Epetra::CompactWYOp::CompactWYOp(
        const Teuchos::RCP<LOCA::GlobalData>& global_data,
	const Teuchos::RCP<const Epetra_Operator>& jacOperator, 
	const Teuchos::RCP<const Epetra_MultiVector>& A_multiVec, 
	const Teuchos::RCP<const Epetra_MultiVector>& Y_x_multiVec,
	const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& Y_p_matrix,
	const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& T_matrix) :
  globalData(global_data),
  label("LOCA::Epetra::CompactWYOp"),
  localMap(Y_x_multiVec->NumVectors(), 0, jacOperator->Comm()),
  J(jacOperator),
  A(A_multiVec),
  Y_x(Y_x_multiVec),
  Y_p(View, localMap, Y_p_matrix->values(), Y_p_matrix->stride(), 
      Y_p_matrix->numCols()),
  T(View, localMap, T_matrix->values(), T_matrix->stride(), 
    T_matrix->numCols()),
  tmpMat1(NULL),
  tmpMV(NULL),
  dblas()
{
}

NOX::Abstract::Group::ReturnType
LOCA::Epetra::Group::computeComplex(double frequency)
{
  std::string callingFunction = 
    "LOCA::Epetra::Group::computeComplex()";

  // We must have the time-dependent interface
  if (userInterfaceTime == Teuchos::null)
    return NOX::Abstract::Group::BadDependency;

#ifdef HAVE_NOX_EPETRAEXT
  if (isValidComplex)
    return NOX::Abstract::Group::Ok;

  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
  NOX::Abstract::Group::ReturnType status;

  // Get Jacobian matrix
  Teuchos::RCP<Epetra_RowMatrix> jac = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(sharedLinearSystem.getObject(this)->getJacobianOperator());
   
  // Create complex matrix
  if (complexMatrix == Teuchos::null) {
    std::vector< std::vector<int> >rowStencil(2);
    std::vector<int> rowIndex;

    rowStencil[0].push_back(0); rowStencil[0].push_back(1);
    rowStencil[1].push_back(-1);  rowStencil[1].push_back(0);
    rowIndex.push_back(0); rowIndex.push_back(1);

    complexMatrix = Teuchos::rcp(new EpetraExt::BlockCrsMatrix(*jac,
							       rowStencil, 
							       rowIndex, 
							       jac->Comm()));

    // Construct global solution vector, the overlap vector, and importer between them
   complexVec = 
     Teuchos::rcp(new EpetraExt::BlockVector(jac->RowMatrixRowMap(), 
					     complexMatrix->RowMap()));  
  }

  // Get mass matrix M
  Teuchos::RCP<Epetra_RowMatrix> mass = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(shiftedSharedLinearSystem->getObject(this)->getJacobianOperator());

  // Compute w*M
  status = computeShiftedMatrix(0.0, frequency);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							   finalStatus,
							   callingFunction);

  // Load w*M in complex matrix
  complexMatrix->LoadBlock(*mass, 1, 0);

  // Compute -w*M
  status = computeShiftedMatrix(0.0, -frequency);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							   finalStatus,
							   callingFunction);

  // Load -w*M in complex matrix
  complexMatrix->LoadBlock(*mass, 0, 1);

  // Compute J
  status = computeJacobian();
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							   finalStatus,
							   callingFunction);

  // Load J in complex matrix
  complexMatrix->LoadBlock(*jac, 0, 0);
  complexMatrix->LoadBlock(*jac, 1, 1);

  // Create linear system
  if (complexSharedLinearSystem == Teuchos::null) {

    NOX::Epetra::Vector nev(complexVec, NOX::Epetra::Vector::CreateView);

    Teuchos::RCP<Teuchos::ParameterList> lsParams = 
      globalData->parsedParams->getSublist("Linear Solver");

    // Create the Linear System
    Teuchos::RCP<NOX::Epetra::Interface::Required> iReq;
    Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac;
    Teuchos::RCP<NOX::Epetra::LinearSystem> complexLinSys = 
      Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(printParams, 
							*lsParams,
							iReq, 
							iJac, 
							complexMatrix, 
							nev));
    complexLinSys->setJacobianOperatorForSolve(complexMatrix);
    complexSharedLinearSystem = 
      Teuchos::rcp(new NOX::SharedObject<NOX::Epetra::LinearSystem, 
		                         LOCA::Epetra::Group>(complexLinSys));
  }

  if (finalStatus == NOX::Abstract::Group::Ok)
    isValidComplex = true;
//.........这里部分代码省略.........

  void IncSVDPOD::Initialize( const Teuchos::RCP< Teuchos::ParameterList >& params,
                              const Teuchos::RCP< const Epetra_MultiVector >& ss,
                              const Teuchos::RCP< RBGen::FileIOHandler< Epetra_Operator > >& fileio ) {

    using Teuchos::rcp;

    // Get the "Reduced Basis Method" sublist.
    Teuchos::ParameterList rbmethod_params = params->sublist( "Reduced Basis Method" );

    // Get the maximum basis size
    maxBasisSize_ = rbmethod_params.get<int>("Max Basis Size");
    TEUCHOS_TEST_FOR_EXCEPTION(maxBasisSize_ < 2,std::invalid_argument,"\"Max Basis Size\" must be at least 2.");

    // Get a filter
    filter_ = rbmethod_params.get<Teuchos::RCP<Filter<double> > >("Filter",Teuchos::null);
    if (filter_ == Teuchos::null) {
      int k = rbmethod_params.get("Rank",(int)5);
      filter_ = rcp( new RangeFilter<double>(LARGEST,k,k) );
    }

    // Get convergence tolerance
    tol_ = rbmethod_params.get<double>("Convergence Tolerance",tol_);

    // Get debugging flag
    debug_ = rbmethod_params.get<bool>("IncSVD Debug",debug_);

    // Get verbosity level
    verbLevel_ = rbmethod_params.get<int>("IncSVD Verbosity Level",verbLevel_);

    // Get an Anasazi orthomanager
    if (rbmethod_params.isType<
          Teuchos::RCP< Anasazi::OrthoManager<double,Epetra_MultiVector> > 
        >("Ortho Manager")
       ) 
    {
      ortho_ = rbmethod_params.get< 
                Teuchos::RCP<Anasazi::OrthoManager<double,Epetra_MultiVector> >
               >("Ortho Manager");
      TEUCHOS_TEST_FOR_EXCEPTION(ortho_ == Teuchos::null,std::invalid_argument,"User specified null ortho manager.");
    }
    else {
      std::string omstr = rbmethod_params.get("Ortho Manager","DGKS");
      if (omstr == "SVQB") {
        ortho_ = rcp( new Anasazi::SVQBOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
      }
      else { // if omstr == "DGKS"
        ortho_ = rcp( new Anasazi::BasicOrthoManager<double,Epetra_MultiVector,Epetra_Operator>() );
      }
    }

    // Lmin,Lmax,Kstart
    lmin_ = rbmethod_params.get("Min Update Size",1);
    TEUCHOS_TEST_FOR_EXCEPTION(lmin_ < 1 || lmin_ >= maxBasisSize_,std::invalid_argument,
                       "Method requires 1 <= min update size < max basis size.");
    lmax_ = rbmethod_params.get("Max Update Size",maxBasisSize_);
    TEUCHOS_TEST_FOR_EXCEPTION(lmin_ > lmax_,std::invalid_argument,"Max update size must be >= min update size.");

    startRank_ = rbmethod_params.get("Start Rank",lmin_);
    TEUCHOS_TEST_FOR_EXCEPTION(startRank_ < 1 || startRank_ > maxBasisSize_,std::invalid_argument,
                       "Starting rank must be in [1,maxBasisSize_)");
    // MaxNumPasses
    maxNumPasses_ = rbmethod_params.get("Maximum Number Passes",maxNumPasses_);
    TEUCHOS_TEST_FOR_EXCEPTION(maxNumPasses_ != -1 && maxNumPasses_ <= 0, std::invalid_argument,
                       "Maximum number passes must be -1 or > 0.");
    // Save the pointer to the snapshot matrix
    TEUCHOS_TEST_FOR_EXCEPTION(ss == Teuchos::null,std::invalid_argument,"Input snapshot matrix cannot be null.");
    A_ = ss;

    // MaxNumCols
    maxNumCols_ = A_->NumVectors();
    maxNumCols_ = rbmethod_params.get("Maximum Number Columns",maxNumCols_);
    TEUCHOS_TEST_FOR_EXCEPTION(maxNumCols_ < A_->NumVectors(), std::invalid_argument,
                       "Maximum number of columns must be at least as many as in the initializing data set.");

    // V locally replicated or globally distributed
    // this must be true for now
    // Vlocal_ = rbmethod_params.get("V Local",Vlocal_);

    // Allocate space for the factorization
    sigma_.reserve( maxBasisSize_ );
    U_ = Teuchos::null;
    V_ = Teuchos::null;
    U_ = rcp( new Epetra_MultiVector(ss->Map(),maxBasisSize_,false) );
    if (Vlocal_) {
      Epetra_LocalMap lclmap(maxNumCols_,0,ss->Comm());
      V_ = rcp( new Epetra_MultiVector(lclmap,maxBasisSize_,false) );
    }
    else {
      Epetra_Map gblmap(maxNumCols_,0,ss->Comm());
      V_ = rcp( new Epetra_MultiVector(gblmap,maxBasisSize_,false) );
    }
    B_ = rcp( new Epetra_SerialDenseMatrix(maxBasisSize_,maxBasisSize_) );
    resNorms_.reserve(maxBasisSize_);

    // clear counters
    numProc_ = 0;
    curNumPasses_ = 0;

    // we are now initialized, albeit with null rank
    isInitialized_ = true;
//.........这里部分代码省略.........

static int run_test(Teuchos::RCP<Epetra_CrsMatrix> matrix,
	  bool verbose,           // display the graph before & after
	  bool contract,          // set global number of partitions to 1/2 num procs
	  int partitioningType,   // hypergraph or graph partitioning, or simple
	  int vertexWeightType,   // use vertex weights?
	  int edgeWeightType,     // use edge/hyperedge weights?
	  int objectType)         // use isorropia's CrsMatrix or CrsGraph
{
  int rc=0, fail = 0;
#ifdef HAVE_EPETRAEXT
  int localProc = 0;
  double balance1, balance2, cutn1, cutn2, cutl1, cutl2;
  double balance3, cutn3, cutl3;
  double cutWgt1, cutWgt2, cutWgt3;
  int numCuts1, numCuts2, numCuts3, valid;
  int numPartitions = 0;
  int keepDenseEdges = 0;
  int numProcs = 1;

#ifdef HAVE_MPI
  const Epetra_MpiComm &Comm = dynamic_cast<const Epetra_MpiComm &>(matrix->Comm());
  localProc = Comm.MyPID();
  numProcs = Comm.NumProc();
#else
  const Epetra_SerialComm &Comm = dynamic_cast<const Epetra_SerialComm &>(matrix->Comm());
#endif

  int numRows = matrix->NumGlobalRows();

  if (numRows < (numProcs * 100)){
    // By default Zoltan throws out dense edges, defined as those
    // whose number of non-zeros exceeds 25% of the number of vertices.
    //
    // If dense edges are thrown out of a small matrix, there may be nothing left.
    keepDenseEdges = 1;
  }

  double myShareBefore = 1.0 / numProcs;
  double myShare = myShareBefore;

  if (contract){
    numPartitions = numProcs / 2;

    if (numPartitions > numRows)
      numPartitions = numRows;

    if (numPartitions > 0){
      if (localProc < numPartitions){
	myShare = 1.0 / numPartitions;
      }
      else{
	myShare = 0.0;
      }
    }
    else{
      contract = 0;
    }
  }

  // If we want Zoltan's or Isorropia's default weights, then we don't
  // need to supply a CostDescriber object to createBalancedCopy,
  // so we get to test the API functions that don't take a CostDescriber.

  bool noCosts = ((vertexWeightType == NO_APPLICATION_SUPPLIED_WEIGHTS) &&
		   (edgeWeightType == NO_APPLICATION_SUPPLIED_WEIGHTS));

  // Test the interface that has no parameters, if possible

  bool noParams =
    ((partitioningType == HYPERGRAPH_PARTITIONING) && // default, so requires no params
     (numPartitions == 0) &&                          // >0 would require a parameter
     (keepDenseEdges == 0));                          // >0 would require a parameter

  // Maps for original object
  const Epetra_Map &sourceRowMap = matrix->RowMap();
  const Epetra_Map &sourceRangeMap = matrix->RangeMap();
//   const Epetra_Map &sourceColMap = matrix->ColMap();
  const Epetra_Map &sourceDomainMap = matrix->DomainMap();

  int numCols = matrix->NumGlobalCols();
  int nMyRows = sourceRowMap.NumMyElements();
  int base = sourceRowMap.IndexBase();

  // Compute vertex and edge weights

  Isorropia::Epetra::CostDescriber costs;

  Teuchos::RCP<Epetra_Vector> vptr;

  Teuchos::RCP<Epetra_CrsMatrix> eptr;

  Teuchos::RCP<Epetra_Vector> hyperEdgeWeights;

  if (edgeWeightType != NO_APPLICATION_SUPPLIED_WEIGHTS){

    if (partitioningType == GRAPH_PARTITIONING){

      // Create graph edge weights.

      eptr = Teuchos::rcp(new Epetra_CrsMatrix(*matrix));
//.........这里部分代码省略.........

// helper routines
bool SplitMatrix2x2(Teuchos::RCP<const Epetra_CrsMatrix> A,
    const Epetra_Map& A11rowmap,
    const Epetra_Map& A22rowmap,
    Teuchos::RCP<Epetra_CrsMatrix>& A11,
    Teuchos::RCP<Epetra_CrsMatrix>& A12,
    Teuchos::RCP<Epetra_CrsMatrix>& A21,
    Teuchos::RCP<Epetra_CrsMatrix>& A22)
{
  if (A==Teuchos::null)
  {
    std::cout << "ERROR: SplitMatrix2x2: A==null on entry" << std::endl;
    return false;
  }

  const Epetra_Comm& Comm   = A->Comm();
  const Epetra_Map&  A22map = A22rowmap;
  const Epetra_Map&  A11map = A11rowmap;

  //----------------------------- create a parallel redundant map of A22map
  std::map<int,int> a22gmap;
  {
    std::vector<int> a22global(A22map.NumGlobalElements());
    int count=0;
    for (int proc=0; proc<Comm.NumProc(); ++proc)
    {
      int length = 0;
      if (proc==Comm.MyPID())
      {
        for (int i=0; i<A22map.NumMyElements(); ++i)
        {
          a22global[count+length] = A22map.GID(i);
          ++length;
        }
      }
      Comm.Broadcast(&length,1,proc);
      Comm.Broadcast(&a22global[count],length,proc);
      count += length;
    }
    if (count != A22map.NumGlobalElements())
    {
      std::cout << "ERROR SplitMatrix2x2: mismatch in dimensions" << std::endl;
      return false;
    }

    // create the map
    for (int i=0; i<count; ++i)
      a22gmap[a22global[i]] = 1;
    a22global.clear();
  }

  //--------------------------------------------------- create matrix A22
  A22 = Teuchos::rcp(new Epetra_CrsMatrix(Copy,A22map,100));
  {
    std::vector<int>    a22gcindices(100);
    std::vector<double> a22values(100);
    for (int i=0; i<A->NumMyRows(); ++i)
    {
      const int grid = A->GRID(i);
      if (A22map.MyGID(grid)==false)
        continue;
      int     numentries;
      double* values;
      int*    cindices;
      int err = A->ExtractMyRowView(i,numentries,values,cindices);
      if (err)
      {
        std::cout << "ERROR: SplitMatrix2x2: A->ExtractMyRowView returned " << err << std::endl;
        return false;
      }
      if (numentries>(int)a22gcindices.size())
      {
        a22gcindices.resize(numentries);
        a22values.resize(numentries);
      }
      int count=0;
      for (int j=0; j<numentries; ++j)
      {
        const int gcid = A->ColMap().GID(cindices[j]);
        // see whether we have gcid in a22gmap
        std::map<int,int>::iterator curr = a22gmap.find(gcid);
        if (curr==a22gmap.end()) continue;
        //std::cout << gcid << " ";
        a22gcindices[count] = gcid;
        a22values[count]    = values[j];
        ++count;
      }
      //std::cout << std::endl; fflush(stdout);
      // add this filtered row to A22
      err = A22->InsertGlobalValues(grid,count,&a22values[0],&a22gcindices[0]);
      if (err<0)
      {
        std::cout << "ERROR: SplitMatrix2x2: A->InsertGlobalValues returned " << err << std::endl;
        return false;
      }

    } //for (int i=0; i<A->NumMyRows(); ++i)
    a22gcindices.clear();
    a22values.clear();
  }
//.........这里部分代码省略.........

void PeridigmNS::Block::createMapsFromGlobalMaps(Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapScalarPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOverlapVectorPointMap,
                                                 Teuchos::RCP<const Epetra_BlockMap> globalOwnedScalarBondMap,
                                                 Teuchos::RCP<const Epetra_Vector>   globalBlockIds,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalNeighborhoodData,
                                                 Teuchos::RCP<const PeridigmNS::NeighborhoodData> globalContactNeighborhoodData)
{
  double* globalBlockIdsPtr;
  globalBlockIds->ExtractView(&globalBlockIdsPtr);

  // Create a list of all the on-processor elements that are part of this block

  vector<int> IDs;
  IDs.reserve(globalOverlapScalarPointMap->NumMyElements()); // upper bound
  vector<int> bondIDs;
  bondIDs.reserve(globalOverlapScalarPointMap->NumMyElements());
  vector<int> bondElementSize;
  bondElementSize.reserve(globalOwnedScalarPointMap->NumMyElements());

  for(int iLID=0 ; iLID<globalOwnedScalarPointMap->NumMyElements() ; ++iLID){
    if(globalBlockIdsPtr[iLID] == blockID) {
      int globalID = globalOwnedScalarPointMap->GID(iLID);
      IDs.push_back(globalID);
    }
  }

  // Record the size of these elements in the bond map
  // Note that if an element has no bonds, it has no entry in the bondMap
  // So, the bond map and the scalar map can have a different number of entries (different local IDs)

  for(int iLID=0 ; iLID<globalOwnedScalarBondMap->NumMyElements() ; ++iLID){
    int globalID = globalOwnedScalarBondMap->GID(iLID);
    int localID = globalOwnedScalarPointMap->LID(globalID);
    if(globalBlockIdsPtr[localID] == blockID){
      bondIDs.push_back(globalID);
      bondElementSize.push_back(globalOwnedScalarBondMap->ElementSize(iLID));
    }
  }

  // Create the owned scalar point map, the owned vector point map, and the owned scalar bond map

  int numGlobalElements = -1;
  int numMyElements = IDs.size();
  int* myGlobalElements = 0;
  if(numMyElements > 0)
    myGlobalElements = &IDs.at(0);
  int elementSize = 1;
  int indexBase = 0;
  ownedScalarPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  elementSize = 3;
  ownedVectorPointMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSize, indexBase, globalOwnedScalarPointMap->Comm()));

  numMyElements = bondElementSize.size();
  myGlobalElements = 0;
  int* elementSizeList = 0;
  if(numMyElements > 0){
    myGlobalElements = &bondIDs.at(0);
    elementSizeList = &bondElementSize.at(0);
  }
  ownedScalarBondMap =
    Teuchos::rcp(new Epetra_BlockMap(numGlobalElements, numMyElements, myGlobalElements, elementSizeList, indexBase, globalOwnedScalarPointMap->Comm()));

  // Create a list of nodes that need to be ghosted (both across material boundaries and across processor boundaries)
  set<int> ghosts;

  // Check the neighborhood list for things that need to be ghosted
  int* const globalNeighborhoodList = globalNeighborhoodData->NeighborhoodList();
  int globalNeighborhoodListIndex = 0;
  for(int iLID=0 ; iLID<globalNeighborhoodData->NumOwnedPoints() ; ++iLID){
    int numNeighbors = globalNeighborhoodList[globalNeighborhoodListIndex++];
    if(globalBlockIdsPtr[iLID] == blockID) {
      for(int i=0 ; i<numNeighbors ; ++i){
        int neighborGlobalID = globalOverlapScalarPointMap->GID( globalNeighborhoodList[globalNeighborhoodListIndex + i] );
        ghosts.insert(neighborGlobalID);
      }
    }
    globalNeighborhoodListIndex += numNeighbors;
  }

  // Check the contact neighborhood list for things that need to be ghosted
  if(!globalContactNeighborhoodData.is_null()){
    int* const globalContactNeighborhoodList = globalContactNeighborhoodData->NeighborhoodList();
    int globalContactNeighborhoodListIndex = 0;
    for(int iLID=0 ; iLID<globalContactNeighborhoodData->NumOwnedPoints() ; ++iLID){
      int numNeighbors = globalContactNeighborhoodList[globalContactNeighborhoodListIndex++];
      if(globalBlockIdsPtr[iLID] == blockID) {
        for(int i=0 ; i<numNeighbors ; ++i){
          int neighborGlobalID = globalOverlapScalarPointMap->GID( globalContactNeighborhoodList[globalContactNeighborhoodListIndex + i] );
          ghosts.insert(neighborGlobalID);
        }
      }
      globalContactNeighborhoodListIndex += numNeighbors;
    }
  }

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

AmesosBTFGlobal_LinearProblem::NewTypeRef
AmesosBTFGlobal_LinearProblem::
operator()( OriginalTypeRef orig )
{
  origObj_ = &orig;

  // Extract the matrix and vectors from the linear problem
  OldRHS_ = Teuchos::rcp( orig.GetRHS(), false );
  OldLHS_ = Teuchos::rcp( orig.GetLHS(), false );
  OldMatrix_ = Teuchos::rcp( dynamic_cast<Epetra_CrsMatrix *>( orig.GetMatrix() ), false );
	
  int nGlobal = OldMatrix_->NumGlobalRows(); 
  int n = OldMatrix_->NumMyRows();

  // Check if the matrix is on one processor.
  int myMatProc = -1, matProc = -1;
  int myPID = OldMatrix_->Comm().MyPID();
  int numProcs = OldMatrix_->Comm().NumProc();

  const Epetra_BlockMap& oldRowMap = OldMatrix_->RowMap();

  // Get some information about the parallel distribution.
  int maxMyRows = 0;
  std::vector<int> numGlobalElem( numProcs );
  OldMatrix_->Comm().GatherAll(&n, &numGlobalElem[0], 1);
  OldMatrix_->Comm().MaxAll(&n, &maxMyRows, 1);

  for (int proc=0; proc<numProcs; proc++) 
  {
    if (OldMatrix_->NumGlobalNonzeros() == OldMatrix_->NumMyNonzeros())
      myMatProc = myPID;
  }

  OldMatrix_->Comm().MaxAll( &myMatProc, &matProc, 1 );

  Teuchos::RCP<Epetra_CrsMatrix> serialMatrix;
  Teuchos::RCP<Epetra_Map> serialMap;	
  if( oldRowMap.DistributedGlobal() && matProc == -1) 
  {
    // The matrix is distributed and needs to be moved to processor zero.
    // Set the zero processor as the master.
    matProc = 0;
    serialMap = Teuchos::rcp( new Epetra_Map( Epetra_Util::Create_Root_Map( OldMatrix_->RowMap(), matProc ) ) );
    
    Epetra_Import serialImporter( *serialMap, OldMatrix_->RowMap() );
    serialMatrix = Teuchos::rcp( new Epetra_CrsMatrix( Copy, *serialMap, 0 ) );
    serialMatrix->Import( *OldMatrix_, serialImporter, Insert );
    serialMatrix->FillComplete();
  }
  else {
    // The old matrix has already been moved to one processor (matProc).
    serialMatrix = OldMatrix_;
  }

  if( debug_ )
  {
    cout << "Original (serial) Matrix:\n";
    cout << *serialMatrix << endl;
  }

  // Obtain the current row and column orderings
  std::vector<int> origGlobalRows(nGlobal), origGlobalCols(nGlobal);
  serialMatrix->RowMap().MyGlobalElements( &origGlobalRows[0] );
  serialMatrix->ColMap().MyGlobalElements( &origGlobalCols[0] );
  
  // Perform reindexing on the full serial matrix (needed for BTF).
  Epetra_Map reIdxMap( serialMatrix->RowMap().NumGlobalElements(), serialMatrix->RowMap().NumMyElements(), 0, serialMatrix->Comm() );
  Teuchos::RCP<EpetraExt::ViewTransform<Epetra_CrsMatrix> > reIdxTrans =
    Teuchos::rcp( new EpetraExt::CrsMatrix_Reindex( reIdxMap ) );
  Epetra_CrsMatrix newSerialMatrix = (*reIdxTrans)( *serialMatrix );
  reIdxTrans->fwd();
  
  // Compute and apply BTF to the serial CrsMatrix and has been filtered by the threshold
  EpetraExt::AmesosBTF_CrsMatrix BTFTrans( threshold_, upperTri_, verbose_, debug_ );
  Epetra_CrsMatrix newSerialMatrixBTF = BTFTrans( newSerialMatrix );
  
  rowPerm_ = BTFTrans.RowPerm();
  colPerm_ = BTFTrans.ColPerm();
  blockPtr_ = BTFTrans.BlockPtr();
  numBlocks_ = BTFTrans.NumBlocks();
 
  if (myPID == matProc && verbose_) {
    bool isSym = true;
    for (int i=0; i<nGlobal; ++i) {
      if (rowPerm_[i] != colPerm_[i]) {
        isSym = false;
        break;
      }
    }
    std::cout << "The BTF permutation symmetry (0=false,1=true) is : " << isSym << std::endl;
  }
  
  // Compute the permutation w.r.t. the original row and column GIDs.
  std::vector<int> origGlobalRowsPerm(nGlobal), origGlobalColsPerm(nGlobal);
  if (myPID == matProc) {
    for (int i=0; i<nGlobal; ++i) {
      origGlobalRowsPerm[i] = origGlobalRows[ rowPerm_[i] ];
      origGlobalColsPerm[i] = origGlobalCols[ colPerm_[i] ];
    }
  }
//.........这里部分代码省略.........