C++ tpetra::MultiVector类代码示例

// ============================================================================
void
BorderingHelpers::
dissect(const Tpetra::MultiVector<double,int,int> & x,
        Tpetra::MultiVector<double,int,int> & xSmall,
        double * lambda
       )
{
#ifndef NDEBUG
  TEUCHOS_ASSERT_EQUALITY(x.NumVectors(), xSmall.NumVectors());
  // Make sure the maps are matching.
  std::shared_ptr<const Tpetra::Map<int,int>> extendedMap =
    nosh::BorderingHelpers::extendMapBy1(xSmall.getMap());
  TEUCHOS_ASSERT(x.getMap().SameAs(*extendedMap));
#endif

  Epetra_Import importer(xSmall.getMap(), x.getMap());

  // Strip off the phase constraint variable.
  xSmall.Import(x, importer, Insert);

  // TODO Check if we need lambda on all procs.
  if (x.getMap().Comm().MyPID() == 0) {
    const int n = x.MyLength();
    for (int k = 0; k < x.NumVectors(); k++)
      lambda[k] = (*(x(k)))[n - 1];
  }

  return;
}

void SingletonFilter<MatrixType>::CreateReducedRHSTempl(const Tpetra::MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& LHS,
                                                        const Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& RHS,
                                                        Tpetra::MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& ReducedRHS)
{
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const RangeScalar > > RHS_ptr = RHS.get2dView();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<const DomainScalar> > LHS_ptr = LHS.get2dView();
  Teuchos::ArrayRCP<Teuchos::ArrayRCP<RangeScalar> >        ReducedRHS_ptr = ReducedRHS.get2dViewNonConst();

  size_t NumVectors = LHS.getNumVectors();

  for (size_t i = 0 ; i < NumRows_ ; ++i)
    for (size_t k = 0 ; k < NumVectors ; ++k)
      ReducedRHS_ptr[k][i] = RHS_ptr[k][InvReorder_[i]];

  for (size_t i = 0 ; i < NumRows_ ; ++i) {
    LocalOrdinal ii = InvReorder_[i];
    size_t Nnz;
    A_->getLocalRowCopy(ii,Indices_(),Values_(),Nnz);

    for (size_t j = 0 ; j < Nnz ; ++j) {
      if (Reorder_[Indices_[j]] == -1) {
        for (size_t k = 0 ; k < NumVectors ; ++k)
          ReducedRHS_ptr[k][i] -= (RangeScalar)Values_[j] * (RangeScalar)LHS_ptr[k][Indices_[j]];
      }
    }
  }
}

void TpetraOperatorWrapper::apply(const Tpetra::MultiVector<ST,LO,GO,NT>& X, Tpetra::MultiVector<ST,LO,GO,NT>& Y,Teuchos::ETransp mode,ST alpha, ST beta) const
{
   if (!useTranspose_)
   {
       // allocate space for each vector
       RCP<Thyra::MultiVectorBase<ST> > tX;
       RCP<Thyra::MultiVectorBase<ST> > tY; 

       tX = Thyra::createMembers(thyraOp_->domain(),X.getNumVectors()); 
       tY = Thyra::createMembers(thyraOp_->range(),X.getNumVectors());

       Thyra::assign(tX.ptr(),0.0);
       Thyra::assign(tY.ptr(),0.0);

       // copy epetra X into thyra X
       mapStrategy_->copyTpetraIntoThyra(X, tX.ptr());
       mapStrategy_->copyTpetraIntoThyra(Y, tY.ptr()); // if this matrix isn't block square, this probably won't work!

       // perform matrix vector multiplication
       thyraOp_->apply(Thyra::NOTRANS,*tX,tY.ptr(),alpha,beta);

       // copy thyra Y into epetra Y
       mapStrategy_->copyThyraIntoTpetra(tY, Y);
   }
   else
   {
       TEUCHOS_ASSERT(false);
   }
}

void 
Chebyshev<MatrixType>::
applyMat (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) const
{
  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
   "Ifpack2::Chebyshev::applyMat(): X.getNumVectors() != Y.getNumVectors().");
  impl_.getMatrix ()->apply (X, Y, mode);
}

