C++ Epetra_CrsGraph类代码示例

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

}

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

//EpetraCrsMatrix_To_TpetraCrsMatrix: copies Epetra_CrsMatrix to its analogous Tpetra_CrsMatrix
Teuchos::RCP<Tpetra_CrsMatrix> Petra::EpetraCrsMatrix_To_TpetraCrsMatrix(const Epetra_CrsMatrix& epetraCrsMatrix_,
                                                               const Teuchos::RCP<const Teuchos::Comm<int> >& commT_)
{
  //get row map of Epetra::CrsMatrix & convert to Tpetra::Map
  auto tpetraRowMap_ = EpetraMap_To_TpetraMap(epetraCrsMatrix_.RowMap(), commT_);

  //get col map of Epetra::CrsMatrix & convert to Tpetra::Map
  auto tpetraColMap_ = EpetraMap_To_TpetraMap(epetraCrsMatrix_.ColMap(), commT_);

  //get CrsGraph of Epetra::CrsMatrix & convert to Tpetra::CrsGraph
  const Epetra_CrsGraph epetraCrsGraph_ = epetraCrsMatrix_.Graph();
  std::size_t maxEntries = epetraCrsGraph_.GlobalMaxNumIndices();
  Teuchos::RCP<Tpetra_CrsGraph> tpetraCrsGraph_ = Teuchos::rcp(new Tpetra_CrsGraph(tpetraRowMap_, tpetraColMap_, maxEntries));

  for (LO i=0; i<epetraCrsGraph_.NumMyRows(); i++) {
     LO NumEntries; LO *Indices;
     epetraCrsGraph_.ExtractMyRowView(i, NumEntries, Indices);
     tpetraCrsGraph_->insertLocalIndices(i, NumEntries, Indices);
  }
  tpetraCrsGraph_->fillComplete();

  //convert Epetra::CrsMatrix to Tpetra::CrsMatrix, after creating Tpetra::CrsMatrix based on above Tpetra::CrsGraph
  Teuchos::RCP<Tpetra_CrsMatrix> tpetraCrsMatrix_ = Teuchos::rcp(new Tpetra_CrsMatrix(tpetraCrsGraph_));
  tpetraCrsMatrix_->setAllToScalar(0.0);

  for (LO i=0; i<epetraCrsMatrix_.NumMyRows(); i++) {
     LO NumEntries; LO *Indices; ST *Values;
     epetraCrsMatrix_.ExtractMyRowView(i, NumEntries, Values, Indices);
     tpetraCrsMatrix_->replaceLocalValues(i, NumEntries, Values, Indices);
  }
  tpetraCrsMatrix_->fillComplete();

  return tpetraCrsMatrix_;

}

//-----------------------------------------------------------------------------
// Function      : Indexor::setupAcceleratedMatrixIndexing
// Purpose       :
// Special Notes :
// Scope         : public
// Creator       : Rob Hoekstra, SNL, Parallel Computational Sciences
// Creation Date : 08/23/02
//-----------------------------------------------------------------------------
bool Indexor::setupAcceleratedMatrixIndexing( const std::string & graph_name )
{
  Epetra_CrsGraph * graph = 0;

  assert( pdsMgr_ != 0 );
  // Never, EVER do work inside an assert argument, or that work will not
  // be done when asserts are disabled.
  graph = pdsMgr_->getMatrixGraph( graph_name );
  assert( graph != 0 );

  int NumRows = graph->NumMyRows();
  matrixIndexMap_.clear();
  matrixIndexMap_.resize( NumRows );

  int NumElements;
  int * Elements;
  for( int i = 0; i < NumRows; ++i )
  {
    graph->ExtractMyRowView( i, NumElements, Elements );
    for( int j = 0; j < NumElements; ++j ) matrixIndexMap_[i][ Elements[j] ] = j;
  }

  accelMatrixIndex_ = true;

  return true;
}

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

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

