C++ Epetra_CrsGraph::FillComplete方法代码示例

//============================================================================
Epetra_CrsGraph* Ifpack_CreateOverlappingCrsMatrix(const Epetra_CrsGraph* Graph,
                                                   const int OverlappingLevel)
{

  if (OverlappingLevel == 0)
    return(0); // All done
  if (Graph->Comm().NumProc() == 1)
    return(0); // All done

  Epetra_CrsGraph* OverlappingGraph;
  Epetra_BlockMap* OverlappingMap;
  OverlappingGraph = const_cast<Epetra_CrsGraph*>(Graph);
  OverlappingMap = const_cast<Epetra_BlockMap*>(&(Graph->RowMap()));

  Epetra_CrsGraph* OldGraph;
  Epetra_BlockMap* OldMap;
  const Epetra_BlockMap* DomainMap = &(Graph->DomainMap());
  const Epetra_BlockMap* RangeMap = &(Graph->RangeMap());

  for (int level = 1; level <= OverlappingLevel ; ++level) {

    OldGraph = OverlappingGraph;
    OldMap = OverlappingMap;

    Epetra_Import* OverlappingImporter;
    OverlappingImporter = const_cast<Epetra_Import*>(OldGraph->Importer());
    OverlappingMap = new Epetra_BlockMap(OverlappingImporter->TargetMap());

    if (level < OverlappingLevel)
      OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap, 0);
    else
      // On last iteration, we want to filter out all columns except
      // those that correspond
      // to rows in the graph.  This assures that our matrix is square
      OverlappingGraph = new Epetra_CrsGraph(Copy, *OverlappingMap,
                                          *OverlappingMap, 0);

    OverlappingGraph->Import(*OldGraph, *OverlappingImporter, Insert);
    if (level < OverlappingLevel)
      OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
    else {
      // Copy last OverlapImporter because we will use it later
      OverlappingImporter = new Epetra_Import(*OverlappingMap, *DomainMap);
      OverlappingGraph->FillComplete(*DomainMap, *RangeMap);
    }

    if (level > 1) {
      delete OldGraph;
      delete OldMap;
    }

    delete OverlappingMap;
    OverlappingGraph->FillComplete();

  }

  return(OverlappingGraph);
}

/*----------------------------------------------------------------------*
 |  make a deep copy of a graph                              m.gee 01/05|
 |  allocate the new graph                                              |
 *----------------------------------------------------------------------*/
Epetra_CrsGraph* ML_NOX::deepcopy_graph(const Epetra_CrsGraph* oldgraph)
{
   int  i,ierr;
   int  nrows = oldgraph->NumMyRows();
   int* nIndicesperRow = new int[nrows];

   for (i=0; i<nrows; i++)
      nIndicesperRow[i] = oldgraph->NumMyIndices(i);
   Epetra_CrsGraph* graph = new Epetra_CrsGraph(Copy,oldgraph->RowMap(),oldgraph->ColMap(),
                                                &(nIndicesperRow[0]));
   delete [] nIndicesperRow;
   nIndicesperRow = 0;
   
   for (i=0; i<nrows; i++)
   {
      int  numIndices;
      int* Indices=0;
      ierr = oldgraph->ExtractMyRowView(i,numIndices,Indices);
      ierr = graph->InsertMyIndices(i,numIndices,Indices);
   }

   graph->FillComplete();
   return graph;
}                                                     

