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

//==============================================================================
BlockCrsMatrix::BlockCrsMatrix(
        const Epetra_CrsGraph & BaseGraph,
        const Epetra_CrsGraph & LocalBlockGraph,
        const Epetra_Comm & GlobalComm  ) 
  : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseGraph, LocalBlockGraph, GlobalComm )) ),
    BaseGraph_( BaseGraph ),
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    RowStencil_int_( ),
    RowIndices_int_( ),
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    RowStencil_LL_( ),
    RowIndices_LL_( ),
#endif
    ROffset_(BlockUtility::CalculateOffset64(BaseGraph.RowMap())),
    COffset_(BlockUtility::CalculateOffset64(BaseGraph.ColMap()))
{
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && LocalBlockGraph.RowMap().GlobalIndicesInt())
    BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_int_, RowStencil_int_);
  else
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesLongLong() && LocalBlockGraph.RowMap().GlobalIndicesLongLong())
    BlockUtility::GenerateRowStencil(LocalBlockGraph, RowIndices_LL_, RowStencil_LL_);
  else
#endif
    throw "EpetraExt::BlockCrsMatrix::BlockCrsMatrix: Error, Global indices unknown.";
}

//==============================================================================
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);
}

//==============================================================================
// Epetra_OffsetIndex constructor from Importer
Epetra_OffsetIndex::Epetra_OffsetIndex( const Epetra_CrsGraph & SourceGraph,
                                        const Epetra_CrsGraph & TargetGraph,
                                        Epetra_Import & Importer )
  : Epetra_Object("Epetra::OffsetIndex"),
    NumSame_(0),
    SameOffsets_(0),
    NumPermute_(0),
    PermuteOffsets_(0),
    NumExport_(0),
    NumRemote_(0),
    RemoteOffsets_(0),
    DataOwned_(true)
{
  NumSame_ = Importer.NumSameIDs();

  NumPermute_ = Importer.NumPermuteIDs();
  int * PermuteLIDs = Importer.PermuteToLIDs();

  NumExport_ = Importer.NumExportIDs();
  int * ExportLIDs = Importer.ExportLIDs();

  NumRemote_ = Importer.NumRemoteIDs();
  int * RemoteLIDs = Importer.RemoteLIDs();

  if(!SourceGraph.RowMap().GlobalIndicesTypeMatch(TargetGraph.RowMap()))
     throw ReportError("Epetra_OffsetIndex::Epetra_OffsetIndex: SourceGraph and TargetGraph global indices type mismatch", -1);
  if(SourceGraph.RowMap().GlobalIndicesInt()) {
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
     GenerateLocalOffsets_<int>( SourceGraph, TargetGraph,
                            PermuteLIDs );

     GenerateRemoteOffsets_<int>( SourceGraph, TargetGraph,
                             ExportLIDs, RemoteLIDs,
                             Importer.Distributor() );
#else
    throw ReportError("Epetra_OffsetIndex::Epetra_OffsetIndex: ERROR, GlobalIndicesInt but no API for it.",-1);
#endif
  }
  else if(SourceGraph.RowMap().GlobalIndicesLongLong()) {
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
     GenerateLocalOffsets_<long long>( SourceGraph, TargetGraph,
                            PermuteLIDs );

     GenerateRemoteOffsets_<long long>( SourceGraph, TargetGraph,
                             ExportLIDs, RemoteLIDs,
                             Importer.Distributor() );
#else
    throw ReportError("Epetra_OffsetIndex::Epetra_OffsetIndex: ERROR, GlobalIndicesLongLong but no API for it.",-1);
#endif
  }
  else
     throw ReportError("Epetra_OffsetIndex::Epetra_OffsetIndex: SourceGraph global indices type unknown", -1);
}

