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

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

}

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

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

//============================================================================
// very simple debugging function that prints on screen the structure
// of the input matrix.
void Ifpack_PrintSparsity_Simple(const Epetra_RowMatrix& A)
{
  int MaxEntries = A.MaxNumEntries();
  std::vector<int> Indices(MaxEntries);
  std::vector<double> Values(MaxEntries);
  std::vector<bool> FullRow(A.NumMyRows());

  cout << "+-";
  for (int j = 0 ; j < A.NumMyRows() ; ++j)
    cout << '-';
  cout << "-+" << endl;

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

    int Length;
    A.ExtractMyRowCopy(i,MaxEntries,Length,
                       &Values[0], &Indices[0]);

    for (int j = 0 ; j < A.NumMyRows() ; ++j)
      FullRow[j] = false;

    for (int j = 0 ; j < Length ; ++j) {
      FullRow[Indices[j]] = true;
    }

    cout << "| ";
    for (int j = 0 ; j < A.NumMyRows() ; ++j) {
      if (FullRow[j])
        cout << '*';
      else
        cout << ' ';
    }
    cout << " |" << endl;
  }

  cout << "+-";
  for (int j = 0 ; j < A.NumMyRows() ; ++j)
    cout << '-';
  cout << "-+" << endl << endl;

}

double Ifpack_FrobeniusNorm(const Epetra_RowMatrix& A)
{
  double MyNorm = 0.0, GlobalNorm;

  std::vector<int> colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

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

    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    for (int j = 0 ; j < Nnz ; ++j) {
      MyNorm += colVal[j] * colVal[j];
    }
  }

  A.Comm().SumAll(&MyNorm,&GlobalNorm,1);

  return(sqrt(GlobalNorm));
}

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

//.........这里部分代码省略.........
    Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
    Epetra_MultiVector YA2(YA1);
    Epetra_MultiVector YB1(YA1);
    Epetra_MultiVector YB2(YA1);
    X.Random();

    bool transA = false;
    A.SetUseTranspose(transA);
    B.SetUseTranspose(transA);
    A.Apply(X,YA1);
    A.Multiply(transA, X, YA2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
    B.Apply(X,YB1);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
    B.Multiply(transA, X, YB2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);

  }
  {// transpose case
    Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
    Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
    Epetra_MultiVector YA2(YA1);
    Epetra_MultiVector YB1(YA1);
    Epetra_MultiVector YB2(YA1);
    X.Random();

    bool transA = true;
    A.SetUseTranspose(transA);
    B.SetUseTranspose(transA);
    A.Apply(X,YA1);
    A.Multiply(transA, X, YA2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
    B.Apply(X,YB1);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
    B.Multiply(transA, X,YB2);
    EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);

  }

  Epetra_Vector diagA(A.RowMatrixRowMap());
  EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
  Epetra_Vector diagB(B.RowMatrixRowMap());
  EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
  EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);

  Epetra_Vector rowA(A.RowMatrixRowMap());
  EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
  Epetra_Vector rowB(B.RowMatrixRowMap());
  EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
  EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);

  Epetra_Vector colA(A.RowMatrixColMap());
  EPETRA_TEST_ERR(A.InvColSums(colA),ierr);
  Epetra_Vector colB(B.RowMatrixColMap());
  EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
  EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr);

  EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
  EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);

  EPETRA_TEST_ERR(A.RightScale(colA),ierr);
  EPETRA_TEST_ERR(B.RightScale(colB),ierr);


  EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
  EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);


  EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf after scaling", verbose), ierr);
  EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne after scaling", verbose),ierr);

  vector<double> valuesA(A.MaxNumEntries());
  vector<int> indicesA(A.MaxNumEntries());
  vector<double> valuesB(B.MaxNumEntries());
  vector<int> indicesB(B.MaxNumEntries());
  return(0);
  for (int i=0; i<A.NumMyRows(); i++) {
    int nA, nB;
    EPETRA_TEST_ERR(A.ExtractMyRowCopy(i, A.MaxNumEntries(), nA, &valuesA[0], &indicesA[0]),ierr);
    EPETRA_TEST_ERR(B.ExtractMyRowCopy(i, B.MaxNumEntries(), nB, &valuesB[0], &indicesB[0]),ierr);
    EPETRA_TEST_ERR(!nA==nB,ierr);
    for (int j=0; j<nA; j++) {
      double curVal = valuesA[j];
      int curIndex = indicesA[j];
      bool notfound = true;
      int jj = 0;
      while (notfound && jj< nB) {
	if (!checkValues(curVal, valuesB[jj])) notfound = false;
	jj++;
      }
      EPETRA_TEST_ERR(notfound, ierr);
      vector<int>::iterator p = find(indicesB.begin(),indicesB.end(),curIndex);  // find curIndex in indicesB
      EPETRA_TEST_ERR(p==indicesB.end(), ierr);
    }

  }
  if (verbose) cout << "RowMatrix Methods check OK" << endl;

  return (ierr);
}

