C++ Epetra_RowMatrix::RowMatrixRowMap方法代码示例

BlockCrsMatrix::BlockCrsMatrix(
        const Epetra_RowMatrix & BaseMatrix,
        const vector< vector<long long> > & RowStencil,
        const vector<long long> & RowIndices,
        const Epetra_Comm & GlobalComm  ) 
  : Epetra_CrsMatrix( Copy, *(BlockUtility::GenerateBlockGraph( BaseMatrix, RowStencil, RowIndices, GlobalComm )) ),
    BaseGraph_( Copy, BaseMatrix.RowMatrixRowMap(), 1 ), //Junk to satisfy constructor
    RowStencil_LL_( RowStencil ),
    RowIndices_LL_( RowIndices ),
    ROffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixRowMap())),
    COffset_(BlockUtility::CalculateOffset64(BaseMatrix.RowMatrixColMap()))
{
}

int RowMatrixToHandle(FILE * handle, const Epetra_RowMatrix & A) {
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  if(A.RowMatrixRowMap().GlobalIndicesInt()) {
    return RowMatrixToHandle<int>(handle, A);
  }
  else
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  if(A.RowMatrixRowMap().GlobalIndicesLongLong()) {
    return RowMatrixToHandle<long long>(handle, A);
  }
  else
#endif
    throw "EpetraExt::RowMatrixToHandle: GlobalIndices type unknown";
}