void BlockUtility::TGenerateRowStencil(const Epetra_CrsGraph& LocalBlockGraph,
                                      std::vector<int_type> RowIndices,
                                      std::vector< std::vector<int_type> >& RowStencil)
{
  // Get row indices
  int NumMyRows = LocalBlockGraph.NumMyRows();
  RowIndices.resize(NumMyRows);
  const Epetra_BlockMap& RowMap = LocalBlockGraph.RowMap();
  RowMap.MyGlobalElements(&RowIndices[0]);

  // Get stencil
  RowStencil.resize(NumMyRows);
  if (LocalBlockGraph.IndicesAreGlobal()) {
    for (int i=0; i<NumMyRows; i++) {
      int_type Row = RowIndices[i];
      int NumCols = LocalBlockGraph.NumGlobalIndices(Row);
      RowStencil[i].resize(NumCols);
      LocalBlockGraph.ExtractGlobalRowCopy(Row, NumCols, NumCols,
                                           &RowStencil[i][0]);
      for (int k=0; k<NumCols; k++)
        RowStencil[i][k] -= Row;
    }
  }
  else {
    for (int i=0; i<NumMyRows; i++) {
      int NumCols = LocalBlockGraph.NumMyIndices(i);
	  std::vector<int> RowStencil_local(NumCols);
      RowStencil[i].resize(NumCols);
      LocalBlockGraph.ExtractMyRowCopy(i, NumCols, NumCols,
                                       &RowStencil_local[0]);
      for (int k=0; k<NumCols; k++)
        RowStencil[i][k] = (int_type) ((int) (LocalBlockGraph.GCID64(RowStencil_local[k]) - RowIndices[i]));
    }
  }
}

int compute_graph_metrics(const Epetra_CrsGraph &graph,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, int &numCuts, double &cutWgt, double &cutn, double &cutl)
{
  const Epetra_BlockMap & rmap = graph.RowMap();
  const Epetra_BlockMap & cmap = graph.ColMap();

  int maxEdges = cmap.NumMyElements();

  std::vector<std::vector<int> > myRows(rmap.NumMyElements());

  if (maxEdges > 0){
    int numEdges = 0;
    int *nborLID  = new int [maxEdges];
    for (int i=0; i<rmap.NumMyElements(); i++){

      graph.ExtractMyRowCopy(i, maxEdges, numEdges, nborLID);
      std::vector<int> cols(numEdges);
      for (int j=0; j<numEdges; j++){
        cols[j] = nborLID[j];
      }
      myRows[i] = cols;
    }
    delete [] nborLID;
  }
 
  return compute_graph_metrics(rmap, cmap, myRows, costs, myGoalWeight,
                               balance, numCuts, cutWgt, cutn, cutl);

}

void BlockUtility::GenerateRowStencil(const Epetra_CrsGraph& LocalBlockGraph,
                                      std::vector<long long> RowIndices,
                                      std::vector< std::vector<long long> >& RowStencil)
{
  if(LocalBlockGraph.RowMap().GlobalIndicesLongLong())
    BlockUtility::TGenerateRowStencil<long long>(LocalBlockGraph, RowIndices, RowStencil);
  else
    throw "EpetraExt::BlockUtility::GenerateRowStencil: Global Indices not long long.";
}

int compute_hypergraph_metrics(const Epetra_CrsGraph &graph,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, double &cutn, double &cutl)  // output
{
  return compute_hypergraph_metrics(graph.RowMap(), graph.ColMap(),
                                     graph.NumGlobalCols(), costs,
                                     myGoalWeight,
                                     balance, cutn, cutl);
}

