C++ Epetra_Vector::Values方法代码示例

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

void EpetraLinearOp::computeAbsRowSum(Epetra_Vector & x) const
{
  TEUCHOS_ASSERT(!is_null(rowMatrix_));

  RCP<Epetra_CrsMatrix> crsMatrix 
    = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(rowMatrix_);

  TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
    Exceptions::OpNotSupported,
    "EpetraLinearOp::computeAbsRowSum(...): wrapped matrix must be of type "
    "Epetra_CrsMatrix for this method. Other operator types are not supported."
    );

  //
  // Put inverse of the sum of absolute values of the ith row of A in x[i].
  // (this is a modified copy of Epetra_CrsMatrix::InvRowSums)
  //

  if (crsMatrix->Filled()) {
    TEUCHOS_TEST_FOR_EXCEPTION(is_null(crsMatrix),
      std::invalid_argument,
      "EpetraLinearOp::computeAbsRowSum(...): Epetra_CrsMatrix must be filled"
    );
  } 
  int i, j;
  x.PutScalar(0.0); // Make sure we sum into a vector of zeros.
  double * xp = (double*)x.Values();
  if (crsMatrix->Graph().RangeMap().SameAs(x.Map()) && crsMatrix->Exporter() != 0) {
    Epetra_Vector x_tmp(crsMatrix->RowMap());
    x_tmp.PutScalar(0.0);
    double * x_tmp_p = (double*)x_tmp.Values();
    for (i=0; i < crsMatrix->NumMyRows(); i++) {
      int      NumEntries = 0;
      double * RowValues  = 0;
      crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
      for (j=0; j < NumEntries; j++)  x_tmp_p[i] += std::abs(RowValues[j]);
    }
    TEUCHOS_TEST_FOR_EXCEPT(0!=x.Export(x_tmp, *crsMatrix->Exporter(), Add)); //Export partial row sums to x.
  }
  else if (crsMatrix->Graph().RowMap().SameAs(x.Map())) {
    for (i=0; i < crsMatrix->NumMyRows(); i++) {
      int      NumEntries = 0;
      double * RowValues  = 0;
      crsMatrix->ExtractMyRowView(i,NumEntries,RowValues);
      double scale = 0.0;
      for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
      xp[i] = scale;
    }
  }
  else { // x.Map different than both crsMatrix->Graph().RowMap() and crsMatrix->Graph().RangeMap()
    TEUCHOS_TEST_FOR_EXCEPT(true); // The map of x must be the RowMap or RangeMap of A.
  }
}