int Ifpack_CrsRiluk::InitAllValues(const Epetra_RowMatrix & OverlapA, int MaxNumEntries) {

  int ierr = 0;
  int i, j;
  int NumIn, NumL, NumU;
  bool DiagFound;
  int NumNonzeroDiags = 0;


  vector<int> InI(MaxNumEntries); // Allocate temp space
  vector<int> LI(MaxNumEntries);
  vector<int> UI(MaxNumEntries);
  vector<double> InV(MaxNumEntries);
  vector<double> LV(MaxNumEntries);
  vector<double> UV(MaxNumEntries);

  bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced

  if (ReplaceValues) {
    L_->PutScalar(0.0); // Zero out L and U matrices
    U_->PutScalar(0.0);
  }

  D_->PutScalar(0.0); // Set diagonal values to zero
  double *DV;
  EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
    

  // First we copy the user's matrix into L and U, regardless of fill level

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

    EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0])); // Get Values and Indices
    
    // Split into L and U (we don't assume that indices are ordered).
    
    NumL = 0; 
    NumU = 0; 
    DiagFound = false;
    
    for (j=0; j< NumIn; j++) {
      int k = InI[j];

      if (k==i) {
	DiagFound = true;
	DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_
      }

      else if (k < 0) {EPETRA_CHK_ERR(-1);} // Out of range

      else if (k < i) {
	LI[NumL] = k;
	LV[NumL] = InV[j];
	NumL++;
      }
      else if (k<NumMyRows()) {
	UI[NumU] = k;
	UV[NumU] = InV[j];
	NumU++;
      }
    }
    
    // Check in things for this row of L and U

    if (DiagFound) NumNonzeroDiags++;
    else DV[i] = Athresh_;

    if (NumL) {
      if (ReplaceValues) {
	EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));
      }
      else {
	EPETRA_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));
      }
    }

    if (NumU) {
      if (ReplaceValues) {
	EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));
      }
      else {
	EPETRA_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));
      }
    }
    
  }

  if (!ReplaceValues) {
    // The domain of L and the range of U are exactly their own row maps (there is no communication).
    // The domain of U and the range of L must be the same as those of the original matrix,
    // However if the original matrix is a VbrMatrix, these two latter maps are translation from
    // a block map to a point map.
    EPETRA_CHK_ERR(L_->FillComplete(L_->RowMatrixColMap(), *L_RangeMap_));
    EPETRA_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));
  }

  // At this point L and U have the values of A in the structure of L and U, and diagonal vector D

  SetValuesInitialized(true);
  SetFactored(false);
//.........这里部分代码省略.........