int check(Epetra_CrsGraph& L, Epetra_CrsGraph& U, Ifpack_IlukGraph& LU,
          int NumGlobalRows1, int NumMyRows1, int LevelFill1, bool verbose) {
  using std::cout;
  using std::endl;

  int i, j;
  int NumIndices, * Indices;
  int NumIndices1, * Indices1;

  bool debug = true;

  Epetra_CrsGraph& L1 = LU.L_Graph();
  Epetra_CrsGraph& U1 = LU.U_Graph();

  // Test entries and count nonzeros

  int Nout = 0;

  for (i=0; i<LU.NumMyRows(); i++) {

    assert(L.ExtractMyRowView(i, NumIndices, Indices)==0);
    assert(L1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
    assert(NumIndices==NumIndices1);
    for (j=0; j<NumIndices1; j++) {
      if (debug &&(Indices[j]!=Indices1[j])) {
        int MyPID = L.RowMap().Comm().MyPID();
        cout << "Proc " << MyPID
             << " Local Row = " << i
             << "  L.Indices["<< j <<"]  = " << Indices[j]
             << " L1.Indices["<< j <<"] = " << Indices1[j] << endl;
      }
      assert(Indices[j]==Indices1[j]);
    }
    Nout += (NumIndices-NumIndices1);

    assert(U.ExtractMyRowView(i, NumIndices, Indices)==0);
    assert(U1.ExtractMyRowView(i, NumIndices1, Indices1)==0);
    assert(NumIndices==NumIndices1);
    for (j=0; j<NumIndices1; j++)  {
      if (debug &&(Indices[j]!=Indices1[j])) {
        int MyPID = L.RowMap().Comm().MyPID();
        cout << "Proc " << MyPID
             << " Local Row = " << i
             << "  U.Indices["<< j <<"]  = " << Indices[j]
             << " U1.Indices["<< j <<"] = " << Indices1[j] << endl;
      }
      assert(Indices[j]==Indices1[j]);
    }
    Nout += (NumIndices-NumIndices1);
  }

  // Test query functions

  int NumGlobalRows = LU.NumGlobalRows();
  if (verbose) cout << "\n\nNumber of Global Rows = " << NumGlobalRows << endl<< endl;

  assert(NumGlobalRows==NumGlobalRows1);

  int NumGlobalNonzeros = LU.NumGlobalNonzeros();
  if (verbose) cout << "\n\nNumber of Global Nonzero entries = "
                    << NumGlobalNonzeros << endl<< endl;

  int NoutG = 0;

  L.RowMap().Comm().SumAll(&Nout, &NoutG, 1);

  assert(NumGlobalNonzeros==L.NumGlobalNonzeros()+U.NumGlobalNonzeros()-NoutG);

  int NumMyRows = LU.NumMyRows();
  if (verbose) cout << "\n\nNumber of Rows = " << NumMyRows << endl<< endl;

  assert(NumMyRows==NumMyRows1);

  int NumMyNonzeros = LU.NumMyNonzeros();
  if (verbose) cout << "\n\nNumber of Nonzero entries = " << NumMyNonzeros << endl<< endl;

  assert(NumMyNonzeros==L.NumMyNonzeros()+U.NumMyNonzeros()-Nout);

  if (verbose) cout << "\n\nLU check OK" << endl<< endl;

  return(0);
}

void show_matrix(const char *txt, const Epetra_CrsGraph &graph, const Epetra_Comm &comm)
{
  int me = comm.MyPID();

  if (comm.NumProc() > 10){
    if (me == 0){
      std::cerr << txt << std::endl;
      std::cerr << "Printed matrix format only works for 10 or fewer processes" << std::endl;
    }
    return;
  }

  const Epetra_BlockMap &rowmap = graph.RowMap();
  const Epetra_BlockMap &colmap = graph.ColMap();

  int myRows = rowmap.NumMyElements();
  int numRows = graph.NumGlobalRows();
  int numCols = graph.NumGlobalCols();
  int base = rowmap.IndexBase();

  if ((numRows > 200) || (numCols > 500)){
    if (me == 0){
      std::cerr << txt << std::endl;
      std::cerr << "show_matrix: problem is too large to display" << std::endl;
    }
    return;
  }

  int *myA = new int [numRows * numCols];
  memset(myA, 0, sizeof(int) * numRows * numCols);

  int *myIndices;

  int *myRowGIDs = rowmap.MyGlobalElements();

  for (int i=0; i< myRows; i++){
    int myRowLID = rowmap.LID(myRowGIDs[i]);

    int numEntries = graph.NumMyIndices(myRowLID);

    if (numEntries > 0){
      int rc = graph.ExtractMyRowView(myRowLID, numEntries, myIndices);
      if (rc){
        std::cerr << txt << std::endl;
        std::cerr << "extract graph error" << std::endl;
        return;
      }

      int *row = myA + (numCols * (myRowGIDs[i] - base));

      for (int j=0; j < numEntries; j++){
        int gid = colmap.GID(myIndices[j]);
        row[gid-base] = me+1;
      }
    }
  }

  printMatrix(txt, myA, NULL, NULL, numRows, numCols, comm);

  delete [] myA;
}

//==============================================================================
int check(Epetra_CrsGraph& A, int NumMyRows1, long long NumGlobalRows1, int NumMyNonzeros1,
	  long long NumGlobalNonzeros1, long long* MyGlobalElements, bool verbose)
{
  (void)MyGlobalElements;
  int ierr = 0;
	int i;
	int j;
	int forierr = 0;
  int NumGlobalIndices;
  int NumMyIndices;
	int* MyViewIndices;
  int MaxNumIndices = A.MaxNumIndices();
  int* MyCopyIndices = new int[MaxNumIndices];
  long long* GlobalCopyIndices = new long long[MaxNumIndices];

  // Test query functions

  int NumMyRows = A.NumMyRows();
  if(verbose) cout << "Number of local Rows = " << NumMyRows << endl;

  EPETRA_TEST_ERR(!(NumMyRows==NumMyRows1),ierr);

  int NumMyNonzeros = A.NumMyNonzeros();
  if(verbose) cout << "Number of local Nonzero entries = " << NumMyNonzeros << endl;

  EPETRA_TEST_ERR(!(NumMyNonzeros==NumMyNonzeros1),ierr);

  long long NumGlobalRows = A.NumGlobalRows64();
  if(verbose) cout << "Number of global Rows = " << NumGlobalRows << endl;

  EPETRA_TEST_ERR(!(NumGlobalRows==NumGlobalRows1),ierr);

  long long NumGlobalNonzeros = A.NumGlobalNonzeros64();
  if(verbose) cout << "Number of global Nonzero entries = " << NumGlobalNonzeros << endl;

  EPETRA_TEST_ERR(!(NumGlobalNonzeros==NumGlobalNonzeros1),ierr);

  // GlobalRowView should be illegal (since we have local indices)

  EPETRA_TEST_ERR(!(A.ExtractGlobalRowView(A.RowMap().MaxMyGID64(), NumGlobalIndices, GlobalCopyIndices)==-2),ierr);

  // Other binary tests

  EPETRA_TEST_ERR(A.NoDiagonal(),ierr);
  EPETRA_TEST_ERR(!(A.Filled()),ierr);
  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MaxMyGID64())),ierr);
  EPETRA_TEST_ERR(!(A.MyGRID(A.RowMap().MinMyGID64())),ierr);
  EPETRA_TEST_ERR(A.MyGRID(1+A.RowMap().MaxMyGID64()),ierr);
  EPETRA_TEST_ERR(A.MyGRID(-1+A.RowMap().MinMyGID64()),ierr);
  EPETRA_TEST_ERR(!(A.MyLRID(0)),ierr);
  EPETRA_TEST_ERR(!(A.MyLRID(NumMyRows-1)),ierr);
  EPETRA_TEST_ERR(A.MyLRID(-1),ierr);
  EPETRA_TEST_ERR(A.MyLRID(NumMyRows),ierr);

  forierr = 0;
  for(i = 0; i < NumMyRows; i++) {
    long long Row = A.GRID64(i);
    A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyIndices);
    A.ExtractMyRowView(i, NumMyIndices, MyViewIndices);
    forierr += !(NumGlobalIndices==NumMyIndices);
    for(j = 1; j < NumMyIndices; j++) EPETRA_TEST_ERR(!(MyViewIndices[j-1]<MyViewIndices[j]),ierr);
    for(j = 0; j < NumGlobalIndices; j++) {
			forierr += !(GlobalCopyIndices[j]==A.GCID64(MyViewIndices[j]));
			forierr += !(A.LCID(GlobalCopyIndices[j])==MyViewIndices[j]);
    }
  }
  EPETRA_TEST_ERR(forierr,ierr);
  forierr = 0;
  for(i = 0; i < NumMyRows; i++) {
    long long Row = A.GRID64(i);
    A.ExtractGlobalRowCopy(Row, MaxNumIndices, NumGlobalIndices, GlobalCopyIndices);
    A.ExtractMyRowCopy(i, MaxNumIndices, NumMyIndices, MyCopyIndices);
    forierr += !(NumGlobalIndices==NumMyIndices);
    for(j = 1; j < NumMyIndices; j++)
			EPETRA_TEST_ERR(!(MyCopyIndices[j-1]<MyCopyIndices[j]),ierr);
    for(j = 0; j < NumGlobalIndices; j++) {
			forierr += !(GlobalCopyIndices[j]==A.GCID64(MyCopyIndices[j]));
			forierr += !(A.LCID(GlobalCopyIndices[j])==MyCopyIndices[j]);
    }

  }
  EPETRA_TEST_ERR(forierr,ierr);

  delete[] MyCopyIndices;
  delete[] GlobalCopyIndices;

  if(verbose) cout << "Rows sorted check OK" << endl;

  return(ierr);
}