//=============================================================================
int Epetra_FastCrsMatrix::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const {
  //
  // This function find y such that Ly = x or Uy = x or the transpose cases.
  //

  if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.

  if ((Upper) && (!UpperTriangular())) EPETRA_CHK_ERR(-2);
  if ((!Upper) && (!LowerTriangular())) EPETRA_CHK_ERR(-3);
  if ((!UnitDiagonal) && (NoDiagonal())) EPETRA_CHK_ERR(-4); // If matrix has no diagonal, we must use UnitDiagonal
  if ((!UnitDiagonal) && (NumMyDiagonals()<NumMyRows_)) EPETRA_CHK_ERR(-5); // Need each row to have a diagonal
      

  int i, j, j0;
  int * NumEntriesPerRow = NumEntriesPerRow_;
  int ** Indices = Indices_;
  double ** Values = Values_;
  int NumMyCols_ = NumMyCols();

  // If upper, point to last row
  if ((Upper && !Trans) || (!Upper && Trans)) {
    NumEntriesPerRow += NumMyRows_-1;
    Indices += NumMyRows_-1;
    Values += NumMyRows_-1;
  }
    
  double *xp = (double*)x.Values();
  double *yp = (double*)y.Values();

  if (!Trans) {

    if (Upper) {

      j0 = 1;
      if (NoDiagonal()) j0--; // Include first term if no diagonal
      for (i=NumMyRows_-1; i >=0; i--) {
	int      NumEntries = *NumEntriesPerRow--;
	int *    RowIndices = *Indices--;
	double * RowValues  = *Values--;
	double sum = 0.0;
	for (j=j0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
	
	if (UnitDiagonal) yp[i] = xp[i] - sum;
	else yp[i] = (xp[i] - sum)/RowValues[0];

      }
    }
    else {
      j0 = 1;
      if (NoDiagonal()) j0--; // Include first term if no diagonal
      for (i=0; i < NumMyRows_; i++) {
	int      NumEntries = *NumEntriesPerRow++ - j0;
	int *    RowIndices = *Indices++;
	double * RowValues  = *Values++;
	double sum = 0.0;
	for (j=0; j < NumEntries; j++) sum += RowValues[j] * yp[RowIndices[j]];
	
	if (UnitDiagonal) yp[i] = xp[i] - sum;
	else yp[i] = (xp[i] - sum)/RowValues[NumEntries];

      }
    }
  }

  // ***********  Transpose case *******************************

  else {

    if (xp!=yp) for (i=0; i < NumMyCols_; i++) yp[i] = xp[i]; // Initialize y for transpose solve
    
    if (Upper) {

      j0 = 1;
      if (NoDiagonal()) j0--; // Include first term if no diagonal
    
      for (i=0; i < NumMyRows_; i++) {
	int      NumEntries = *NumEntriesPerRow++;
	int *    RowIndices = *Indices++;
	double * RowValues  = *Values++;
	if (!UnitDiagonal) yp[i] = yp[i]/RowValues[0];
	for (j=j0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
      }
    }
    else {

      j0 = 1;
      if (NoDiagonal()) j0--; // Include first term if no diagonal
    
      for (i=NumMyRows_-1; i >= 0; i--) {
	int      NumEntries = *NumEntriesPerRow-- - j0;
	int *    RowIndices = *Indices--;
	double * RowValues  = *Values--;
	if (!UnitDiagonal)  yp[i] = yp[i]/RowValues[NumEntries];
	for (j=0; j < NumEntries; j++) yp[RowIndices[j]] -= RowValues[j] * yp[i];
      }
    }

  }
  UpdateFlops(2*NumGlobalNonzeros64());
//.........这里部分代码省略.........

//=============================================================================
int Epetra_FastCrsMatrix::Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const {
  //
  // This function forms the product y = A * x or y = A' * x
  //

  int i, j;
  double * xp = (double*)x.Values();
  double *yp = (double*)y.Values();
  int NumMyCols_ = NumMyCols();


  if (!TransA) {

    // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
    if (Importer()!=0) {
      if (ImportVector_!=0) {
	if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
      }
      if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
      ImportVector_->Import(x, *Importer(), Insert);
      xp = (double*)ImportVector_->Values();
    }

    // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
    if (Exporter()!=0) {
      if (ExportVector_!=0) {
	if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
      }
      if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
      yp = (double*)ExportVector_->Values();
    }

    // Do actual computation

    for (i=0; i < NumMyRows_; i++) {
      int      NumEntries = *NumEntriesPerRow++;
      int *    RowIndices = *Indices++;
      double * RowValues  = *Values++;
      double sum = 0.0;
      for (j=0; j < NumEntries; j++) sum += RowValues[j] * xp[RowIndices[j]];

      yp[i] = sum;

    }
    if (Exporter()!=0) y.Export(*ExportVector_, *Exporter(), Add); // Fill y with Values from export vector
  }

  else { // Transpose operation


    // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors

    if (Exporter()!=0) {
      if (ExportVector_!=0) {
	if (ExportVector_->NumVectors()!=1) { delete ExportVector_; ExportVector_= 0;}
      }
      if (ExportVector_==0) ExportVector_ = new Epetra_MultiVector(RowMap(),1); // Create Export vector if needed
      ExportVector_->Import(x, *Exporter(), Insert);
      xp = (double*)ExportVector_->Values();
    }

    // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
    if (Importer()!=0) {
      if (ImportVector_!=0) {
	if (ImportVector_->NumVectors()!=1) { delete ImportVector_; ImportVector_= 0;}
      }
      if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(ColMap(),1); // Create import vector if needed
      yp = (double*)ImportVector_->Values();
    }

    // Do actual computation

    for (i=0; i < NumMyCols_; i++) yp[i] = 0.0; // Initialize y for transpose multiply
        
    for (i=0; i < NumMyRows_; i++) {
      int      NumEntries = *NumEntriesPerRow++;
      int *    RowIndices = *Indices++;
      double * RowValues  = *Values++;
      for (j=0; j < NumEntries; j++) yp[RowIndices[j]] += RowValues[j] * xp[i];
    }
    if (Importer()!=0) y.Export(*ImportVector_, *Importer(), Add); // Fill y with Values from export vector
  }

  UpdateFlops(2*NumGlobalNonzeros64());
  return(0);
}