int compute_graph_metrics(const Epetra_RowMatrix &matrix,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, int &numCuts, double &cutWgt, double &cutn, double &cutl)
{

  const Epetra_Map &rmap = matrix.RowMatrixRowMap();
  const Epetra_Map &cmap = matrix.RowMatrixColMap();

  int maxEdges = matrix.MaxNumEntries();

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

  if (maxEdges > 0){
    int numEdges = 0;
    int *nborLID  = new int [maxEdges];
    double *tmp = new double [maxEdges];

    for (int i=0; i<rmap.NumMyElements(); i++){

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

}

void BlockCrsMatrix::LoadBlock(const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
{
  if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
	return TLoadBlock<int>(BaseMatrix, Row, Col);
  else
    throw "EpetraExt::BlockCrsMatrix::LoadBlock: Global indices not int";
}

void BlockCrsMatrix::SumIntoGlobalBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int Row, const int Col)
{
  if(Epetra_CrsMatrix::RowMatrixRowMap().GlobalIndicesInt() && BaseMatrix.RowMatrixRowMap().GlobalIndicesInt())
	return TSumIntoGlobalBlock<int>(alpha, BaseMatrix, Row, Col);
  else
    throw "EpetraExt::BlockCrsMatrix::SumIntoGlobalBlock: Global indices not int";
}

int writeRowMatrix(FILE * handle, const Epetra_RowMatrix & A) {

  long long numRows_LL = A.NumGlobalRows64();
  if(numRows_LL > std::numeric_limits<int>::max())
    throw "EpetraExt::writeRowMatrix: numRows_LL > std::numeric_limits<int>::max()";

  int numRows = static_cast<int>(numRows_LL);
  Epetra_Map rowMap = A.RowMatrixRowMap();
  Epetra_Map colMap = A.RowMatrixColMap();
  const Epetra_Comm & comm = rowMap.Comm();
  long long ioffset = 1 - rowMap.IndexBase64(); // Matlab indices start at 1
  long long joffset = 1 - colMap.IndexBase64(); // Matlab indices start at 1
  if (comm.MyPID()!=0) {
    if (A.NumMyRows()!=0) {EPETRA_CHK_ERR(-1);}
    if (A.NumMyCols()!=0) {EPETRA_CHK_ERR(-1);}
  }
  else {
    if (numRows!=A.NumMyRows()) {EPETRA_CHK_ERR(-1);}
    Epetra_SerialDenseVector values(A.MaxNumEntries());
    Epetra_IntSerialDenseVector indices(A.MaxNumEntries());
    for (int i=0; i<numRows; i++) {
      long long I = rowMap.GID64(i) + ioffset;
      int numEntries;
      if (A.ExtractMyRowCopy(i, values.Length(), numEntries, 
			     values.Values(), indices.Values())!=0) {EPETRA_CHK_ERR(-1);}
      for (int j=0; j<numEntries; j++) {
	long long J = colMap.GID64(indices[j]) + joffset;
	double val = values[j];
	fprintf(handle, "%lld %lld %22.16e\n", I, J, val);
      }
    }
  }
  return(0);
}

// ============================================================================
void EpetraExt::XMLWriter::
Write(const std::string& Label, const Epetra_RowMatrix& Matrix)
{
  TEUCHOS_TEST_FOR_EXCEPTION(IsOpen_ == false, std::logic_error,
                     "No file has been opened");

  int Rows = Matrix.NumGlobalRows();
  int Cols = Matrix.NumGlobalRows();
  int Nonzeros = Matrix.NumGlobalNonzeros();

  if (Comm_.MyPID() == 0)
  {
    std::ofstream of(FileName_.c_str(), std::ios::app);
    of << "<PointMatrix Label=\"" << Label << '"'
      << " Rows=\"" << Rows << '"'
      << " Columns=\"" << Cols<< '"'
      << " Nonzeros=\"" << Nonzeros << '"'
      << " Type=\"double\" StartingIndex=\"0\">" << std::endl;
  }

  int Length = Matrix.MaxNumEntries();
  std::vector<int> Indices(Length);
  std::vector<double> Values(Length);

  for (int iproc = 0; iproc < Comm_.NumProc(); iproc++)
  {
    if (iproc == Comm_.MyPID())
    {
      std::ofstream of(FileName_.c_str(), std::ios::app);
      of.precision(15);

      for (int i = 0; i < Matrix.NumMyRows(); ++i)
      {
        int NumMyEntries;
        Matrix.ExtractMyRowCopy(i, Length, NumMyEntries, &Values[0], &Indices[0]);

        int GRID = Matrix.RowMatrixRowMap().GID(i);

        for (int j = 0; j < NumMyEntries; ++j)
          of << GRID << " " << Matrix.RowMatrixColMap().GID(Indices[j])
             << " " << std::setiosflags(std::ios::scientific) << Values[j] << std::endl;
      }
      of.close();
    }
    Comm_.Barrier();
  }

  if (Comm_.MyPID() == 0)
  {
    std::ofstream of(FileName_.c_str(), std::ios::app);
    of << "</PointMatrix>" << std::endl;
    of.close();
  }
}

int RowMatrixToMatrixMarketFile( const char *filename, const Epetra_RowMatrix & A, 
				 const char * matrixName,
				 const char *matrixDescription, 
				 bool writeHeader) {
  long long M = A.NumGlobalRows64();
  long long N = A.NumGlobalCols64();
  long long nz = A.NumGlobalNonzeros64();

  FILE * handle = 0;

  if (A.RowMatrixRowMap().Comm().MyPID()==0) { // Only PE 0 does this section

    handle = fopen(filename,"w");
    if (!handle) {EPETRA_CHK_ERR(-1);}
    MM_typecode matcode;
    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_coordinate(&matcode);
    mm_set_real(&matcode);

    if (writeHeader==true) { // Only write header if requested (true by default)
    
      if (mm_write_banner(handle, matcode)!=0) {EPETRA_CHK_ERR(-1);}
      
      if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
      if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
      
      if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {EPETRA_CHK_ERR(-1);}
    }
  }
    
  if (RowMatrixToHandle(handle, A)!=0) {EPETRA_CHK_ERR(-1);}// Everybody calls this routine

  if (A.RowMatrixRowMap().Comm().MyPID()==0) // Only PE 0 opened a file
    if (fclose(handle)!=0) {EPETRA_CHK_ERR(-1);}
  return(0);
}

static int make_my_A(const Epetra_RowMatrix &matrix, int *myA, const Epetra_Comm &comm)
{
  int me = comm.MyPID();

  const Epetra_Map &rowmap = matrix.RowMatrixRowMap();
  const Epetra_Map &colmap = matrix.RowMatrixColMap();

  int myRows = matrix.NumMyRows();
  int numRows = matrix.NumGlobalRows();
  int numCols = matrix.NumGlobalCols();
  int base = rowmap.IndexBase();
  int maxRow = matrix.MaxNumEntries();

  memset(myA, 0, sizeof(int) * numRows * numCols);

  int *myIndices = new int [maxRow];
  double *tmp = new double [maxRow];

  int rowLen = 0;

  for (int i=0; i< myRows; i++){

    int rc = matrix.ExtractMyRowCopy(i, maxRow, rowLen, tmp, myIndices);

    if (rc){
      if (me == 0){
        std::cout << "Error in make_my_A" << std::endl;
      }
       return 1;
    }

    int *row = myA + (numCols * (rowmap.GID(i) - base));

    for (int j=0; j < rowLen; j++){

      int colGID = colmap.GID(myIndices[j]);
      
      row[colGID - base] = me + 1;
    }
  }

  if (maxRow){
    delete [] myIndices;
    delete [] tmp;
  }
  return 0;
}

int compute_hypergraph_metrics(const Epetra_RowMatrix &matrix,
            Isorropia::Epetra::CostDescriber &costs,
            double &myGoalWeight,
            double &balance, double &cutn, double &cutl)  // output
{
  const Epetra_BlockMap &rmap = 
    static_cast<const Epetra_BlockMap &>(matrix.RowMatrixRowMap());

  const Epetra_BlockMap &cmap = 
    static_cast<const Epetra_BlockMap &>(matrix.RowMatrixColMap());

  return compute_hypergraph_metrics(rmap, cmap,
                                     matrix.NumGlobalCols(),
                                     costs,
                                     myGoalWeight,
                                     balance, cutn, cutl);

}

int Amesos_Scalapack::NumericFactorization() {
  
  if( debug_ == 1 ) std::cout << "Entering `NumericFactorization()'" << std::endl;
  
  NumNumericFact_++;
  
  iam_ = Comm().MyPID();
  
  Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
  const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
  NumGlobalElements_ = OriginalMap.NumGlobalElements();

  NumGlobalNonzeros_ = RowMatrixA->NumGlobalNonzeros();
  
  RedistributeA();
  ConvertToScalapack();
  
  return PerformNumericFactorization( );
}

void BlockCrsMatrix::TSumIntoBlock(double alpha, const Epetra_RowMatrix & BaseMatrix, const int_type Row, const int_type Col)
{
  std::vector<int_type>& RowIndices_ = TRowIndices<int_type>();
  std::vector< std::vector<int_type> >& RowStencil_ = TRowStencil<int_type>();
  int_type RowOffset = RowIndices_[(std::size_t)Row] * ROffset_;
  int_type ColOffset = (RowIndices_[(std::size_t)Row] + RowStencil_[(std::size_t)Row][(std::size_t)Col]) * COffset_;

//  const Epetra_CrsGraph & BaseGraph = BaseMatrix.Graph();
  const Epetra_BlockMap & BaseMap = BaseMatrix.RowMatrixRowMap();
  const Epetra_BlockMap & BaseColMap = BaseMatrix.RowMatrixColMap();

  // This routine copies entries of a BaseMatrix into big  BlockCrsMatrix
  // It performs the following operation on the global IDs row-by-row
  // this->val[i+rowOffset][j+ColOffset] = BaseMatrix.val[i][j]

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

  for (int i=0; i<BaseMap.NumMyElements(); i++) {
    BaseMatrix.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices_local[0] );

    // Convert to BlockMatrix Global numbering scheme
    for( int l = 0; l < NumIndices; ++l ) {
       Indices_global[l] = ColOffset +  (int_type) BaseColMap.GID64(Indices_local[l]);
       Values[l] *= alpha;
    }

    int_type BaseRow = (int_type) BaseMap.GID64(i);
    ierr = this->SumIntoGlobalValues(BaseRow + RowOffset, NumIndices, &Values[0], &Indices_global[0]); 
    if (ierr != 0) std::cout << "WARNING BlockCrsMatrix::SumIntoBlock SumIntoGlobalValues err = " << ierr <<
	    "\n\t  Row " << BaseRow + RowOffset << "Col start" << Indices_global[0] << std::endl;

  }
}

//.........这里部分代码省略.........
      int MyRow;
      
      double *RowValues;
      int *ColIndices;
      int MaxNumEntries = FatOut_->MaxNumEntries();
      
      std::vector<int>ColIndicesV(MaxNumEntries);
      std::vector<double>RowValuesV(MaxNumEntries);
      
      int NumMyElements = FatOut_->NumMyRows() ; 
      for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
	EPETRA_CHK_ERR( FatOut_->
			ExtractMyRowView( MyRow, NzThisRow, RowValues, 
					  ColIndices ) != 0 ) ;
	//
	//  The following eight lines are just a sanity check on MyRow:
	//
	int MyGlobalRow =  FatOut_->GRID( MyRow );
	assert( MyGlobalRow%npcol_ == mypcol_ ) ;   // I should only own rows belonging to my processor column
	int MyTrueRow = MyGlobalRow/npcol_ ;  // This is the original row
	int UniformRows =  ( MyTrueRow / ( nprow_ * nb_ ) ) * nb_ ; 
	int AllExcessRows = MyTrueRow - UniformRows * nprow_ ; 
	int OurExcessRows =  AllExcessRows - ( myprow_ * nb_ ) ; 
	
	if (  MyRow != UniformRows + OurExcessRows ) {
	  std::cout << " iam _ = " << iam_ 
	       << " MyGlobalRow = " << MyGlobalRow 
	       << " MyTrueRow = " << MyTrueRow 
	       << " UniformRows = " << UniformRows 
	       << " AllExcessRows = " << AllExcessRows 
	       << " OurExcessRows = " << OurExcessRows 
	       << " MyRow = " << MyRow  << std::endl ;  
	}
	
	assert( OurExcessRows >= 0 &&  OurExcessRows < nb_ );
	assert( MyRow == UniformRows + OurExcessRows ) ; 
	
	for ( int j = 0; j < NzThisRow; j++ ) { 
	  assert(  FatOut_->RowMatrixColMap().GID( ColIndices[j] ) ==
		   FatOut_->GCID( ColIndices[j] ) );
	  
	  int MyGlobalCol =  FatOut_->GCID( ColIndices[j] );
	  assert( (MyGlobalCol/nb_)%npcol_ == mypcol_ ) ; 
	  int UniformCols =  ( MyGlobalCol / ( npcol_ * nb_ ) ) * nb_ ; 
	  int AllExcessCols = MyGlobalCol - UniformCols * npcol_ ; 
	  int OurExcessCols =  AllExcessCols - ( mypcol_ * nb_ ) ; 
	  assert( OurExcessCols >= 0 &&  OurExcessCols < nb_ );
	  int MyCol = UniformCols + OurExcessCols ; 
	  
	  DenseA_[ MyCol * lda_ + MyRow ] = RowValues[j] ; 
	}
      }
      
    } else { 
      
      int NumMyElements = ScaLAPACK1DMatrix_->NumMyRows() ; 
      
      assert( NumGlobalElements_ ==ScaLAPACK1DMatrix_->NumGlobalRows());
      assert( NumGlobalElements_ ==ScaLAPACK1DMatrix_->NumGlobalCols());
      DenseA_.resize( NumGlobalElements_ * NumMyElements ) ;
      for ( int i = 0 ; i < (int)DenseA_.size() ; i++ ) DenseA_[i] = 0 ; 
      
      int NzThisRow ;
      int MyRow;
      
      double *RowValues;
      int *ColIndices;
      int MaxNumEntries = ScaLAPACK1DMatrix_->MaxNumEntries();
      
      assert( DescA_[8] == lda_ ) ; //  Double check Lda
      std::vector<int>ColIndicesV(MaxNumEntries);
      std::vector<double>RowValuesV(MaxNumEntries);
      
      for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
	EPETRA_CHK_ERR( ScaLAPACK1DMatrix_->
			ExtractMyRowView( MyRow, NzThisRow, RowValues, 
					  ColIndices ) != 0 ) ;
	
	for ( int j = 0; j < NzThisRow; j++ ) { 
	  DenseA_[ ( ScaLAPACK1DMatrix_->RowMatrixColMap().GID( ColIndices[j] ) ) 
		   + MyRow * NumGlobalElements_ ] = RowValues[j] ; 
	}
      }
      //
      //  Create a map to allow us to redistribute the vectors X and B 
      //
      Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
      const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
      assert( NumGlobalElements_ == OriginalMap.NumGlobalElements() ) ;
      int NumMyElements_ = 0 ;
      if (iam_==0) NumMyElements_ = NumGlobalElements_;
      
      if (VectorMap_) { delete VectorMap_ ; VectorMap_ = 0 ; } 
      VectorMap_ = new Epetra_Map( NumGlobalElements_, NumMyElements_, 0, Comm() );
    }
  }
  ConTime_ += Time_->ElapsedTime();
  
  return 0;
}   

