C++ Epetra_CrsMatrix类代码示例

bool probing_test(Epetra_CrsMatrix & in_mat, bool build_list){
  const Epetra_CrsGraph & graph=in_mat.Graph();

  Teuchos::ParameterList main, zoltan;

  if(build_list){
//     zoltan.set("DISTANCE","2");
    main.set("ZOLTAN",zoltan);
  }

  Isorropia::Epetra::Prober prober(Teuchos::rcp<const Epetra_CrsGraph>(&graph,false),main);
  Epetra_CrsMatrix out_mat(Copy,graph);
  int rv=prober.probe(in_mat,out_mat);
  if(rv!=0) {printf("ERROR: probing failed\n");return false;}
  
#ifdef HAVE_EPETRAEXT
  EpetraExt::MatrixMatrix::Add(in_mat,false,1,out_mat,-1);
  double nrm=out_mat.NormInf()/in_mat.NormInf();
  if(!in_mat.Comm().MyPID())
    printf("diff norm = %22.16e\n",nrm);
  if(nrm < 1e-12) return true;
  else return false;
#endif
  return true;
}

/*!
 * \brief Solve.
 */
void JacobiSolver::iterate( const int max_iters, const double tolerance )
{
    // Extract the linear problem.
    Epetra_CrsMatrix *A = 
	dynamic_cast<Epetra_CrsMatrix*>( d_linear_problem->GetMatrix() );
    Epetra_Vector *x = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetLHS() );
    const Epetra_Vector *b = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetRHS() );

    // Setup the residual.
    Epetra_Map row_map = A->RowMap();
    Epetra_Vector residual( row_map );

    // Iterate.
    Epetra_CrsMatrix H = buildH( A );
    Epetra_Vector temp_vec( row_map );
    d_num_iters = 0;
    double residual_norm = 1.0;
    double b_norm;
    b->NormInf( &b_norm );
    double conv_crit = b_norm*tolerance;
    while ( residual_norm > conv_crit && d_num_iters < max_iters )
    {
	H.Apply( *x, temp_vec );
	x->Update( 1.0, temp_vec, 1.0, *b, 0.0 );

	A->Apply( *x, temp_vec );
	residual.Update( 1.0, *b, -1.0, temp_vec, 0.0 );

	residual.NormInf( &residual_norm );
	++d_num_iters;
    }
}

//==========================================================================
int Ifpack_CrsRiluk::InitValues(const Epetra_CrsMatrix & A) {

  UserMatrixIsCrs_ = true;

  if (!Allocated()) AllocateCrs();

  Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (Epetra_CrsMatrix *) &A, false );

  if (IsOverlapped_) {
  
    OverlapA = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph()) );
    EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));
    EPETRA_CHK_ERR(OverlapA->FillComplete());
  }
  
  // Get Maximun Row length
  int MaxNumEntries = OverlapA->MaxNumEntries();

  // Set L range map and U domain map
  U_DomainMap_ = Teuchos::rcp( &(A.DomainMap()), false );
  L_RangeMap_ = Teuchos::rcp( &(A.RangeMap()), false );
  // Do the rest using generic Epetra_RowMatrix interface

  EPETRA_CHK_ERR(InitAllValues(*OverlapA, MaxNumEntries));

  return(0);
}

// Simple Power method algorithm
double power_method(const Epetra_CrsMatrix& A) {  
  // variable needed for iteration
  double lambda = 0.0;
  int niters = A.RowMap().NumGlobalElements()*10;
  double tolerance = 1.0e-10;
  // Create vectors
  Epetra_Vector q(A.RowMap());
  Epetra_Vector z(A.RowMap());
  Epetra_Vector resid(A.RowMap());
  // Fill z with random Numbers
  z.Random();
  // variable needed for iteration
  double normz;
  double residual = 0;
  int iter = 0;
  while (iter==0 || (iter < niters && residual > tolerance)) {
    z.Norm2(&normz); // Compute 2-norm of z
    q.Scale(1.0/normz, z);
    A.Multiply(false, q, z); // Compute z = A*q
    q.Dot(z, &lambda); // Approximate maximum eigenvalue
    if (iter%10==0 || iter+1==niters) {
      // Compute A*q - lambda*q every 10 iterations
      resid.Update(1.0, z, -lambda, q, 0.0);
      resid.Norm2(&residual);
      if (q.Map().Comm().MyPID()==0)
	std::cout << "Iter = " << iter << "  Lambda = " << lambda 
	     << "  Two-norm of A*q - lambda*q = " 
	     << residual << std::endl;
    } 
    iter++;
  }
  return(lambda);
}