int DoCopyRowMatrix(mxArray* matlabA, int& valueCount, const Epetra_RowMatrix& A) {
  //cout << "doing DoCopyRowMatrix\n";
  int ierr = 0;
  int numRows = A.NumGlobalRows();
  //cout << "numRows: " << numRows << "\n";
  Epetra_Map rowMap = A.RowMatrixRowMap();
  Epetra_Map colMap = A.RowMatrixColMap();
  int minAllGID = rowMap.MinAllGID();

  const Epetra_Comm & comm = rowMap.Comm();
  //cout << "did global setup\n";
  if (comm.MyPID()!=0) {
    if (A.NumMyRows()!=0) ierr = -1;
    if (A.NumMyCols()!=0) ierr = -1;
  }
  else {
	// declare and get initial values of all matlabA pointers
	double* matlabAvaluesPtr = mxGetPr(matlabA);
	int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
	int* matlabArowIndicesPtr = mxGetIr(matlabA);

	// set all matlabA pointers to the proper offset
	matlabAvaluesPtr += valueCount;
	matlabArowIndicesPtr += valueCount;

    if (numRows!=A.NumMyRows()) ierr = -1;
    Epetra_SerialDenseVector values(A.MaxNumEntries());
    Epetra_IntSerialDenseVector indices(A.MaxNumEntries());
    //cout << "did proc0 setup\n";
    for (int i=0; i<numRows; i++) {
	  //cout << "extracting a row\n";
	  int I = rowMap.GID(i);
      int numEntries = 0;
      if (A.ExtractMyRowCopy(i, values.Length(), numEntries, 
	  		     values.Values(), indices.Values())) return(-1);
	  matlabAcolumnIndicesPtr[I - minAllGID] = valueCount;  // set the starting index of column I
	  double* serialValuesPtr = values.Values();
      for (int j=0; j<numEntries; j++) {
		int J = colMap.GID(indices[j]);
		*matlabAvaluesPtr = *serialValuesPtr++;
		*matlabArowIndicesPtr = J;
		// increment matlabA pointers
		matlabAvaluesPtr++;
		matlabArowIndicesPtr++;
		valueCount++;
      }
    }
    //cout << "proc0 row extraction for this chunck is done\n";
  }

/*
  if (comm.MyPID() == 0) {
  cout << "printing matlabA pointers\n";
	double* matlabAvaluesPtr = mxGetPr(matlabA);
	int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
	int* matlabArowIndicesPtr = mxGetIr(matlabA);
  for(int i=0; i < numRows; i++) {
	for(int j=0; j < A.MaxNumEntries(); j++) {
	  cout << "*matlabAvaluesPtr: " << *matlabAvaluesPtr++ << " *matlabAcolumnIndicesPtr: " << *matlabAcolumnIndicesPtr++ << " *matlabArowIndicesPtr" << *matlabArowIndicesPtr++ << "\n";
	}
  }
  
  cout << "done printing matlabA pointers\n";
  }
  */
  
  int ierrGlobal;
  comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
  return(ierrGlobal);
}