void Epetra_OffsetIndex::GenerateRemoteOffsets_( const Epetra_CrsGraph & SourceGraph,
                                                 const Epetra_CrsGraph & TargetGraph,
                                                 const int * ExportLIDs,
                                                 const int * RemoteLIDs,
                                                 Epetra_Distributor & Distor )
{
  int numProcs = SourceGraph.RowMap().Comm().NumProc();
  if (numProcs < 2) {
    return;
  }

  const int GlobalMaxNumIndices = SourceGraph.GlobalMaxNumIndices();

  int NumIndices;
  /* "Indices" appears to be unused -- [email protected]
  int * Indices = 0;
  if( GlobalMaxNumIndices>0 ) Indices = new int[GlobalMaxNumIndices];
  */

  //Pack Source Rows
  int * Sizes = 0;
  if( NumExport_ > 0 ) Sizes = new int[NumExport_];
  int TotalSize = 0;
  for( int i = 0; i < NumExport_; ++i ) {
    Sizes[i] = SourceGraph.NumMyIndices(ExportLIDs[i]) + 1;
    TotalSize += Sizes[i];
  }

  int_type * SourceArray = new int_type[TotalSize+1];
  int Loc = 0;
  for( int i = 0; i < NumExport_; ++i ) {
    int_type GID = (int_type) SourceGraph.GRID64(ExportLIDs[i]);
    SourceArray[Loc] = Sizes[i]-1;
    SourceGraph.ExtractGlobalRowCopy( GID,
                                      GlobalMaxNumIndices,
                                      NumIndices,
                                      &(SourceArray[Loc+1]) );
    Loc += Sizes[i];
  }

  //Push to Target
  char * cRecvArray = 0;
  int_type * RecvArray = 0;
  int RecvArraySize = 0;

  Distor.Do( reinterpret_cast<char *>(SourceArray),
             (int)sizeof(int_type),
             Sizes,
             RecvArraySize,
             cRecvArray );
  RecvArray = reinterpret_cast<int_type*>(cRecvArray);

  //Construct RemoteOffsets
  if( NumRemote_ > 0 ) RemoteOffsets_ = new int*[NumRemote_];
  for( int i = 0; i < NumRemote_; ++i ) RemoteOffsets_[i] = 0;

  Loc = 0;
  for( int i = 0; i < NumRemote_; ++i ) {
    NumIndices = (int) RecvArray[Loc];
    RemoteOffsets_[i] = new int[NumIndices];
    ++Loc;
    int FLoc = 0;
    int Start = 0;
    for( int j = 0; j < NumIndices; ++j ) {
      Start = FLoc;
      if( TargetGraph.FindGlobalIndexLoc(RemoteLIDs[i],RecvArray[Loc],Start,FLoc) )
        RemoteOffsets_[i][j] = FLoc;
      else
        RemoteOffsets_[i][j] = -1;
      ++Loc;
    }
  }

  /* "Indices" appears to be unused -- [email protected]
  if( GlobalMaxNumIndices>0 ) delete [] Indices;
  */

  if( Sizes ) delete [] Sizes;
  if( SourceArray ) delete [] SourceArray;
  if( RecvArraySize ) delete [] cRecvArray;
}