C++ MultiVector::putScalar方法代码示例

void ReorderFilter<MatrixType>::
apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
       Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
       Teuchos::ETransp mode,
       scalar_type alpha,
       scalar_type beta) const
{
  typedef Teuchos::ScalarTraits<scalar_type> STS;

  // Note: This isn't AztecOO compliant.  But neither was Ifpack's version.
  // Note: The localized maps mean the matvec is trivial (and has no import)
  TEUCHOS_TEST_FOR_EXCEPTION(
    X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
    "Ifpack2::ReorderFilter::apply: X.getNumVectors() != Y.getNumVectors().");

  const scalar_type zero = STS::zero ();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const scalar_type> > x_ptr = X.get2dView();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<scalar_type> > y_ptr = Y.get2dViewNonConst();

  Y.putScalar (zero);
  const size_t NumVectors = Y.getNumVectors ();

  for (size_t i = 0; i < A_->getNodeNumRows (); ++i) {
    size_t Nnz;
    // Use this class's getrow to make the below code simpler
    getLocalRowCopy (i, Indices_ (), Values_ (), Nnz);
    if (mode == Teuchos::NO_TRANS) {
      for (size_t j = 0; j < Nnz; ++j) {
        for (size_t k = 0; k < NumVectors; ++k) {
          y_ptr[k][i] += Values_[j] * x_ptr[k][Indices_[j]];
        }
      }
    }
    else if (mode == Teuchos::TRANS) {
      for (size_t j = 0; j < Nnz; ++j) {
        for (size_t k = 0; k < NumVectors; ++k) {
          y_ptr[k][Indices_[j]] += Values_[j] * x_ptr[k][i];
        }
      }
    }
    else { //mode==Teuchos::CONJ_TRANS
      for (size_t j = 0; j < Nnz; ++j) {
        for (size_t k = 0; k < NumVectors; ++k) {
          y_ptr[k][Indices_[j]] += STS::conjugate(Values_[j]) * x_ptr[k][i];
        }
      }
    }
  }
}

    /*! \brief This routine takes the Tpetra::MultiVector \c x and applies the operator
      to it resulting in the Tpetra::MultiVector \c y, which is returned.
      \note It is expected that any problem with applying this operator to \c x will be
      indicated by an std::exception being thrown.
    */
    void Apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x, Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& y, ETrans trans = NOTRANS ) const {

      TEUCHOS_TEST_FOR_EXCEPTION(trans!=NOTRANS, MueLuOpFailure,
                         "Belos::MueLuOp::Apply, transpose mode != NOTRANS not supported by MueLu preconditionners.");

      //FIXME InitialGuessIsZero currently does nothing in MueLu::Hierarchy.Iterate(), but it matters for AMGX
      y.putScalar(0.0);

#ifdef HAVE_MUELU_AMGX
      if (!AMGX_.is_null())
        AMGX_->apply(x, y);
#endif

      if (!Hierarchy_.is_null()) {
        Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> & temp_x = const_cast<Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> &>(x);

        const Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tX(rcpFromRef(temp_x));
        Xpetra::TpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tY(rcpFromRef(y));
        Hierarchy_->Iterate(tX, tY, 1, true);
      }
    }

void SingletonFilter<MatrixType>::apply(const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X, 
				       Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
				       Teuchos::ETransp mode, 
				       Scalar alpha,
				       Scalar beta) const
{  
  // Note: This isn't AztecOO compliant.  But neither was Ifpack's version.

  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
			     "Ifpack2::SingletonFilter::apply ERROR: X.getNumVectors() != Y.getNumVectors().");

  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> > x_ptr = X.get2dView();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >       y_ptr = Y.get2dViewNonConst();

  Y.putScalar(zero);
  size_t NumVectors = Y.getNumVectors();


  for (size_t i = 0 ; i < NumRows_ ; ++i) {
    size_t Nnz;
    // Use this class's getrow to make the below code simpler
    getLocalRowCopy(i,Indices_(),Values_(),Nnz);
    if (mode==Teuchos::NO_TRANS){
      for (size_t j = 0 ; j < Nnz ; ++j) 
	for (size_t k = 0 ; k < NumVectors ; ++k)
	  y_ptr[k][i] += Values_[j] * x_ptr[k][Indices_[j]];      
    }
    else if (mode==Teuchos::TRANS){
      for (size_t j = 0 ; j < Nnz ; ++j) 
	for (size_t k = 0 ; k < NumVectors ; ++k)
	  y_ptr[k][Indices_[j]] += Values_[j] * x_ptr[k][i];
    }
    else { //mode==Teuchos::CONJ_TRANS
      for (size_t j = 0 ; j < Nnz ; ++j) 
	for (size_t k = 0 ; k < NumVectors ; ++k)
	  y_ptr[k][Indices_[j]] += Teuchos::ScalarTraits<Scalar>::conjugate(Values_[j]) * x_ptr[k][i];
    }
  }
}