//.........这里部分代码省略.........
    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;
    
    nb_ = grid_nb_;
    
    for( int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
      
      RowMatrixA->ExtractMyRowCopy( LocalRow, 
				    MaxNumIndices,
				    NumIndices, 
				    &MatrixValues[0],
				    &ColumnIndices[0] );
      
      for (int i=0; i<npcol_; i++ )  FatRowPtrs[i] = 0 ; 

      //
      //  Deal the individual matrix entries out to the row owned by
      //  the process to which this matrix entry will belong.
      //
      for( int i=0 ; i<NumIndices ; ++i ) {
	int GlobalCol = MyGlobalColumns[ ColumnIndices[i] ];
	int pcol_i = pcolnum( GlobalCol, nb_, npcol_ ) ;
	if ( FatRowPtrs[ pcol_i ]+1 >= FatColumnIndices[ pcol_i ].size() ) {
	  FatColumnIndices[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
	  FatMatrixValues[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
	}
	FatColumnIndices[pcol_i][FatRowPtrs[pcol_i]] = GlobalCol ;
	FatMatrixValues[pcol_i][FatRowPtrs[pcol_i]] = MatrixValues[i];
	
	FatRowPtrs[ pcol_i ]++;
      }
      
      //
      //  Insert each of the npcol_ rows individually
      //
      for ( int pcol_i = 0 ; pcol_i < npcol_ ; pcol_i++ ) { 
	FatIn.InsertGlobalValues( MyGlobalElements[LocalRow]*npcol_ + pcol_i, 
				  FatRowPtrs[ pcol_i ],
				  &FatMatrixValues[ pcol_i ][0], 
				  &FatColumnIndices[ pcol_i ][0] );

//=============================================================================
int Amesos_Mumps::ConvertToTriplet(const bool OnlyValues)
{

  Epetra_RowMatrix* ptr;
  if (Comm().NumProc() == 1)
    ptr = &Matrix();
  else {
    ptr = &RedistrMatrix(true);
  }

  ResetTimer();
  
#ifdef EXTRA_DEBUG_INFO
  Epetra_CrsMatrix* Eptr = dynamic_cast<Epetra_CrsMatrix*>( ptr );
  if ( ptr->NumGlobalNonzeros() < 300 ) SetICNTL(4,3 );  // Enable more debug info for small matrices
  if ( ptr->NumGlobalNonzeros() < 42 && Eptr ) { 
      std::cout << " Matrix = " << std::endl ; 
      Eptr->Print( std::cout ) ; 
  } else {
      assert( Eptr );
  }
#endif

  Row.resize(ptr->NumMyNonzeros());
  Col.resize(ptr->NumMyNonzeros());
  Val.resize(ptr->NumMyNonzeros());

  int MaxNumEntries = ptr->MaxNumEntries();
  std::vector<int> Indices;
  std::vector<double> Values;
  Indices.resize(MaxNumEntries);
  Values.resize(MaxNumEntries);

  int count = 0;

  for (int i = 0; i < ptr->NumMyRows() ; ++i) {

    int GlobalRow = ptr->RowMatrixRowMap().GID(i);

    int NumEntries = 0;
    int ierr;
    ierr = ptr->ExtractMyRowCopy(i, MaxNumEntries,
				   NumEntries, &Values[0],
				   &Indices[0]);
    AMESOS_CHK_ERR(ierr);

    for (int j = 0 ; j < NumEntries ; ++j) {
      if (OnlyValues == false) {
	Row[count] = GlobalRow + 1;
	Col[count] = ptr->RowMatrixColMap().GID(Indices[j]) + 1;
      }
      
      // MS // Added on 15-Mar-05.
      if (AddToDiag_ && Indices[j] == i)
        Values[j] += AddToDiag_;

      Val[count] = Values[j];
      count++;
    }
  }

  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
  
  assert (count <= ptr->NumMyNonzeros());

  return(0);
}

int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
                   const int NumPDEEqns)
{

  int NumMyRows = A.NumMyRows();
  long long NumGlobalRows = A.NumGlobalRows64();
  long long NumGlobalCols = A.NumGlobalCols64();
  long long MyBandwidth = 0, GlobalBandwidth;
  long long MyLowerNonzeros = 0, MyUpperNonzeros = 0;
  long long GlobalLowerNonzeros, GlobalUpperNonzeros;
  long long MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
  long long MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
  double MyMin, MyAvg, MyMax;
  double GlobalMin, GlobalAvg, GlobalMax;
  long long GlobalStorage;

  bool verbose = (A.Comm().MyPID() == 0);

  GlobalStorage = sizeof(int*) * NumGlobalRows +
    sizeof(int) * A.NumGlobalNonzeros64() +
    sizeof(double) * A.NumGlobalNonzeros64();

  if (verbose) {
    print();
    Ifpack_PrintLine();
    print<const char*>("Label", A.Label());
    print<long long>("Global rows", NumGlobalRows);
    print<long long>("Global columns", NumGlobalCols);
    print<long long>("Stored nonzeros", A.NumGlobalNonzeros64());
    print<long long>("Nonzeros / row", A.NumGlobalNonzeros64() / NumGlobalRows);
    print<double>("Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
  }

  long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
  long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
  long long NumMyDirichletRows = 0, NumGlobalDirichletRows;

  std::vector<int> colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  Epetra_Vector Diag(A.RowMatrixRowMap());
  Epetra_Vector RowSum(A.RowMatrixRowMap());
  Diag.PutScalar(0.0);
  RowSum.PutScalar(0.0);

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

    long long GRID = A.RowMatrixRowMap().GID64(i);
    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    if (Nnz == 0)
      NumMyEmptyRows++;

    if (Nnz == 1)
      NumMyDirichletRows++;

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

      double v = colVal[j];
      if (v < 0) v = -v;
      if (colVal[j] != 0.0)
        NumMyActualNonzeros++;

      long long GCID = A.RowMatrixColMap().GID64(colInd[j]);

      if (GCID != GRID)
        RowSum[i] += v;
      else
        Diag[i] = v;

      if (GCID < GRID)
        MyLowerNonzeros++;
      else if (GCID > GRID)
        MyUpperNonzeros++;
      long long b = GCID - GRID;
      if (b < 0) b = -b;
      if (b > MyBandwidth)
        MyBandwidth = b;
    }

    if (Diag[i] > RowSum[i])
      MyDiagonallyDominant++;

    if (Diag[i] >= RowSum[i])
      MyWeaklyDiagonallyDominant++;

    RowSum[i] += Diag[i];
  }

  // ======================== //
  // summing up global values //
  // ======================== //

  A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
  A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
  A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
  A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
  A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    fprintf(fp,"%s","/pnum { 72 div 2.54 mul 20 string ");
    fprintf(fp,"%s","cvs print ( ) print} def\n");
    fprintf(fp,"%s","/Cshow {dup stringwidth pop -2 div 0 rmoveto show} def\n");

    /* we leave margins etc. in cm so it is easy to modify them if
       needed by editing the output file */
    fprintf(fp,"%s","gsave\n");
    if (ltit != 0) {
      fprintf(fp,"/Helvetica findfont %e cm scalefont setfont\n",
              fnstit);
      fprintf(fp,"%f cm %f cm moveto\n",
              xtit,ytit);
      fprintf(fp,"(%s) Cshow\n", title);
      fprintf(fp,"%f cm %f cm translate\n",
              lrmrgn,botmrgn);
    }
    fprintf(fp,"%f cm %d div dup scale \n",
            siz, (int) m);
    /* draw a frame around the matrix */

    fprintf(fp,"%f setlinewidth\n",
            frlw);
    fprintf(fp,"%s","newpath\n");
    fprintf(fp,"%s","0 0 moveto ");
    if (square) {
      printf("------------------- %d\n", (int) m);
      fprintf(fp,"%d %d lineto\n",
              (int) m, 0);
      fprintf(fp,"%d %d lineto\n",
              (int) m, (int) m);
      fprintf(fp,"%d %d lineto\n",
              0, (int) m);
    }
    else {
      fprintf(fp,"%d %d lineto\n",
              (int) nc, 0);
      fprintf(fp,"%d %d lineto\n",
              (int) nc, (int) nr);
      fprintf(fp,"%d %d lineto\n",
              0, (int) nr);
    }
    fprintf(fp,"%s","closepath stroke\n");

    /* plotting loop */

    fprintf(fp,"%s","1 1 translate\n");
    fprintf(fp,"%s","0.8 setlinewidth\n");
    fprintf(fp,"%s","/p {moveto 0 -.40 rmoveto \n");
    fprintf(fp,"%s","           0  .80 rlineto stroke} def\n");

    fclose(fp);
  }

  int MaxEntries = A.MaxNumEntries();
  std::vector<int> Indices(MaxEntries);
  std::vector<double> Values(MaxEntries);

  for (int pid = 0 ; pid < NumProc ; ++pid) {

    if (pid == MyPID) {

      fp = fopen(FileName,"a");
      if( fp == NULL ) {
        fprintf(stderr,"%s","ERROR\n");
        exit(EXIT_FAILURE);
      }

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

        if (i % NumPDEEqns) continue;

        int Nnz;
        A.ExtractMyRowCopy(i,MaxEntries,Nnz,&Values[0],&Indices[0]);

        long long grow = A.RowMatrixRowMap().GID64(i);

        for (int j = 0 ; j < Nnz ; ++j) {
          int col = Indices[j];
          if (col % NumPDEEqns == 0) {
            long long gcol = A.RowMatrixColMap().GID64(Indices[j]);
            grow /= NumPDEEqns;
            gcol /= NumPDEEqns;
            fprintf(fp,"%lld %lld p\n",
                    gcol, NumGlobalRows - grow - 1);
          }
        }
      }

      fprintf(fp,"%s","%end of data for this process\n");

      if( pid == NumProc - 1 )
        fprintf(fp,"%s","showpage\n");

      fclose(fp);
    }
    Comm.Barrier();
  }

  return(0);
}