C++ Epetra_RowMatrix类代码示例

void show_matrix(const char *txt, const Epetra_RowMatrix &matrix, const Epetra_Comm &comm)
{
  int me = comm.MyPID();
  if (comm.NumProc() > 10){
    if (me == 0){
      std::cout << txt << std::endl;
      std::cout << "Printed matrix format only works for 10 or fewer processes" << std::endl;
    }
    return;
  }

  int numRows = matrix.NumGlobalRows();
  int numCols = matrix.NumGlobalCols();

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

  make_my_A(matrix, myA, comm);

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

  delete [] myA;
}

void Test(const string what, Epetra_RowMatrix& A)
{
  T Prec(&A);

  bool UseTranspose = true;

  IFPACK_CHK_ERRV(Prec.Initialize());
  IFPACK_CHK_ERRV(Prec.Compute());
  IFPACK_CHK_ERRV(Prec.SetUseTranspose(UseTranspose));

  Epetra_MultiVector LHS_exact(A.OperatorDomainMap(), 2);
  Epetra_MultiVector LHS(A.OperatorDomainMap(), 2);
  Epetra_MultiVector RHS(A.OperatorRangeMap(), 2);

  LHS_exact.Random(); LHS.PutScalar(0.0);

  A.Multiply(UseTranspose, LHS_exact, RHS);

  Prec.ApplyInverse(RHS, LHS);

  LHS.Update(1.0, LHS_exact, -1.0);
  double norm[2];

  LHS.Norm2(norm);
  norm[0] += norm[1];

  if (norm[0] > 1e-5)
  {
    cout << what << ": Test failed: norm = " << norm[0] << endl;
    exit(EXIT_FAILURE);
  }

  cout << what << ": Test passed: norm = " << norm[0] << endl;
}

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

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

  Epetra_Map map = A.RowMatrixRowMap();
  const Epetra_Comm & comm = map.Comm();
  int numProc = comm.NumProc();

  if (numProc==1 || !A.Map().DistributedGlobal())
    writeRowMatrix(handle, A);
  else {
    int numRows = map.NumMyElements();
    
    Epetra_Map allGidsMap((int_type) -1, numRows, (int_type) 0,comm);
    
    typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
    for (int i=0; i<numRows; i++) allGids[i] = (int_type) map.GID64(i);
    
    // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
    int numChunks = numProc;
    int stripSize = allGids.GlobalLength64()/numChunks;
    int remainder = allGids.GlobalLength64()%numChunks;
    int curStart = 0;
    int curStripSize = 0;
    typename Epetra_GIDTypeSerialDenseVector<int_type>::impl importGidList;
    if (comm.MyPID()==0) 
      importGidList.Size(stripSize+1); // Set size of vector to max needed
    for (int i=0; i<numChunks; i++) {
      if (comm.MyPID()==0) { // Only PE 0 does this part
	curStripSize = stripSize;
	if (i<remainder) curStripSize++; // handle leftovers
	for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
	curStart += curStripSize;
      }
      // The following import map will be non-trivial only on PE 0.
      if (comm.MyPID()>0) assert(curStripSize==0);
      Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
      Epetra_Import gidImporter(importGidMap, allGidsMap);
      typename Epetra_GIDTypeVector<int_type>::impl importGids(importGidMap);
      if (importGids.Import(allGids, gidImporter, Insert)!=0) {EPETRA_CHK_ERR(-1); }

      // importGids now has a list of GIDs for the current strip of matrix rows.
      // Use these values to build another importer that will get rows of the matrix.

      // The following import map will be non-trivial only on PE 0.
      Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), map.IndexBase64(), comm);
      Epetra_Import importer(importMap, map);
      Epetra_CrsMatrix importA(Copy, importMap, 0);
      if (importA.Import(A, importer, Insert)!=0) {EPETRA_CHK_ERR(-1); }
      if (importA.FillComplete(A.OperatorDomainMap(), importMap)!=0) {EPETRA_CHK_ERR(-1);}

      // Finally we are ready to write this strip of the matrix to ostream
      if (writeRowMatrix(handle, importA)!=0) {EPETRA_CHK_ERR(-1);}
    }
  }
  return(0);
}

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

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 power_method(bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0,
		 Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
{
	
  // Fill z with random Numbers
  Epetra_Vector z(z0);
	
  // variable needed for iteration
  double normz, residual;

  int ierr = 1;
	
  for(int iter = 0; iter < niters; iter++) {
		z.Norm2(&normz); // Compute 2-norm of z
		q.Scale(1.0/normz, z);
		A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
		q.Dot(z, lambda); // Approximate maximum eigenvaluE
		if(iter%100==0 || iter+1==niters) {
			resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
			resid.Norm2(&residual);
			if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda
											 << "  Residual of A*q - lambda*q = " << residual << endl;
		}
		if(residual < tolerance) {
			ierr = 0;
			break;
		}
	}
  return(ierr);
}

void ModeLaplace1DQ1::memoryInfo() const {

  int myPid = MyComm.MyPID();

  Epetra_RowMatrix *Mat = dynamic_cast<Epetra_RowMatrix *>(M);
  if ((myPid == 0) && (Mat)) {
    cout << " Total number of nonzero entries in mass matrix      = ";
    cout.width(15);
    cout << Mat->NumGlobalNonzeros() << endl;
    double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
    memSize += 2*Mat->NumGlobalRows()*sizeof(int);
    cout << " Memory requested for mass matrix per processor      = (EST) ";
    cout.precision(2);
    cout.width(6);
    cout.setf(ios::fixed, ios::floatfield);
    cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
    cout << endl;
  }

  Mat = dynamic_cast<Epetra_RowMatrix *>(K);
  if ((myPid == 0) && (Mat)) {
    cout << " Total number of nonzero entries in stiffness matrix = ";
    cout.width(15);
    cout << Mat->NumGlobalNonzeros() << endl;
    double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
    memSize += 2*Mat->NumGlobalRows()*sizeof(int);
    cout << " Memory requested for stiffness matrix per processor = (EST) ";
    cout.precision(2);
    cout.width(6);
    cout.setf(ios::fixed, ios::floatfield);
    cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
    cout << endl;
  }

}

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 EpetraExt::scaleModelFuncFirstDerivOp(
  const Epetra_Vector *invVarScaling,
  const Epetra_Vector *fwdFuncScaling,
  Epetra_Operator *funcDerivOp,
  bool *didScaling
  )
{
  TEUCHOS_TEST_FOR_EXCEPT(0==funcDerivOp);
  TEUCHOS_TEST_FOR_EXCEPT(0==didScaling);
  *didScaling = false; // Assume not scaled to start
  Epetra_RowMatrix *funcDerivRowMatrix
    = dynamic_cast<Epetra_RowMatrix*>(funcDerivOp);
  if (funcDerivRowMatrix) {
    if (fwdFuncScaling)
      funcDerivRowMatrix->LeftScale(*fwdFuncScaling);
    if (invVarScaling)
      funcDerivRowMatrix->RightScale(*invVarScaling);
    *didScaling = true;
    // Note that above I do the func scaling before the var scaling since that
    // is the same order it was done for W in Thyra::EpetraModelEvaluator
  }
}

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

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;

  }
}