void
Redistributor<Node>::redistribute_reverse(const  ::Tpetra::MultiVector<double,int,int,Node> & input_vector,  ::Tpetra::MultiVector<double,int,int,Node> & output_vector)
{
  if (!created_importer_) {
    create_importer(input_vector.Map());
  }

  // Export using the importer
  output_vector.Export(input_vector, *importer_, ::Tpetra::INSERT);

}

void BorderedOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node >::apply(
     const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node >& X,
           Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node >& Y,
     Teuchos::ETransp mode, 
     Scalar coefAx, 
     Scalar coefY ) const 
{
  //bool opHasTrans = A_->hasTransposeApply();
  //TEUCHOS_TEST_FOR_EXCEPTION( mode  &&  !opHasTrans, std::runtime_error,
  //"Ifpack2::BorderedOperator::apply() ERROR: The operator does not implement transpose.");
  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
     "Ifpack2::BorderedOperator::apply() ERROR: X.getNumVectors() != Y.getNumVectors().");
  A_->apply(X, Y, mode, coefAx, coefY );
}

void Hiptmair<MatrixType>::
applyHiptmairSmoother(const Tpetra::MultiVector<typename MatrixType::scalar_type,
                      typename MatrixType::local_ordinal_type,
                      typename MatrixType::global_ordinal_type,
                      typename MatrixType::node_type>& X,
                      Tpetra::MultiVector<typename MatrixType::scalar_type,
                      typename MatrixType::local_ordinal_type,
                      typename MatrixType::global_ordinal_type,
                      typename MatrixType::node_type>& Y) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
    global_ordinal_type, node_type> MV;
  const scalar_type ZERO = STS::zero ();
  const scalar_type ONE = STS::one ();

  RCP<MV> res1 = rcp (new MV (A_->getRowMap (), X.getNumVectors ()));
  RCP<MV> vec1 = rcp (new MV (A_->getRowMap (), X.getNumVectors ()));
  RCP<MV> res2 = rcp (new MV (PtAP_->getRowMap (), X.getNumVectors ()));
  RCP<MV> vec2 = rcp (new MV (PtAP_->getRowMap (), X.getNumVectors ()));

  if (preOrPost_ == "pre" || preOrPost_ == "both") {
    // apply initial relaxation to primary space
    A_->apply (Y, *res1);
    res1->update (ONE, X, -ONE);
    vec1->putScalar (ZERO);
    ifpack2_prec1_->apply (*res1, *vec1);
    Y.update (ONE, *vec1, ONE);
  }

  // project to auxiliary space and smooth
  A_->apply (Y, *res1);
  res1->update (ONE, X, -ONE);
  P_->apply (*res1, *res2, Teuchos::TRANS);
  vec2->putScalar (ZERO);
  ifpack2_prec2_->apply (*res2, *vec2);
  P_->apply (*vec2, *vec1, Teuchos::NO_TRANS);
  Y.update (ONE,*vec1,ONE);

  if (preOrPost_ == "post" || preOrPost_ == "both") {
    // smooth again on primary space
    A_->apply (Y, *res1);
    res1->update (ONE, X, -ONE);
    vec1->putScalar (ZERO);
    ifpack2_prec1_->apply (*res1, *vec1);
    Y.update (ONE, *vec1, ONE);
  }
}

void
OverlappingRowMatrix<MatrixType>::
exportMultiVector (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &OvX,
                   Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> &X,
                   Tpetra::CombineMode CM)
{
  X.doExport (OvX, *Importer_, CM);
}

void
Chebyshev<MatrixType>::
applyMat (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) const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
    "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");

  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
  TEUCHOS_TEST_FOR_EXCEPTION(
    A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::applyMat: The input "
    "matrix A is null.  Please call setMatrix() with a nonnull input matrix "
    "before calling this method.");

  A->apply (X, Y, mode);
}

 void
 apply(
     const Tpetra::MultiVector<double,int,int> & X,
     Tpetra::MultiVector<double,int,int> & Y,
     Teuchos::ETransp mode = Teuchos::NO_TRANS,
     double alpha = Teuchos::ScalarTraits<double>::one(),
     double beta = Teuchos::ScalarTraits<double>::zero()
     ) const
 {
   for (size_t k = 0; k < Y.getNumVectors(); k++) {
     const auto x_data = X.getData(k);
     const auto x0_data = x0_.getData();
     auto y_data = Y.getDataNonConst(k);
     for (size_t i = 0; i < y_data.size(); i++) {
       y_data[i] = 2 * x0_data[i] * x_data[i];
     }
   }
   return;
 }