int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
{   
  if (verbose) {
    cout << "******************* four_quads ***********************"<<endl;
  }

  //This function assembles a matrix representing a finite-element mesh
  //of four 2-D quad elements. There are 9 nodes in the problem. The
  //same problem is assembled no matter how many processors are being used
  //(within reason). It may not work if more than 9 processors are used.
  //
  //  *------*------*
  // 6|     7|     8|
  //  | E2   | E3   |
  //  *------*------*
  // 3|     4|     5|
  //  | E0   | E1   |
  //  *------*------*
  // 0      1      2
  //
  //Nodes are denoted by * with node-numbers below and left of each node.
  //E0, E1 and so on are element-numbers.
  //
  //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
  //for each element. Thus, the coefficient value at position 0,0 should end up
  //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
  //
  //Depending on the number of processors being used, the locations of the
  //specific matrix positions (in terms of which processor owns them) will vary.
  //
  
  int numProcs = Comm.NumProc();
  
  int numNodes = 9;
  int numElems = 4;
  int numNodesPerElem = 4;

  int blockSize = 1;
  int indexBase = 0;

  //Create a map using epetra-defined linear distribution.
  Epetra_BlockMap map(numNodes, blockSize, indexBase, Comm);

  Epetra_CrsGraph* graph = NULL;

  int* nodes = new int[numNodesPerElem];
  int i, j, k, err = 0;

  if (preconstruct_graph) {
    graph = new Epetra_CrsGraph(Copy, map, 1);

    //we're going to fill the graph with indices, but remember it will only
    //accept indices in rows for which map.MyGID(row) is true.

    for(i=0; i<numElems; ++i) {
      switch(i) {
      case 0:
        nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
        break;
      case 1:
        nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
        break;
      case 2:
        nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
        break;
      case 3:
        nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
        break;
      }

      for(j=0; j<numNodesPerElem; ++j) {
        if (map.MyGID(nodes[j])) {
          err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
                                           nodes);
          if (err<0) return(err);
        }
      }
    }

    EPETRA_CHK_ERR( graph->FillComplete() );
  }

  Epetra_FEVbrMatrix* A = NULL;

  if (preconstruct_graph) {
    A = new Epetra_FEVbrMatrix(Copy, *graph);
  }
  else {
    A = new Epetra_FEVbrMatrix(Copy, map, 1);
  }

  //EPETRA_CHK_ERR( A->PutScalar(0.0) );

  double* values_1d = new double[numNodesPerElem*numNodesPerElem];
  double** values_2d = new double*[numNodesPerElem];

  for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;

  int offset = 0;
  for(i=0; i<numNodesPerElem; ++i) {
//.........这里部分代码省略.........

Epetra_CrsGraph * BlockUtility::TGenerateBlockGraph(
        const Epetra_CrsGraph & BaseGraph,
        const Epetra_CrsGraph & LocalBlockGraph,
        const Epetra_Comm & GlobalComm )
{
  const Epetra_BlockMap & BaseRowMap = BaseGraph.RowMap();
  const Epetra_BlockMap & BaseColMap = BaseGraph.ColMap();
  int_type ROffset = BlockUtility::TCalculateOffset<int_type>(BaseRowMap);
  (void) ROffset; // Silence "unused variable" compiler warning.
  int_type COffset = BlockUtility::TCalculateOffset<int_type>(BaseColMap);

  //Get Base Global IDs
  const Epetra_BlockMap & BlockRowMap = LocalBlockGraph.RowMap();
  const Epetra_BlockMap & BlockColMap = LocalBlockGraph.ColMap();

  int NumBlockRows = BlockRowMap.NumMyElements();
  vector<int_type> RowIndices(NumBlockRows);
  BlockRowMap.MyGlobalElements(&RowIndices[0]);

  int Size = BaseRowMap.NumMyElements();

  Epetra_Map *GlobalRowMap =
    GenerateBlockMap(BaseRowMap, BlockRowMap, GlobalComm);


  int MaxIndices = BaseGraph.MaxNumIndices();
  vector<int_type> Indices(MaxIndices);

  Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy,
                               dynamic_cast<Epetra_BlockMap&>(*GlobalRowMap),
                               0 );

  int NumBlockIndices, NumBaseIndices;
  int *BlockIndices, *BaseIndices;
  for( int i = 0; i < NumBlockRows; ++i )
  {
    LocalBlockGraph.ExtractMyRowView(i, NumBlockIndices, BlockIndices);

    for( int j = 0; j < Size; ++j )
    {
      int_type GlobalRow = (int_type) GlobalRowMap->GID64(j+i*Size);

      BaseGraph.ExtractMyRowView( j, NumBaseIndices, BaseIndices );
      for( int k = 0; k < NumBlockIndices; ++k )
      {
        int_type ColOffset = (int_type) BlockColMap.GID64(BlockIndices[k]) * COffset;

        for( int l = 0; l < NumBaseIndices; ++l )
          Indices[l] = (int_type) BaseGraph.GCID64(BaseIndices[l]) + ColOffset;

        GlobalGraph->InsertGlobalIndices( GlobalRow, NumBaseIndices, &Indices[0] );
      }
    }
  }

  const Epetra_BlockMap & BaseDomainMap = BaseGraph.DomainMap();
  const Epetra_BlockMap & BaseRangeMap = BaseGraph.RangeMap();
  const Epetra_BlockMap & BlockDomainMap = LocalBlockGraph.DomainMap();
  const Epetra_BlockMap & BlockRangeMap = LocalBlockGraph.RangeMap();

  Epetra_Map *GlobalDomainMap =
    GenerateBlockMap(BaseDomainMap, BlockDomainMap, GlobalComm);
  Epetra_Map *GlobalRangeMap =
    GenerateBlockMap(BaseRangeMap, BlockRangeMap, GlobalComm);

  GlobalGraph->FillComplete(*GlobalDomainMap, *GlobalRangeMap);

  delete GlobalDomainMap;
  delete GlobalRangeMap;
  delete GlobalRowMap;

  return GlobalGraph;
}