int Amesos_Scalapack::RedistributeA( ) {

  if( debug_ == 1 ) std::cout << "Entering `RedistributeA()'" << std::endl;

  Time_->ResetStartTime();
  
  Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
  EPETRA_CHK_ERR( RowMatrixA == 0 ) ; 

  const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
  int NumberOfProcesses = Comm().NumProc() ; 

  //
  //  Compute a uniform distribution as ScaLAPACK would want it
  //    MyFirstElement - The first element which this processor would have
  //    NumExpectedElemetns - The number of elements which this processor would have
  //

  int NumRows_ = RowMatrixA->NumGlobalRows() ; 
  int NumColumns_  = RowMatrixA->NumGlobalCols() ; 
  if ( MaxProcesses_ > 0 ) {
    NumberOfProcesses = EPETRA_MIN( NumberOfProcesses, MaxProcesses_ ) ; 
  }
  else {
    int ProcessNumHeuristic = (1+NumRows_/200)*(1+NumRows_/200);
    NumberOfProcesses = EPETRA_MIN( NumberOfProcesses,  ProcessNumHeuristic );
  }
  
  if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:171" << std::endl;
  //
  // Create the ScaLAPACK data distribution.
  // The TwoD data distribution is created in a completely different
  // manner and is not transposed (whereas the SaLAPACK 1D data
  // distribution was transposed) 
  //
  if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:163" << std::endl;
  Comm().Barrier(); 
  if ( TwoD_distribution_ ) { 
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:166" << std::endl;
    Comm().Barrier(); 
    npcol_ = EPETRA_MIN( NumberOfProcesses, 
			 EPETRA_MAX ( 2, (int) sqrt( NumberOfProcesses * 0.5 ) ) ) ; 
    nprow_ = NumberOfProcesses / npcol_ ;

    //
    //  Create the map for FatA - our first intermediate matrix
    //
    int NumMyElements = RowMatrixA->RowMatrixRowMap().NumMyElements() ;
    std::vector<int> MyGlobalElements( NumMyElements );
    RowMatrixA->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;

    int NumMyColumns = RowMatrixA->RowMatrixColMap().NumMyElements() ;
    std::vector<int> MyGlobalColumns( NumMyColumns );
    RowMatrixA->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;

    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:194" << std::endl;

    std::vector<int> MyFatElements( NumMyElements * npcol_ ); 
    
    for( int LocalRow=0; LocalRow<NumMyElements; LocalRow++ ) {
      for (int i = 0 ; i < npcol_; i++ ){
	MyFatElements[LocalRow*npcol_+i] = MyGlobalElements[LocalRow]*npcol_+i;
      } 
    }
    
    Epetra_Map FatInMap( npcol_*NumRows_, NumMyElements*npcol_, 
			 &MyFatElements[0], 0, Comm() ); 
    
    //
    //  Create FatIn, our first intermediate matrix
    //
    Epetra_CrsMatrix FatIn( Copy, FatInMap, 0 );
    
    
    std::vector<std::vector<int> > FatColumnIndices(npcol_,std::vector<int>(1));
    std::vector<std::vector<double> > FatMatrixValues(npcol_,std::vector<double>(1));
    std::vector<int> FatRowPtrs(npcol_);  // A FatRowPtrs[i] = the number 
    // of entries in local row LocalRow*npcol_ + i 
    
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:219" << std::endl;
    //
    mypcol_ = iam_%npcol_;
    myprow_ = (iam_/npcol_)%nprow_;
    if ( iam_ >= nprow_ * npcol_ ) {
      myprow_ = nprow_; 
      mypcol_ = npcol_; 
    }
    //  Each row is split into npcol_ rows, with each of the 
    //  new rows containing only those elements belonging to 
    //  its process column (in the ScaLAPACK 2D process grid)
    //
    int MaxNumIndices = RowMatrixA->MaxNumEntries();
    int NumIndices;
    std::vector<int> ColumnIndices(MaxNumIndices);
    std::vector<double> MatrixValues(MaxNumIndices); 
    
    if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:232 NumMyElements = " 
			    << NumMyElements 
			    << std::endl;
    
//.........这里部分代码省略.........

//=============================================================================
int Amesos_Dscpack::PerformSymbolicFactorization()
{
  ResetTimer(0);
  ResetTimer(1);

  MyPID_    = Comm().MyPID();
  NumProcs_ = Comm().NumProc();
  
  Epetra_RowMatrix *RowMatrixA = Problem_->GetMatrix();
  if (RowMatrixA == 0)
    AMESOS_CHK_ERR(-1);

  const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ;
  const Epetra_MpiComm& comm1   = dynamic_cast<const Epetra_MpiComm &> (Comm());
  int numrows                   = RowMatrixA->NumGlobalRows();
  int numentries                = RowMatrixA->NumGlobalNonzeros();

  Teuchos::RCP<Epetra_CrsGraph> Graph;

  Epetra_CrsMatrix* CastCrsMatrixA = 
    dynamic_cast<Epetra_CrsMatrix*>(RowMatrixA); 

  if (CastCrsMatrixA)
  {
    Graph = Teuchos::rcp(const_cast<Epetra_CrsGraph*>(&(CastCrsMatrixA->Graph())), false);
  }
  else
  {
    int MaxNumEntries = RowMatrixA->MaxNumEntries();
    Graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, OriginalMap, MaxNumEntries));

    std::vector<int>    Indices(MaxNumEntries);
    std::vector<double> Values(MaxNumEntries);

    for (int i = 0 ; i < RowMatrixA->NumMyRows() ; ++i)
    {
      int NumEntries;
      RowMatrixA->ExtractMyRowCopy(i, MaxNumEntries, NumEntries,
                                   &Values[0], &Indices[0]);

      for (int j = 0 ; j < NumEntries ; ++j)
        Indices[j] = RowMatrixA->RowMatrixColMap().GID(Indices[j]);

      int GlobalRow = RowMatrixA->RowMatrixRowMap().GID(i);
      Graph->InsertGlobalIndices(GlobalRow, NumEntries, &Indices[0]);
    }

    Graph->FillComplete();
  }

  //
  //  Create a replicated map and graph 
  //
  std::vector<int> AllIDs( numrows ) ; 
  for ( int i = 0; i < numrows ; i++ ) AllIDs[i] = i ; 

  Epetra_Map      ReplicatedMap( -1, numrows, &AllIDs[0], 0, Comm());
  Epetra_Import   ReplicatedImporter(ReplicatedMap, OriginalMap);
  Epetra_CrsGraph ReplicatedGraph( Copy, ReplicatedMap, 0 ); 

  AMESOS_CHK_ERR(ReplicatedGraph.Import(*Graph, ReplicatedImporter, Insert));
  AMESOS_CHK_ERR(ReplicatedGraph.FillComplete());

  //
  //  Convert the matrix to Ap, Ai
  //
  std::vector <int> Replicates(numrows);
  std::vector <int> Ap(numrows + 1);
  std::vector <int> Ai(EPETRA_MAX(numrows, numentries));

  for( int i = 0 ; i < numrows; i++ ) Replicates[i] = 1; 
  
  int NumEntriesPerRow ;
  int *ColIndices = 0 ;
  int Ai_index = 0 ; 
  for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
    AMESOS_CHK_ERR( ReplicatedGraph.ExtractMyRowView( MyRow, NumEntriesPerRow, ColIndices ) );
    Ap[MyRow] = Ai_index ; 
    for ( int j = 0; j < NumEntriesPerRow; j++ ) { 
      Ai[Ai_index] = ColIndices[j] ; 
      Ai_index++;
    }
  }
  assert( Ai_index == numentries ) ; 
  Ap[ numrows ] = Ai_index ; 
  
  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);

  ResetTimer(0);

  //
  //  Call Dscpack Symbolic Factorization
  //  
  int OrderCode = 2;
  std::vector<double> MyANonZ;
  
  NumLocalNonz = 0 ; 
  GlobalStructNewColNum = 0 ; 
  GlobalStructNewNum = 0 ;  
//.........这里部分代码省略.........