void
OverlappingRowMatrix<MatrixType>::
apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
       Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &Y,
       Teuchos::ETransp mode,
       scalar_type alpha,
       scalar_type beta) const
{
  using Teuchos::ArrayRCP;
  using Teuchos::as;
  typedef scalar_type RangeScalar;
  typedef scalar_type DomainScalar;

  typedef Teuchos::ScalarTraits<RangeScalar> STRS;
  TEUCHOS_TEST_FOR_EXCEPTION(
    X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
    "Ifpack2::OverlappingRowMatrix::apply: The input X and the output Y must "
    "have the same number of columns.  X.getNumVectors() = "
    << X.getNumVectors() << " != Y.getNumVectors() = " << Y.getNumVectors()
    << ".");

  // FIXME (mfh 13 July 2013) This would be a good candidate for a
  // Kokkos local parallel operator implementation.  That would
  // obviate the need for getting views of the data and make the code
  // below a lot simpler.

  const RangeScalar zero = STRS::zero ();
  ArrayRCP<ArrayRCP<const DomainScalar> > x_ptr = X.get2dView();
  ArrayRCP<ArrayRCP<RangeScalar> >        y_ptr = Y.get2dViewNonConst();
  Y.putScalar(zero);
  size_t NumVectors = Y.getNumVectors();

  const size_t numMyRowsA = A_->getNodeNumRows ();
  for (size_t i = 0; i < numMyRowsA; ++i) {
    size_t Nnz;
    // Use this class's getrow to make the below code simpler
    A_->getLocalRowCopy (i, Indices_ (),Values_ (), Nnz);
    if (mode == Teuchos::NO_TRANS) {
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][i] += as<RangeScalar> (Values_[j]) *
            as<RangeScalar> (x_ptr[k][Indices_[j]]);
    }
    else if (mode == Teuchos::TRANS){
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][Indices_[j]] += as<RangeScalar> (Values_[j]) *
            as<RangeScalar> (x_ptr[k][i]);
    }
    else { // mode == Teuchos::CONJ_TRANS
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][Indices_[j]] +=
            STRS::conjugate (as<RangeScalar> (Values_[j])) *
            as<RangeScalar> (x_ptr[k][i]);
    }
  }

  const size_t numMyRowsB = ExtMatrix_->getNodeNumRows ();
  for (size_t i = 0 ; i < numMyRowsB ; ++i) {
    size_t Nnz;
    // Use this class's getrow to make the below code simpler
    ExtMatrix_->getLocalRowCopy (i, Indices_ (), Values_ (), Nnz);
    if (mode == Teuchos::NO_TRANS) {
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][numMyRowsA+i] += as<RangeScalar> (Values_[j]) *
            as<RangeScalar> (x_ptr[k][Indices_[j]]);
    }
    else if (mode == Teuchos::TRANS) {
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][numMyRowsA+Indices_[j]] += as<RangeScalar> (Values_[j]) *
            as<RangeScalar> (x_ptr[k][i]);
    }
    else { // mode == Teuchos::CONJ_TRANS
      for (size_t j = 0; j < Nnz; ++j)
        for (size_t k = 0; k < NumVectors; ++k)
          y_ptr[k][numMyRowsA+Indices_[j]] +=
            STRS::conjugate (as<RangeScalar> (Values_[j])) *
            as<RangeScalar> (x_ptr[k][i]);
    }
  }
}