int test_azoo_conv_anorm(Epetra_CrsMatrix& A,
                         Epetra_Vector& x,
                         Epetra_Vector& b,
                         bool verbose)
{
  if (verbose) {
    cout << "testing AztecOO with AZ_conv = AZ_Anorm" << endl;
  }

  Epetra_Vector soln_Anorm(x), soln_none(x), vec1(x), rhs(x);

  //We'll put large numbers in a vector and use that to generate an rhs
  //which has norm much larger than the infinity-norm of the matrix.
  vec1.PutScalar(1000.0);
  soln_Anorm.PutScalar(0.0);
  soln_none.PutScalar(0.0);

  A.Multiply(false, vec1, rhs);

  AztecOO azoo(&A, &soln_Anorm, &rhs);
  azoo.SetAztecOption(AZ_conv, AZ_Anorm);
  azoo.SetAztecOption(AZ_solver, AZ_cg);
  //azoo.SetAztecOption(AZ_scaling, AZ_sym_diag);

  azoo.Iterate(30, 1.e-5);

  AztecOO azoo1(&A, &soln_none, &rhs);
  azoo1.SetAztecOption(AZ_conv, AZ_rhs);
  azoo1.SetAztecOption(AZ_solver, AZ_cg);

  azoo1.Iterate(30, 1.e-5);

  double rhsnorm = 0.0;
  rhs.Norm2(&rhsnorm);

  double anorm = A.NormInf();

  double rnrm_anorm = resid2norm(A, soln_Anorm, rhs);
  double rnrm_rhs = resid2norm(A, soln_none, rhs);

  //we expect the ratio rnrm_anorm/rnrm_rhs to roughly equal
  //the ratio anorm/rhsnorm, since rnrm_anorm is the residual norm
  //obtained for the solve that used Anorm scaling to determine
  //convergence, and rnrm_rhs is the residual norm obtained for
  //the solve that used rhs scaling.

  double ratio1 = rnrm_anorm/rnrm_rhs;
  double ratio2 = anorm/rhsnorm;

  cout << "ratio1: " << ratio1 << ", ratio2: " << ratio2 << endl;
  if (std::abs(ratio1 - ratio2) > 1.e-1) {
    if (verbose) {
      cout << "anorm: " << anorm << ", rhsnorm: " << rhsnorm
       << "rnrm_anorm: " << rnrm_anorm << ", rnrm_rhs: " << rnrm_rhs<<endl;
    }
    return(-1);
  }

  return(0);
}

int InitMatValues( const Epetra_CrsMatrix& newA, Epetra_CrsMatrix* A )
{
  int numMyRows = newA.NumMyRows();
  int maxNum = newA.MaxNumEntries();
  int numIn;
  int *idx = 0;
  double *vals = 0;

  idx = new int[maxNum];
  vals = new double[maxNum];

  // For each row get the values and indices, and replace the values in A.
  for (int i=0; i<numMyRows; ++i) {

    // Get the values and indices from newA.
    EPETRA_CHK_ERR( newA.ExtractMyRowCopy(i, maxNum, numIn, vals, idx) );

    // Replace the values in A
    EPETRA_CHK_ERR( A->ReplaceMyValues(i, numIn, vals, idx) );

  }

  // Clean up.
  delete [] idx;
  delete [] vals;

  return 0;
}

Epetra_CrsMatrix* create_and_fill_crs_matrix(const Epetra_Map& emap)
{
  int localproc = emap.Comm().MyPID();
  int local_n = emap.NumMyElements();
  int global_n = emap.NumGlobalElements();
  int myFirstGlobalRow = localproc*local_n;
  int globalCols[3];
  double values[3];
  Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, emap, 3);

  for(int i=0; i<local_n; ++i) {
    int globalRow = myFirstGlobalRow +i;

    int numcols = 0;
    if (globalRow > 0) {
      globalCols[numcols] = globalRow-1;
      values[numcols++] = -1.0;
    }

    globalCols[numcols] = globalRow;
    values[numcols++] = 4.0;

    if (globalRow < global_n-1) {
      globalCols[numcols] = globalRow+1;
      values[numcols++] = -1.0;
    }

    A->InsertGlobalValues(globalRow, numcols, values, globalCols);
  }

  A->FillComplete();

  return A;
}

