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

//==============================================================================
int Ifpack_CrsRiluk::BlockGraph2PointGraph(const Epetra_CrsGraph & BG, Epetra_CrsGraph & PG, bool Upper) {

  if (!BG.IndicesAreLocal()) {EPETRA_CHK_ERR(-1);} // Must have done FillComplete on BG

  int * ColFirstPointInElementList = BG.RowMap().FirstPointInElementList();
  int * ColElementSizeList = BG.RowMap().ElementSizeList();
  if (BG.Importer()!=0) {
    ColFirstPointInElementList = BG.ImportMap().FirstPointInElementList();
    ColElementSizeList = BG.ImportMap().ElementSizeList();
  }

  int Length = (BG.MaxNumIndices()+1) * BG.ImportMap().MaxMyElementSize();
  vector<int> tmpIndices(Length);

  int BlockRow, BlockOffset, NumEntries;
  int NumBlockEntries;
  int * BlockIndices;

  int NumMyRows_tmp = PG.NumMyRows();

  for (int i=0; i<NumMyRows_tmp; i++) {
    EPETRA_CHK_ERR(BG.RowMap().FindLocalElementID(i, BlockRow, BlockOffset));
    EPETRA_CHK_ERR(BG.ExtractMyRowView(BlockRow, NumBlockEntries, BlockIndices));

    int * ptr = &tmpIndices[0]; // Set pointer to beginning of buffer

    int RowDim = BG.RowMap().ElementSize(BlockRow);
    NumEntries = 0;

    // This next line make sure that the off-diagonal entries in the block diagonal of the 
    // original block entry matrix are included in the nonzero pattern of the point graph
    if (Upper) {
      int jstart = i+1;
      int jstop = EPETRA_MIN(NumMyRows_tmp,i+RowDim-BlockOffset);
      for (int j= jstart; j< jstop; j++) {*ptr++ = j; NumEntries++;}
    }

    for (int j=0; j<NumBlockEntries; j++) {
      int ColDim = ColElementSizeList[BlockIndices[j]];
      NumEntries += ColDim;
      assert(NumEntries<=Length); // Sanity test
      int Index = ColFirstPointInElementList[BlockIndices[j]];
      for (int k=0; k < ColDim; k++) *ptr++ = Index++;
    }

    // This next line make sure that the off-diagonal entries in the block diagonal of the 
    // original block entry matrix are included in the nonzero pattern of the point graph
    if (!Upper) {
      int jstart = EPETRA_MAX(0,i-RowDim+1);
      int jstop = i;
      for (int j = jstart; j < jstop; j++) {*ptr++ = j; NumEntries++;}
    }

    EPETRA_CHK_ERR(PG.InsertMyIndices(i, NumEntries, &tmpIndices[0]));
  }

  SetAllocated(true);

  return(0);
}

/*----------------------------------------------------------------------*
 |  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;
}                                                     

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

//==============================================================================
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);
//.........这里部分代码省略.........