Epetra_CrsGraph * BlockUtility::TGenerateBlockGraph(
        const Epetra_RowMatrix & BaseMatrix,
        const vector< vector<int_type> > & RowStencil,
        const vector<int_type> & RowIndices,
        const Epetra_Comm & GlobalComm )
{

  const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
  const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();
  int_type BaseIndex = (int_type) BaseMap.IndexBase64();
  int_type Offset = BlockUtility::TCalculateOffset<int_type>(BaseMap);

  //Get Base Global IDs
  int NumBlockRows = RowIndices.size();
  int Size = BaseMap.NumMyElements();
  int TotalSize = NumBlockRows * Size;
  vector<int_type> GIDs(Size);
  BaseMap.MyGlobalElements( &GIDs[0] );

  vector<int_type> GlobalGIDs( TotalSize );
  for( int i = 0; i < NumBlockRows; ++i )
  {
    for( int j = 0; j < Size; ++j )
      GlobalGIDs[i*Size+j] = GIDs[j] + RowIndices[i] * Offset;
  }

  int_type GlobalSize;
  int_type TotalSize_int_type = TotalSize;
  GlobalComm.SumAll( &TotalSize_int_type, &GlobalSize, 1 );

  Epetra_Map GlobalMap( GlobalSize, TotalSize, &GlobalGIDs[0], BaseIndex, GlobalComm );

  int MaxIndices = BaseMatrix.MaxNumEntries();
  vector<int> Indices_local(MaxIndices);
  vector<int_type> Indices_global(MaxIndices);
  vector<double> Values(MaxIndices);
  int NumIndices;

  Epetra_CrsGraph * GlobalGraph = new Epetra_CrsGraph( Copy,
                               dynamic_cast<Epetra_BlockMap&>(GlobalMap),
                               0 );

  for( int i = 0; i < NumBlockRows; ++i )
  {
    int StencilSize = RowStencil[i].size();
    for( int j = 0; j < Size; ++j )
    {
      int_type GlobalRow = (int_type) GlobalMap.GID64(j+i*Size);

      BaseMatrix.ExtractMyRowCopy( j, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );
      for( int l = 0; l < NumIndices; ++l ) Indices_global[l] = (int_type) BaseColMap.GID64(Indices_local[l]);

      for( int k = 0; k < StencilSize; ++k )
      {
        int_type ColOffset = (RowIndices[i]+RowStencil[i][k]) * Offset;
        if( k > 0 ) ColOffset -= (RowIndices[i]+RowStencil[i][k-1]) * Offset;

        for( int l = 0; l < NumIndices; ++l )
          Indices_global[l] += ColOffset;

        GlobalGraph->InsertGlobalIndices( GlobalRow, NumIndices, &Indices_global[0] );
      }
    }
  }

  GlobalGraph->FillComplete();

  return GlobalGraph;
}