/*----------------------------------------------------------------------*
 |                                                           m.gee 11/05|
 *----------------------------------------------------------------------*/
Epetra_CrsMatrix* ML_NOX::StripZeros(Epetra_CrsMatrix& in, double eps)
{
  Epetra_CrsMatrix* out = new Epetra_CrsMatrix(Copy,in.RowMap(),200);
  for (int lrow=0; lrow<in.NumMyRows(); ++lrow)
  {
    int grow = in.GRID(lrow); 
    if (grow<0) { cout << "ERROR: grow<0 \n"; exit(0); }
    int numentries;
    int* lindices;
    double* values;
    int err  = in.ExtractMyRowView(lrow,numentries,values,lindices);
    if (err) { cout << "ExtractMyRowView returned " << err << endl; exit(0);}
    for (int j=0; j<numentries; ++j)
    {
      int lcol = lindices[j];  
      int gcol = in.GCID(lcol); 
      if (gcol<0) { cout << "ERROR: gcol<0 \n"; exit(0); }
      if (abs(values[j])<eps && gcol != grow)
        continue;
      int err = out->InsertGlobalValues(grow,1,&values[j],&gcol);
      if (err) { cout << "InsertGlobalValues returned " << err << endl; exit(0);}
    }
  }
  out->FillComplete();
  out->OptimizeStorage();
  return out;
}

/* ml_epetra_data_pack::setup - This function does the setup phase for ML_Epetra, pulling
   key parameters from the Teuchos list, and calling the aggregation routines
   Parameters:
   N       - Number of unknowns [I]
   rowind  - Row indices of matrix (CSC format) [I]
   colptr  - Column indices of matrix (CSC format) [I]
   vals    - Nonzero values of matrix (CSC format) [I]
   Returns: IS_TRUE if setup was succesful, IS_FALSE otherwise
*/
int ml_epetra_data_pack::setup(int N,int* rowind,int* colptr, double* vals){
  int i,j;
  int *rnz;
  
  /* Nonzero counts for Epetra */
  rnz=new int[N];
  for(i=0;i<N;i++) rnz[i]=rowind[i+1]-rowind[i];  
  
  /* Epetra Setup */
  Comm= new Epetra_SerialComm;
  Map=new Epetra_Map(N,0,*Comm);
  A=new Epetra_CrsMatrix(Copy,*Map,rnz);
  
  /* Do the matrix assembly */
  for(i=0;i<N;i++)
    for(j=colptr[i];j<colptr[i+1];j++)
      A->InsertGlobalValues(rowind[j],1,&vals[j],&i);
  //NTS: Redo with block assembly, remembering to transpose
  A->FillComplete();

  /* Allocate Memory */
  Prec=new MultiLevelPreconditioner(*A, *List,false);  
  
  /* Build the Heirarchy */
  Prec->ComputePreconditioner();

  /* Pull Operator Complexity */
  operator_complexity = Prec->GetML_Aggregate()->operator_complexity / Prec->GetML_Aggregate()->fine_complexity;

  /* Cleanup */
  delete rnz;
  
  return IS_TRUE;
}/*end setup*/

// B here is the "reduced" matrix.  Square matrices w/ Row=Domain=Range only.
double test_with_matvec_reduced_maps(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, const Epetra_Map & Bfullmap){
  const Epetra_Map & Amap  = A.DomainMap();
  Epetra_Vector Xa(Amap), Ya(Amap), Diff(Amap);
  const Epetra_Map *Bmap  = Bfullmap.NumMyElements() > 0 ? &B.DomainMap() : 0;
  Epetra_Vector *Xb = Bmap ? new Epetra_Vector(*Bmap) : 0;
  Epetra_Vector *Yb = Bmap ? new Epetra_Vector(*Bmap) : 0;

  Epetra_Vector Xb_alias(View,Bfullmap, Bmap ? Xb->Values(): 0);
  Epetra_Vector Yb_alias(View,Bfullmap, Bmap ? Yb->Values(): 0);

  Epetra_Import Ximport(Bfullmap,Amap);

  // Set the input vector
  Xa.SetSeed(24601);
  Xa.Random();
  Xb_alias.Import(Xa,Ximport,Insert);

  // Do the multiplies
  A.Apply(Xa,Ya);
  if(Bmap) B.Apply(*Xb,*Yb);

  // Check solution
  Epetra_Import Yimport(Amap,Bfullmap);
  Diff.Import(Yb_alias,Yimport,Insert);


  Diff.Update(-1.0,Ya,1.0);
  double norm;
  Diff.Norm2(&norm);

  delete Xb; delete Yb;
  return norm;
}