void DiagonalFilter<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
{
  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
  A_->apply(X,Y,mode,alpha,beta);
  Y.elementWiseMultiply(one,*val_,X,one);
}

void IdentitySolver<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;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
                              global_ordinal_type, node_type> MV;

  TEUCHOS_TEST_FOR_EXCEPTION(
    ! isComputed (), std::runtime_error,
    "Ifpack2::IdentitySolver::apply: 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.");

  // "Identity solver" does what it says: it's the identity operator.
  // We have to Export if the domain and range Maps are not the same.
  // Otherwise, this operator would be a permutation, not the identity.
  if (export_.is_null ()) {
    Y.update (alpha, X, beta);
  }
  else {
    if (alpha == STS::one () && beta == STS::zero ()) { // the common case
      Y.doExport (X, *export_, Tpetra::REPLACE);
    }
    else {
      // We know that the domain and range Maps are compatible.  First
      // bring X into the range Map via Export.  Then compute in place
      // in Y.
      MV X_tmp (Y.getMap (), Y.getNumVectors ());
      X_tmp.doExport (X, *export_, Tpetra::REPLACE);
      Y.update (alpha, X_tmp, beta);
    }
  }
  ++numApply_;
}

void
Chebyshev<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
{
  const std::string timerName ("Ifpack2::Chebyshev::apply");
  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
  if (timer.is_null ()) {
    timer = Teuchos::TimeMonitor::getNewCounter (timerName);
  }

  // Start timing here.
  {
    Teuchos::TimeMonitor timeMon (*timer);

    // compute() calls initialize() if it hasn't already been called.
    // Thus, we only need to check isComputed().
    TEUCHOS_TEST_FOR_EXCEPTION(
      ! isComputed (), std::runtime_error,
      "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
      "you may call apply().");
    TEUCHOS_TEST_FOR_EXCEPTION(
      X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
      "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
      "columns.  X.getNumVectors() = " << X.getNumVectors() << " != "
      << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
    applyImpl (X, Y, mode, alpha, beta);
  }
  ++NumApply_;

  // timer->totalElapsedTime() returns the total time over all timer
  // calls.  Thus, we use = instead of +=.
  ApplyTime_ = timer->totalElapsedTime ();
}

void SparseContainer<MatrixType,InverseType>::
applyImpl (const Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& X,
           Tpetra::MultiVector<InverseScalar,InverseLocalOrdinal,InverseGlobalOrdinal,InverseNode>& Y,
           Teuchos::ETransp mode,
           InverseScalar alpha,
           InverseScalar beta) const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    Inverse_->getDomainMap ()->getNodeNumElements () != X.getLocalLength (),
    std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
    "operator and X have incompatible dimensions (" <<
    Inverse_->getDomainMap ()->getNodeNumElements () << " resp. "
    << X.getLocalLength () << ").  Please report this bug to "
    "the Ifpack2 developers.");
  TEUCHOS_TEST_FOR_EXCEPTION(
    Inverse_->getRangeMap ()->getNodeNumElements () != Y.getLocalLength (),
    std::logic_error, "Ifpack2::SparseContainer::apply: Inverse_ "
    "operator and Y have incompatible dimensions (" <<
    Inverse_->getRangeMap ()->getNodeNumElements () << " resp. "
    << Y.getLocalLength () << ").  Please report this bug to "
    "the Ifpack2 developers.");

  Inverse_->apply (X, Y, mode, alpha, beta);
}

void Diagonal<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
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    ! isComputed (), std::runtime_error, "Ifpack2::Diagonal::apply: You "
    "must first call compute() before you may call apply().  Once you have "
    "called compute(), you need not call it again unless the values in the "
    "matrix have changed, or unless you have called setMatrix().");

  Y.elementWiseMultiply (alpha, *inverseDiag_, X, beta);
  ++numApply_;
}