void LocalSparseTriangularSolver<MatrixType>::
apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type,
         global_ordinal_type, node_type>& X,
       Tpetra::MultiVector<scalar_type, local_ordinal_type,
         global_ordinal_type, node_type>& Y,
       Teuchos::ETransp mode,
       scalar_type alpha,
       scalar_type beta) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  typedef scalar_type ST;
  typedef Teuchos::ScalarTraits<ST> STS;
  const char prefix[] = "Ifpack2::LocalSparseTriangularSolver::apply: ";
  if (! out_.is_null ()) {
    *out_ << ">>> DEBUG " << prefix;
    if (A_crs_.is_null ()) {
      *out_ << "A_crs_ is null!" << std::endl;
    }
    else {
      const std::string uplo = A_crs_->isUpperTriangular () ? "U" :
        (A_crs_->isLowerTriangular () ? "L" : "N");
      const std::string trans = (mode == Teuchos::CONJ_TRANS) ? "C" :
        (mode == Teuchos::TRANS ? "T" : "N");
      const std::string diag =
        (A_crs_->getNodeNumDiags () < A_crs_->getNodeNumRows ()) ? "U" : "N";
      *out_ << "uplo=\"" << uplo
            << "\", trans=\"" << trans
            << "\", diag=\"" << diag << "\"" << std::endl;
    }
  }

  TEUCHOS_TEST_FOR_EXCEPTION
    (! isComputed (), std::runtime_error, prefix << "If compute() has not yet "
     "been called, or if you have changed the matrix via setMatrix(), you must "
     "call compute() before you may call this method.");
  // If isComputed() is true, it's impossible for the matrix to be
  // null, or for it not to be a Tpetra::CrsMatrix.
  TEUCHOS_TEST_FOR_EXCEPTION
    (A_.is_null (), std::logic_error, prefix << "A_ is null.  "
     "Please report this bug to the Ifpack2 developers.");
  TEUCHOS_TEST_FOR_EXCEPTION
    (A_crs_.is_null (), std::logic_error, prefix << "A_crs_ is null.  "
     "Please report this bug to the Ifpack2 developers.");
  // However, it _is_ possible that the user called resumeFill() on
  // the matrix, after calling compute().  This is NOT allowed.
  TEUCHOS_TEST_FOR_EXCEPTION
    (! A_crs_->isFillComplete (), std::runtime_error, "If you call this "
     "method, the matrix must be fill complete.  It is not.  This means that "
     " you must have called resumeFill() on the matrix before calling apply(). "
     "This is NOT allowed.  Note that this class may use the matrix's data in "
     "place without copying it.  Thus, you cannot change the matrix and expect "
     "the solver to stay the same.  If you have changed the matrix, first call "
     "fillComplete() on it, then call compute() on this object, before you call"
     " apply().  You do NOT need to call setMatrix, as long as the matrix "
     "itself (that is, its address in memory) is the same.");

  auto G = A_->getGraph ();
  TEUCHOS_TEST_FOR_EXCEPTION
    (G.is_null (), std::logic_error, prefix << "A_ and A_crs_ are nonnull, "
     "but A_'s RowGraph G is null.  "
     "Please report this bug to the Ifpack2 developers.");
  auto importer = G->getImporter ();
  auto exporter = G->getExporter ();

  if (! importer.is_null ()) {
    if (X_colMap_.is_null () || X_colMap_->getNumVectors () != X.getNumVectors ()) {
      X_colMap_ = rcp (new MV (importer->getTargetMap (), X.getNumVectors ()));
    }
    else {
      X_colMap_->putScalar (STS::zero ());
    }
    // See discussion of Github Issue #672 for why the Import needs to
    // use the ZERO CombineMode.  The case where the Export is
    // nontrivial is likely never exercised.
    X_colMap_->doImport (X, *importer, Tpetra::ZERO);
  }
  RCP<const MV> X_cur = importer.is_null () ? rcpFromRef (X) :
    Teuchos::rcp_const_cast<const MV> (X_colMap_);

  if (! exporter.is_null ()) {
    if (Y_rowMap_.is_null () || Y_rowMap_->getNumVectors () != Y.getNumVectors ()) {
      Y_rowMap_ = rcp (new MV (exporter->getSourceMap (), Y.getNumVectors ()));
    }
    else {
      Y_rowMap_->putScalar (STS::zero ());
    }
    Y_rowMap_->doExport (Y, *importer, Tpetra::ADD);
  }
  RCP<MV> Y_cur = exporter.is_null () ? rcpFromRef (Y) : Y_rowMap_;

  localApply (*X_cur, *Y_cur, mode, alpha, beta);

  if (! exporter.is_null ()) {
    Y.putScalar (STS::zero ());
    Y.doExport (*Y_cur, *exporter, Tpetra::ADD);
  }

  ++numApply_;
//.........这里部分代码省略.........