void Epetra_OffsetIndex::GenerateLocalOffsets_( const Epetra_CrsGraph & SourceGraph,
                                                const Epetra_CrsGraph & TargetGraph,
                                                const int * PermuteLIDs )
{
  const int GlobalMaxNumSourceIndices = SourceGraph.GlobalMaxNumIndices();

  int NumSourceIndices;
  int_type * SourceIndices = 0;
  if( GlobalMaxNumSourceIndices>0 ) SourceIndices = new int_type[GlobalMaxNumSourceIndices];

  //setup Same Offsets
  SameOffsets_ = new int*[NumSame_];
  for( int i = 0; i < NumSame_; ++i ) SameOffsets_[i] = 0;

  for( int i = 0; i < NumSame_; ++i ) {
    int_type GID = (int_type) SourceGraph.GRID64(i);
    SourceGraph.ExtractGlobalRowCopy( GID,
                                      GlobalMaxNumSourceIndices,
                                      NumSourceIndices,
                                      SourceIndices );

    if( NumSourceIndices > 0 ) SameOffsets_[i] = new int[NumSourceIndices];

    int Loc = 0;
    int Start = 0;
    for( int j = 0; j < NumSourceIndices; ++j ) {
      Start = Loc;
      if( TargetGraph.FindGlobalIndexLoc(i,SourceIndices[j],Start,Loc) )
        SameOffsets_[i][j] = Loc;
      else
        SameOffsets_[i][j] = -1;
    }
  }

  //do also for permuted ids
  PermuteOffsets_ = new int*[NumPermute_];
  for( int i = 0; i < NumPermute_; ++i ) PermuteOffsets_[i] = 0;

  for( int i = 0; i < NumPermute_; ++i ) {
    int_type GID = (int_type) SourceGraph.GRID64(PermuteLIDs[i]);
    SourceGraph.ExtractGlobalRowCopy( GID,
                                      GlobalMaxNumSourceIndices,
                                      NumSourceIndices,
                                      SourceIndices );

    if( NumSourceIndices > 0 ) PermuteOffsets_[i] = new int[NumSourceIndices];

    int Loc = 0;
    int Start = 0;
    for( int j = 0; j < NumSourceIndices; ++j ) {
      Start = Loc;
      if( TargetGraph.FindGlobalIndexLoc(PermuteLIDs[i],SourceIndices[j],Start,Loc) )
        PermuteOffsets_[i][j] = Loc;
      else
        PermuteOffsets_[i][j] = -1;
    }
  }

  if( GlobalMaxNumSourceIndices>0 ) delete [] SourceIndices;
}

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

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