int four_quads(const Epetra_Comm& Comm, bool preconstruct_graph, bool verbose)
{
  if (verbose) {
    cout << "******************* four_quads ***********************"<<endl;
  }

  //This function assembles a matrix representing a finite-element mesh
  //of four 2-D quad elements. There are 9 nodes in the problem. The
  //same problem is assembled no matter how many processors are being used
  //(within reason). It may not work if more than 9 processors are used.
  //
  //  *------*------*
  // 6|     7|     8|
  //  | E2   | E3   |
  //  *------*------*
  // 3|     4|     5|
  //  | E0   | E1   |
  //  *------*------*
  // 0      1      2
  //
  //Nodes are denoted by * with node-numbers below and left of each node.
  //E0, E1 and so on are element-numbers.
  //
  //Each processor will contribute a sub-matrix of size 4x4, filled with 1's,
  //for each element. Thus, the coefficient value at position 0,0 should end up
  //being 1.0*numProcs, the value at position 4,4 should be 1.0*4*numProcs, etc.
  //
  //Depending on the number of processors being used, the locations of the
  //specific matrix positions (in terms of which processor owns them) will vary.
  //

  int numProcs = Comm.NumProc();

  int numNodes = 9;
  int numElems = 4;
  int numNodesPerElem = 4;

  int indexBase = 0;

  //Create a map using epetra-defined linear distribution.
  Epetra_Map map(numNodes, indexBase, Comm);

  Epetra_CrsGraph* graph = NULL;

  int* nodes = new int[numNodesPerElem];
  int i, j, err = 0;

  if (preconstruct_graph) {
    graph = new Epetra_CrsGraph(Copy, map, 1);

    //we're going to fill the graph with indices, but remember it will only
    //accept indices in rows for which map.MyGID(row) is true.

    for(i=0; i<numElems; ++i) {
      switch(i) {
      case 0:
	nodes[0] = 0; nodes[1] = 1; nodes[2] = 4; nodes[3] = 3;
	break;
      case 1:
	nodes[0] = 1; nodes[1] = 2; nodes[2] = 5; nodes[3] = 4;
	break;
      case 2:
	nodes[0] = 3; nodes[1] = 4; nodes[2] = 7; nodes[3] = 6;
	break;
      case 3:
	nodes[0] = 4; nodes[1] = 5; nodes[2] = 8; nodes[3] = 7;
	break;
      }

      for(j=0; j<numNodesPerElem; ++j) {
	if (map.MyGID(nodes[j])) {
	  err = graph->InsertGlobalIndices(nodes[j], numNodesPerElem,
					   nodes);
	  if (err<0) return(err);
	}
      }
    }

    EPETRA_CHK_ERR( graph->FillComplete() );
  }

  Epetra_FECrsMatrix* A = NULL;

  if (preconstruct_graph) {
    A = new Epetra_FECrsMatrix(Copy, *graph);
  }
  else {
    A = new Epetra_FECrsMatrix(Copy, map, 1);
  }

  EPETRA_CHK_ERR( A->PutScalar(0.0) );

  double* values_1d = new double[numNodesPerElem*numNodesPerElem];
  double** values_2d = new double*[numNodesPerElem];

  for(i=0; i<numNodesPerElem*numNodesPerElem; ++i) values_1d[i] = 1.0;

  int offset = 0;
  for(i=0; i<numNodesPerElem; ++i) {
    values_2d[i] = &(values_1d[offset]);
//.........这里部分代码省略.........

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

  int nRows = orig.NumMyRows();
  int nCols = orig.NumMyCols();

  const Epetra_BlockMap & RowMap = orig.RowMap();

  int numIndices;
  int * Indices;

  Epetra_CrsGraph * TransposeGraph = 0;

  if( !ignoreNonLocalCols_ && orig.DistributedGlobal() )
  {
    std::vector<int> TransNumNZ( nCols, 0 );
    for( int i = 0; i < nRows; ++i )
    {
      orig.ExtractMyRowView( i, numIndices, Indices );
      for( int j = 0; j < numIndices; ++j ) ++TransNumNZ[ Indices[j] ];
    }

    std::vector< std::vector<int> > TransIndices( nCols );
    for( int i = 0; i < nCols; ++i )
      if( TransNumNZ[i] )
      {
        TransIndices[i].resize( TransNumNZ[i] );
        TransNumNZ[i] = 0;
      }

    for( int i = 0; i < nRows; ++i )
    {
      orig.ExtractMyRowView( i, numIndices, Indices );
      for( int j = 0; j < numIndices; ++j )
        TransIndices[ Indices[j] ][ TransNumNZ[ Indices[j] ]++ ] = i;
    }

    Epetra_CrsGraph SharedTransGraph( View, orig.ImportMap(), RowMap, &TransNumNZ[0] );
    for( int i = 0; i < nCols; ++i )
      if( TransNumNZ[i] ) SharedTransGraph.InsertMyIndices( i, TransNumNZ[i], &TransIndices[i][0] );
    SharedTransGraph.FillComplete();

    TransposeGraph = new Epetra_CrsGraph( Copy, RowMap, 0 );
    Epetra_Export Exporter( orig.ImportMap(), RowMap ); 
    TransposeGraph->Export( SharedTransGraph, Exporter, Add );
    TransposeGraph->FillComplete();
  }
  else
  {
    std::vector<int> TransNumNZ( nRows, 0 );
    for( int i = 0; i < nRows; ++i )
    {
      orig.ExtractMyRowView( i, numIndices, Indices );
      for( int j = 0; j < numIndices; ++j )
        if( Indices[j] < nRows ) ++TransNumNZ[ Indices[j] ];
    }

    std::vector< std::vector<int> > TransIndices( nRows );
    for( int i = 0; i < nRows; ++i )
      if( TransNumNZ[i] )
      {
        TransIndices[i].resize( TransNumNZ[i] );
        TransNumNZ[i] = 0;
      }

    for( int i = 0; i < nRows; ++i )
    {
      orig.ExtractMyRowView( i, numIndices, Indices );
      for( int j = 0; j < numIndices; ++j )
        if( Indices[j] < nRows ) TransIndices[ Indices[j] ][ TransNumNZ[ Indices[j] ]++ ] = i;
    }

    TransposeGraph = new Epetra_CrsGraph( Copy, RowMap, RowMap, &TransNumNZ[0] );

    for( int i = 0; i < nRows; ++i )
      if( TransNumNZ[i] ) TransposeGraph->InsertMyIndices( i, TransNumNZ[i], &TransIndices[i][0] );

    TransposeGraph->FillComplete();
  }

  newObj_ = TransposeGraph;

  return *TransposeGraph;
}

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

  int err;

  //Setup Load Balance Object
  float version;
  char * dummy = 0;
  Zoltan::LoadBalance LB( 0, &dummy, &version );
  err = LB.Create( dynamic_cast<const Epetra_MpiComm&>(orig.Comm()).Comm() );
  if( err == ZOLTAN_OK ) err = LB.Set_Param( "LB_METHOD", "GRAPH" );
#ifdef HAVE_LIBPARMETIS
  if( err == ZOLTAN_OK ) err = LB.Set_Param( "GRAPH_PACKAGE", "PARMETIS" );
  if( err == ZOLTAN_OK ) err = LB.Set_Param( "PARMETIS_METHOD", partitionMethod_ );
#endif

  //Setup Query Object
  CrsGraph_Transpose transposeTransform;
  Epetra_CrsGraph & TransGraph = transposeTransform( orig );
  ZoltanQuery Query( orig, &TransGraph );
  if( err == ZOLTAN_OK ) err = LB.Set_QueryObject( &Query );

  if( err != ZOLTAN_OK )
  { cout << "Setup of Zoltan Load Balancing Objects FAILED!\n"; exit(0); }

  //Generate Load Balance
  int changes;
  int num_gid_entries, num_lid_entries;
  int num_import;
  ZOLTAN_ID_PTR import_global_ids, import_local_ids;
  int * import_procs;
  int num_export;
  ZOLTAN_ID_PTR export_global_ids, export_local_ids;
  int * export_procs;

  orig.Comm().Barrier();
  err = LB.Balance( &changes,
                     &num_gid_entries, &num_lid_entries,
                     &num_import, &import_global_ids, &import_local_ids, &import_procs,
                     &num_export, &export_global_ids, &export_local_ids, &export_procs );
  LB.Evaluate( 1, 0, 0, 0, 0, 0, 0 );
  orig.Comm().Barrier();

  //Generate New Element List
  int numMyElements = orig.RowMap().NumMyElements();
  vector<int> elementList( numMyElements );
  orig.RowMap().MyGlobalElements( &elementList[0] );

  int newNumMyElements = numMyElements - num_export + num_import;
  vector<int> newElementList( newNumMyElements );

  set<int> gidSet;
  for( int i = 0; i < num_export; ++i )
    gidSet.insert( export_global_ids[i] );

  //Add unmoved indices to new list
  int loc = 0;
  for( int i = 0; i < numMyElements; ++i )
    if( !gidSet.count( elementList[i] ) )
      newElementList[loc++] = elementList[i];
  
  //Add imports to end of list
  for( int i = 0; i < num_import; ++i )
    newElementList[loc+i] = import_global_ids[i];

  //Free Zoltan Data
  if( err == ZOLTAN_OK )
    err = LB.Free_Data( &import_global_ids, &import_local_ids, &import_procs,
                         &export_global_ids, &export_local_ids, &export_procs );

  //Create Import Map
  NewRowMap_ = new Epetra_Map( orig.RowMap().NumGlobalElements(),
                                           newNumMyElements,
                                           &newElementList[0],
                                           orig.RowMap().IndexBase(),
                                           orig.RowMap().Comm() );

  //Create Importer
  Epetra_Import Importer( *NewRowMap_, orig.RowMap() );

  //Create New Graph
  Epetra_CrsGraph * NewGraph = new Epetra_CrsGraph( Copy, *NewRowMap_, 0 );
  NewGraph->Import( orig, Importer, Insert );
  NewGraph->FillComplete();

  Zoltan::LoadBalance LB2( 0, &dummy, &version );
  err = LB2.Create( dynamic_cast<const Epetra_MpiComm&>(orig.Comm()).Comm() );
  if( err == ZOLTAN_OK ) err = LB2.Set_Param( "LB_METHOD", "GRAPH" );
#ifdef HAVE_LIBPARMETIS
  if( err == ZOLTAN_OK ) err = LB2.Set_Param( "GRAPH_PACKAGE", "PARMETIS" );
  if( err == ZOLTAN_OK ) err = LB2.Set_Param( "PARMETIS_METHOD", partitionMethod_ );
#endif
  CrsGraph_Transpose transTrans;
  Epetra_CrsGraph & trans2 = transTrans( *NewGraph );
  ZoltanQuery query( *NewGraph, &trans2 );
  if( err == ZOLTAN_OK ) err = LB2.Set_QueryObject( &query );
  //err = LB2.Balance( &changes,
//.........这里部分代码省略.........

//==============================================================================
int checkSharedOwnership(Epetra_Comm& Comm, bool verbose) {
	// check to make sure each function returns 1 when it should
	// check to make sure each function doesn't return 1 when it shouldn't
	int ierr = 0;

	// initialize Map
	const int NumMyElements = 10;
	const int IndexBase = 0;
	Epetra_Map Map1((long long) -1, NumMyElements, IndexBase, Comm);
	// initialize Graphs
	const int NumIndicesPerRow = 5;
	Epetra_CrsGraph * SoleOwner = new Epetra_CrsGraph(Copy, Map1, Map1, NumIndicesPerRow);
	Epetra_CrsGraph SharedOrig(Copy, Map1, Map1, NumIndicesPerRow);
	Epetra_CrsGraph SharedOwner(SharedOrig);
	// arrays used by Insert & Remove
	Epetra_IntSerialDenseVector array1(2);
	array1[0] = NumIndicesPerRow / 2;
	array1[1] = array1[0] + 1;
	Epetra_LongLongSerialDenseVector array2(NumIndicesPerRow);
	for(int i = 0; i < NumIndicesPerRow; i++)
		array2[i] = i;
	// output variables (declaring them here lets us comment out indiv. tests)
	int soleOutput, sharedOutput;

	// InsertMyIndices
	if(verbose) cout << "InsertMyIndices..." << endl;
	soleOutput = SoleOwner->InsertMyIndices(0, 2, array1.Values());
	sharedOutput = SharedOwner.InsertMyIndices(0, 2, array1.Values());
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
	if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// RemoveMyIndices (#0)
	if(verbose) cout << "RemoveMyIndices(#0)..." << endl;
	soleOutput = SoleOwner->RemoveMyIndices(0);
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);

  EPETRA_TEST_ERR(!(SoleOwner->NumMyIndices(0)==0),ierr);
  if (ierr != 0) cout << "tests FAILED" << std::endl;

	soleOutput = SoleOwner->InsertMyIndices(0, 2, array1.Values());
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);

	// SortIndices
	//if(verbose) cout << "SortIndices..." << endl;
	//soleOutput = SoleOwner.SortIndices();
	//sharedOutput = SharedOwner.SortIndices();
	//EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	//EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
	//if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// RemoveRedundantIndices
	//if(verbose) cout << "RemoveRedundantIndices..." << endl;
	//SoleOwner.InsertGlobalIndices(0, 1, array1.Values());
	//SharedOwner.InsertGlobalIndices(0, 1, array1.Values());
	//soleOutput = SoleOwner.RemoveRedundantIndices();
	//sharedOutput = SharedOwner.RemoveRedundantIndices();
	//EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	//EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
	//if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// FillComplete (#1)
	if(verbose) cout << "FillComplete..." << endl;
	soleOutput = SoleOwner->FillComplete();
	sharedOutput = SharedOwner.FillComplete();
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
	if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// OptimizeStorage
	if(verbose) cout << "OptimizeStorage..." << endl;
	soleOutput = SoleOwner->OptimizeStorage();
	sharedOutput = SharedOwner.OptimizeStorage();
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == 0), ierr);
	if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// RemoveMyIndices (#1)
	if(verbose) cout << "RemoveMyIndices..." << endl;
	soleOutput = SoleOwner->RemoveMyIndices(0, 1, &array1[1]);
	sharedOutput = SharedOwner.RemoveMyIndices(0, 1, &array1[1]);
	EPETRA_TEST_ERR(!(soleOutput == -1), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == -1), ierr);
	if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// RemoveMyIndices (#2)
	if(verbose) cout << "RemoveMyIndices(#2)..." << endl;
	soleOutput = SoleOwner->RemoveMyIndices(0);
	sharedOutput = SharedOwner.RemoveMyIndices(0);
	EPETRA_TEST_ERR(!(soleOutput == -1), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == -1), ierr);
	if(verbose && ierr > 0) cout << "soleOutput = " << soleOutput << " sharedOutput = " << sharedOutput << endl;

	// FillComplete (#2)
	if(verbose) cout << "FillComplete(#2)..." << endl;
	soleOutput = SoleOwner->FillComplete(SoleOwner->DomainMap(), SoleOwner->RangeMap());
	sharedOutput = SharedOwner.FillComplete(SharedOwner.DomainMap(), SharedOwner.RangeMap());
	EPETRA_TEST_ERR(!(soleOutput == 0), ierr);
	EPETRA_TEST_ERR(!(sharedOutput == 1), ierr);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    while( !Finished )
    {
      DeeperFound = false;
      NarrowerFound = false;

      for( int i = 0; i < Leaves.size(); ++i )
      {

        BFT * TestBFT = new BFT( AdjList, Leaves[i], MinWidth, TooWide );

        if( TooWide )
          delete TestBFT;
        else
        {
          if( TestBFT->Width() < MinWidth ) MinWidth = TestBFT->Width();

          if( TestBFT->Depth() > Diameter )
          {
            delete BestBFT;
            Diameter = TestBFT->Depth();
            BestWidth = TestBFT->Width();
            BestBFT = TestBFT;
            DeeperFound = true;
            NarrowerFound = false;
          }
          else if( (TestBFT->Depth()==Diameter) && (TestBFT->Width()<BestWidth) )
          {
            delete BestBFT;
            BestWidth = TestBFT->Width();
            BestBFT = TestBFT;
            NarrowerFound = true;
          }
          else delete TestBFT;
        }
      }

      if( DeeperFound )
        BestBFT->NonNeighborLeaves( Leaves, AdjList, testLeafWidth_ );
      else if( NarrowerFound )
        Finished = true;
      else Finished = true;
    }
  }

  //std::cout << "\nSymmRCM:\n";
  //std::cout << "----------------------------\n";
  //std::cout << " Depth: " << BestBFT->Depth() << std::endl;
  //std::cout << " Width: " << BestBFT->Width() << std::endl;
  //std::cout << "----------------------------\n\n";

  std::vector<int> RCM;
  BestBFT->ReverseVector( RCM );
  for( int i = 0; i < NumNodes; ++i )
    RCM[i] = orig.RowMap().GID( RCM[i] );

  //Generate New Row Map
  RCMMap_ = new Epetra_Map( orig.RowMap().NumGlobalElements(),
                                        NumNodes,
                                        &RCM[0],
                                        orig.RowMap().IndexBase(),
                                        orig.RowMap().Comm() );

  //Generate New Col Map
  if( RCMMap_->DistributedGlobal() )
  {
    std::vector<int> colIndices = RCM;
    const Epetra_BlockMap & origColMap = orig.ColMap();

    if( origColMap.NumMyElements() > RCMMap_->NumMyElements() )
    {
      for( int i = RCMMap_->NumMyElements(); i < origColMap.NumMyElements(); ++i )
        colIndices.push_back( origColMap.GID(i) );
    }

    RCMColMap_ = new Epetra_Map( orig.ColMap().NumGlobalElements(),
                                 colIndices.size(),
                                 &colIndices[0],
                                 orig.ColMap().IndexBase(),
                                 orig.ColMap().Comm() );
  }
  else
    RCMColMap_ = RCMMap_;

  //Create New Graph
  Epetra_Import Importer( *RCMMap_, orig.RowMap() );
  Epetra_CrsGraph * RCMGraph = new Epetra_CrsGraph( Copy, *RCMMap_, *RCMColMap_, 0 );
  RCMGraph->Import( orig, Importer, Insert );
  RCMGraph->FillComplete();

/*
  std::cout << "origGraph\n";
  std::cout << orig;
  std::cout << "RCMGraph\n";
  std::cout << *RCMGraph;
*/

  newObj_ = RCMGraph;
  
  return *RCMGraph;
}