//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_;

}

int build_matrix_unfused(const Epetra_CrsMatrix & SourceMatrix, Epetra_Export & RowExporter, Epetra_CrsMatrix *&A){
  int rv=0;
  rv=A->Export(SourceMatrix, RowExporter, Insert);
  if(rv) {cerr<<"build_matrix_unfused: Export failed"<<endl; return rv;}

  rv=A->FillComplete(SourceMatrix.DomainMap(), SourceMatrix.RangeMap());
  return rv;
}

int runHypreTest(Epetra_CrsMatrix &A){
  
  Epetra_Vector X(A.RowMap());
  EPETRA_CHK_ERR(X.Random());
  Epetra_Vector Y(A.RowMap());
  EPETRA_CHK_ERR(A.Multiply(false, X, Y));
  
  return 0;
}

shared_ptr<Epetra_CrsMatrix> sparseCholesky(const Epetra_CrsMatrix &mat) {
  // Note: we assume the matrix mat is symmetric and positive-definite
  size_t size = mat.NumGlobalCols();
  if (mat.NumGlobalRows() != size)
    throw std::invalid_argument("sparseCholesky(): matrix must be square");

  int *rowOffsets = 0;
  int *colIndices = 0;
  double *values = 0;
  mat.ExtractCrsDataPointers(rowOffsets, colIndices, values);

  Epetra_SerialComm comm;
  Epetra_LocalMap rowMap(static_cast<int>(size), 0 /* index_base */, comm);
  Epetra_LocalMap columnMap(static_cast<int>(size), 0 /* index_base */, comm);
  shared_ptr<Epetra_CrsMatrix> result = boost::make_shared<Epetra_CrsMatrix>(
      Copy, rowMap, columnMap, mat.GlobalMaxNumEntries());

  arma::Mat<double> localMat;
  arma::Mat<double> localCholesky;
  std::vector<bool> processed(size, false);
  for (size_t r = 0; r < size; ++r) {
    if (processed[r])
      continue;
    int localSize = rowOffsets[r + 1] - rowOffsets[r];
    localMat.set_size(localSize, localSize);
    localMat.fill(0.);
    localCholesky.set_size(localSize, localSize);
    for (int s = 0; s < localSize; ++s) {
      int row = colIndices[rowOffsets[r] + s];
      for (int c = 0; c < localSize; ++c) {
        int col = colIndices[rowOffsets[row] + c];
        if (col != colIndices[rowOffsets[r] + c])
          throw std::invalid_argument("sparseCholesky(): matrix is not "
                                      "block-diagonal");
        localMat(s, c) = values[rowOffsets[row] + c];
      }
    }
    assert(arma::norm(localMat - localMat.t(), "fro") <
           1e-12 * arma::norm(localMat, "fro"));
    localCholesky = arma::chol(localMat); // localCholesky: U
    for (int s = 0; s < localSize; ++s) {
      int row = colIndices[rowOffsets[r] + s];
      processed[row] = true;
#ifndef NDEBUG
      int errorCode =
#endif
          result->InsertGlobalValues(row, s + 1 /* number of values */,
                                     localCholesky.colptr(s),
                                     colIndices + rowOffsets[r]);
      assert(errorCode == 0);
    }
  }
  result->FillComplete(columnMap, rowMap);

  return result;
}

//==============================================================================
Epetra_FastCrsMatrix::Epetra_FastCrsMatrix(const Epetra_CrsMatrix & Matrix, bool UseFloats) 
  : CrsMatrix_(Matrix),
    Values_(0),
    NumMyRows_(Matrix.NumMyRows()),
    NumMyNonzeros(Matrix.NumMyNonzeros()),
    ImportVector_(0),
    ExportVector_(0),
    CV_(Copy)
{
  if (!CrsMatrix_.Filled()) throw CrsMatrix_.ReportError("Input matrix must have called FillComplete()", -1);
  Allocate(UseFloats);
}