//-----------------------------------------------------------------------------
// Function      : Indexor::matrixGlobalToLocal
// Purpose       :
// Special Notes :
// Scope         : public
// Creator       : Rob Hoekstra, SNL, Parallel Computational Sciences
// Creation Date : 08/23/02
//-----------------------------------------------------------------------------
bool Indexor::matrixGlobalToLocal( const std::string & graph_name,
                                         const std::vector<int> & gids,
                                         std::vector< std::vector<int> > & stamp )
{
  Epetra_CrsGraph * graph = 0;

  assert( pdsMgr_ != 0 );
  // Never, EVER do work inside an assert argument, or that work will not
  // be done when asserts are disabled.
  graph = pdsMgr_->getMatrixGraph( graph_name );
  assert( graph != 0 );

  int numRows = stamp.size();

  int numElements;
  int * elements;

  if( accelMatrixIndex_ )
  {
    for( int i = 0; i < numRows; ++i )
    {
      int RowLID = graph->LRID(gids[i]);
      int NumCols = stamp[i].size();
      for( int j = 0; j < NumCols; ++j )
      {
        int lid = graph->LCID(stamp[i][j]);
        stamp[i][j] = matrixIndexMap_[RowLID][lid];
      }
    }
  }
  else
  {
    for( int i = 0; i < numRows; ++i )
    {
      graph->ExtractMyRowView( graph->LRID(gids[i]), numElements, elements );

      std::map<int,int> indexToOffsetMap;
      for( int j = 0; j < numElements; ++j ) indexToOffsetMap[ elements[j] ] = j;

      int numCols = stamp[i].size();
      for( int j = 0; j < numCols; ++j )
      {
        int lid = graph->LCID(stamp[i][j]);
//        assert( indexToOffsetMap.count(lid) );
        stamp[i][j] = indexToOffsetMap[lid];
      }
    }
  }

  return true;
}

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

  Teuchos::RCP<Epetra_CrsGraph> BlockAdjacencyGraph::compute( Epetra_CrsGraph& B, int nbrr, std::vector<int>&r, std::vector<double>& weights, bool verbose)
  {
    // Check if the graph is on one processor.
    int myMatProc = -1, matProc = -1;
    int myPID = B.Comm().MyPID();
    for (int proc=0; proc<B.Comm().NumProc(); proc++)
      {
	if (B.NumGlobalEntries() == B.NumMyEntries())
	  myMatProc = myPID;
      }
    B.Comm().MaxAll( &myMatProc, &matProc, 1 );
    
    if( matProc == -1)
      { cout << "FAIL for Global!  All CrsGraph entries must be on one processor!\n"; abort(); }
    
    int i= 0, j = 0, k, l = 0, p, pm, q = -1, ns;
    int tree_height;
    int error = -1;    /* error detected, possibly a problem with the input */
    int nrr;           /* number of rows in B */
    int nzM = 0;       /* number of edges in graph */
    int m = 0;         /* maximum number of nonzeros in any block row of B */
    int* colstack = 0; /* stack used to process each block row */
    int* bstree = 0;   /* binary search tree */
    std::vector<int> Mi, Mj, Mnum(nbrr+1,0);
    nrr = B.NumMyRows();
    if ( matProc == myPID && verbose )
      std::printf(" Matrix Size = %d      Number of Blocks = %d\n",nrr, nbrr);
    else
      nrr = -1;     /* Prevent processor from doing any computations */
    bstree = csr_bst(nbrr);  /* 0 : nbrr-1 */
    tree_height = ceil31log2(nbrr) + 1;
    error = -1;

    l = 0; j = 0; m = 0;
    for( i = 0; i < nrr; i++ ){
      if( i >= r[l+1] ){
	++l;                 /* new block row */
	m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */
	j = B.NumGlobalIndices(i);
      }else{
	j += B.NumGlobalIndices(i);
      }
    }
    /* one more time for the final block */
     m = EPETRA_MAX(m,j) ;   /* nonzeros in block row */

    colstack = (int*) malloc( EPETRA_MAX(m,1) * sizeof(int) );
    // The compressed graph is actually computed twice,
    // due to concerns about memory limitations.  First, 
    // without memory allocation, just nzM is computed.  
    // Next Mj is allocated. Then, the second time, the
    // arrays are actually populated.
    nzM = 0; q = -1; l = 0;
    int * indices;
    int numEntries;
    for( i = 0; i <= nrr; i++ ){
      if( i >= r[l+1] ){
	if( q > 0 ) std::qsort(colstack,q+1,sizeof(int),compare_ints); /* sort stack */
	if( q >= 0 ) ns = 1; /* l, colstack[0] M */
	for( j=1; j<=q ; j++ ){ /* delete copies */
	  if( colstack[j] > colstack[j-1] ) ++ns;
	}
	nzM += ns; /*M->p[l+1] = M->p[l] + ns;*/
	++l;
	q = -1;
      }
      if( i < nrr ){
	B.ExtractMyRowView( i, numEntries, indices );
	for( k = 0; k < numEntries; k++){
	  j = indices[k];  ns = 0; p = 0;
	  while( (r[bstree[p]] > j)  ||  (j >= r[bstree[p]+1])  ){
	    if( r[bstree[p]] > j){
	      p = 2*p+1;
	    }else{
	      if( r[bstree[p]+1] <= j) p = 2*p+2;
	    }
	    ++ns;
	    if( p > nbrr || ns > tree_height ) {
	      error = j;
	      std::printf("error: p %d  nbrr %d  ns %d %d\n",p,nbrr,ns,j); break;
	    }
	  }
	  colstack[++q] = bstree[p];
	}
	//if( error >-1 ){ std::printf("%d\n",error); break; }
        // p > nbrr is a fatal error that is ignored
      }
    }
    
    if ( matProc == myPID && verbose )
      std::printf("nzM =  %d \n", nzM );
    Mi.resize( nzM );
    Mj.resize( nzM );
    q = -1; l = 0; pm = -1;
    for( i = 0; i <= nrr; i++ ){
      if( i >= r[l+1] ){
	if( q > 0 ) std::qsort(colstack,q+1,sizeof(colstack[0]),compare_ints); /* sort stack */
	if( q >= 0 ){
	  Mi[++pm] = l;
	  Mj[pm] = colstack[0];
//.........这里部分代码省略.........

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