C++ MultiVectorBase::domain方法代码示例

void DefaultMultiVectorLinearOpWithSolve<Scalar>::applyImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &XX,
  const Ptr<MultiVectorBase<Scalar> > &YY,
  const Scalar alpha,
  const Scalar beta
  ) const
{

  using Teuchos::dyn_cast;
  typedef DefaultMultiVectorProductVector<Scalar> MVPV;

  const Ordinal numCols = XX.domain()->dim();

  for (Ordinal col_j = 0; col_j < numCols; ++col_j) {

    const RCP<const VectorBase<Scalar> > x = XX.col(col_j);
    const RCP<VectorBase<Scalar> > y = YY->col(col_j);
 
    RCP<const MultiVectorBase<Scalar> >
      X = dyn_cast<const MVPV>(*x).getMultiVector().assert_not_null();
    RCP<MultiVectorBase<Scalar> >
      Y = dyn_cast<MVPV>(*y).getNonconstMultiVector().assert_not_null();
    
    Thyra::apply( *lows_.getConstObj(), M_trans, *X, Y.ptr(), alpha, beta );

  }

}

THYRA_DEPRECATED
SolveStatus<Scalar>
solveTranspose(
  const LinearOpWithSolveBase<Scalar> &A,
  const EConj conj,
  const MultiVectorBase<Scalar> &B,
  MultiVectorBase<Scalar> *X,
  const SolveCriteria<Scalar> *solveCriteria = NULL
  )
{
  typedef SolveCriteria<Scalar> SC;
  typedef BlockSolveCriteria<Scalar> BSC;
  typedef SolveStatus<Scalar> BSS;
  SC defaultSolveCriteria;
  BSC blockSolveCriteria[1];
  BSS blockSolveStatus[1];
  blockSolveCriteria[0] = BSC(
    solveCriteria ? *solveCriteria : defaultSolveCriteria,
    B.domain()->dim());
  A.solveTranspose(
    conj,B,X,1,
    blockSolveCriteria,
    blockSolveStatus
    );
  return blockSolveStatus[0];
}

Teuchos::Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>
Thyra::norms_inf( const MultiVectorBase<Scalar>& V )
{
  typedef typename ScalarTraits<Scalar>::magnitudeType ScalarMag;
  Array<ScalarMag> norms(V.domain()->dim());
  Thyra::norms_inf<Scalar>(V, norms());
  return norms;
}

void DefaultMultipliedLinearOp<Scalar>::applyImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &X,
  const Ptr<MultiVectorBase<Scalar> > &Y,
  const Scalar alpha,
  const Scalar beta
  ) const
{
  using Teuchos::rcpFromPtr;
  using Teuchos::rcpFromRef;
#ifdef TEUCHOS_DEBUG
  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
    "DefaultMultipliedLinearOp<Scalar>::apply(...)", *this, M_trans, X, &*Y
    );
#endif // TEUCHOS_DEBUG  
  const int nOps = Ops_.size();
  const Ordinal m = X.domain()->dim();
  if( real_trans(M_trans)==NOTRANS ) {
    //
    // Y = alpha * M * X + beta*Y
    // =>
    // Y = alpha * op(Op[0]) * op(Op[1]) * ... * op(Op[numOps-1]) * X + beta*Y
    //
    RCP<MultiVectorBase<Scalar> > T_kp1, T_k; // Temporary propagated between loops 
    for( int k = nOps-1; k >= 0; --k ) {
      RCP<MultiVectorBase<Scalar> > Y_k;
      RCP<const MultiVectorBase<Scalar> > X_k;
      if(k==0) Y_k = rcpFromPtr(Y);  else Y_k = T_k = createMembers(getOp(k)->range(), m);
      if(k==nOps-1) X_k = rcpFromRef(X); else X_k = T_kp1;
      if( k > 0 )
        Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
      else
        Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
      T_kp1 = T_k;
    }
  }
  else {
    //
    // Y = alpha * M' * X + beta*Y
    // =>
    // Y = alpha * Op[numOps-1]' * Op[1]' * ... * Op[0]' * X + beta * Y
    //
    RCP<MultiVectorBase<Scalar> > T_km1, T_k; // Temporary propagated between loops 
    for( int k = 0; k <= nOps-1; ++k ) {
      RCP<MultiVectorBase<Scalar> >         Y_k;
      RCP<const MultiVectorBase<Scalar> >   X_k;
      if(k==nOps-1) Y_k = rcpFromPtr(Y);  else Y_k = T_k = createMembers(getOp(k)->domain(), m);
      if(k==0) X_k = rcpFromRef(X); else X_k = T_km1;
      if( k < nOps-1 )
        Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr());
      else
        Thyra::apply(*getOp(k), M_trans, *X_k, Y_k.ptr(), alpha, beta);
      T_km1 = T_k;
    }
  }
}

void LinearOpScalarProd<Scalar>::scalarProdsImpl(
  const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y,
  const ArrayView<Scalar> &scalarProds_out
  ) const
{
  Teuchos::RCP<MultiVectorBase<Scalar> >
    T = createMembers(Y.range() ,Y.domain()->dim());
  Thyra::apply(*op_, NOTRANS,Y, T.ptr());
  dots(X, *T, scalarProds_out);
}

void DefaultColumnwiseMultiVector<Scalar>::applyImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &X,
  const Ptr<MultiVectorBase<Scalar> > &Y,
  const Scalar alpha,
  const Scalar beta
  ) const
{
#ifdef TEUCHOS_DEBUG
  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
    "MultiVectorBase<Scalar>::apply()", *this, M_trans, X, &*Y);
#endif
  const Ordinal nc = this->domain()->dim();
  const Ordinal m = X.domain()->dim();
  for (Ordinal col_j = 0; col_j < m; ++col_j) {
    const RCP<const VectorBase<Scalar> > x_j = X.col(col_j);
    const RCP<VectorBase<Scalar> > y_j = Y->col(col_j);
    // y_j *= beta
    Vt_S(y_j.ptr(), beta);
    // y_j += alpha*op(M)*x_j
    if(M_trans == NOTRANS) {
      //
      // y_j += alpha*M*x_j = alpha*M.col(0)*x_j(0) + ... + alpha*M.col(nc-1)*x_j(nc-1)
      //
      // Extract an explicit view of x_j
      RTOpPack::ConstSubVectorView<Scalar> x_sub_vec;               
      x_j->acquireDetachedView(Range1D(), &x_sub_vec);
      // Loop through and add the multiple of each column
      for (Ordinal j = 0; j < nc; ++j )
        Vp_StV( y_j.ptr(), Scalar(alpha*x_sub_vec(j)), *this->col(j) );
      // Release the view of x
      x_j->releaseDetachedView(&x_sub_vec);
    }
    else {
      //
      //                        [ alpha*dot(M.col(0),x_j)    ]
      // y_j += alpha*M^T*x_j = [ alpha*dot(M.col(1),x_j)    ]
      //                        [ ...                        ]
      //                        [ alpha*dot(M.col(nc-1),x_j) ]
      //
      // Extract an explicit view of y_j
      RTOpPack::SubVectorView<Scalar> y_sub_vec;               
      y_j->acquireDetachedView(Range1D(), &y_sub_vec);
      // Loop through and add to each element in y_j
      for (Ordinal j = 0; j < nc; ++j )
        y_sub_vec(j) += alpha*dot(*this->col(j), *x_j);
      // Commit explicit view of y
      y_j->commitDetachedView(&y_sub_vec);
    }
  }
}

void Thyra::reductions( const MultiVectorBase<Scalar>& V, const NormOp &op,
  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
{
  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
  const int m = V.domain()->dim();
  Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
  Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
  for( int kc = 0; kc < m; ++kc ) {
    rcp_op_targs[kc] = op.reduct_obj_create();
    op_targs[kc] = rcp_op_targs[kc].ptr();
  }
  applyOp<Scalar>(op, tuple(ptrInArg(V)),
    ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
    op_targs );
  for( int kc = 0; kc < m; ++kc ) {
    norms[kc] = op(*op_targs[kc]);
  }
}

void DefaultDiagonalLinearOp<Scalar>::applyImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &X,
  const Ptr<MultiVectorBase<Scalar> > &Y,
  const Scalar alpha,
  const Scalar beta
  ) const
{
  typedef Teuchos::ScalarTraits<Scalar> ST;

#ifdef TEUCHOS_DEBUG
  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
    "DefaultDiagonalLinearOp<Scalar>::apply(...)",*this, M_trans, X, &*Y
    );
#endif // TEUCHOS_DEBUG

  // Y = beta * Y

  if( beta != ST::one() ) scale<Scalar>(beta, Y);

  // Y += alpha *op(M) * X

  const Ordinal m = X.domain()->dim();

  for (Ordinal col_j = 0; col_j < m; ++col_j) {
    const RCP<const VectorBase<Scalar> > x = X.col(col_j);
    const RCP<VectorBase<Scalar> > y = Y->col(col_j);
    if (ST::isComplex) {
      if ( M_trans==NOTRANS || M_trans==TRANS ) {
        ele_wise_prod( alpha, *diag_.getConstObj(), *x, y.ptr() );
      }
      else {
        ele_wise_conj_prod( alpha, *diag_.getConstObj(), *x, y.ptr() );
      }
    }
    else {
      ele_wise_prod( alpha, *diag_.getConstObj(), *x, y.ptr() );
    }
  }

}

SolveStatus<Scalar>
DefaultMultiVectorLinearOpWithSolve<Scalar>::solveImpl(
  const EOpTransp transp,
  const MultiVectorBase<Scalar> &BB,
  const Ptr<MultiVectorBase<Scalar> > &XX,
  const Ptr<const SolveCriteria<Scalar> > solveCriteria
  ) const
{

  using Teuchos::dyn_cast;
  using Teuchos::outArg;
  using Teuchos::inOutArg;
  typedef DefaultMultiVectorProductVector<Scalar> MVPV;

  const Ordinal numCols = BB.domain()->dim();

  SolveStatus<Scalar> overallSolveStatus;
  accumulateSolveStatusInit(outArg(overallSolveStatus));
  
  for (Ordinal col_j = 0; col_j < numCols; ++col_j) {

    const RCP<const VectorBase<Scalar> > b = BB.col(col_j);
    const RCP<VectorBase<Scalar> > x = XX->col(col_j);

    RCP<const MultiVectorBase<Scalar> >
      B = dyn_cast<const MVPV>(*b).getMultiVector().assert_not_null();
    RCP<MultiVectorBase<Scalar> >
      X = dyn_cast<MVPV>(*x).getNonconstMultiVector().assert_not_null();

    const SolveStatus<Scalar> solveStatus =
      Thyra::solve(*lows_.getConstObj(), transp, *B, X.ptr(), solveCriteria);

    accumulateSolveStatus(
      SolveCriteria<Scalar>(), // Never used
      solveStatus, inOutArg(overallSolveStatus) );

  }
  
  return overallSolveStatus;

}

void doExplicitMultiVectorAdjoint(
  const MultiVectorBase<Scalar>& mvIn, MultiVectorBase<Scalar>* mvTransOut
  )
{
  typedef Teuchos::ScalarTraits<Scalar> ST;
#ifdef TEUCHOS_DEBUG
  TEST_FOR_EXCEPT(0==mvTransOut);
  THYRA_ASSERT_VEC_SPACES("doExplicitMultiVectorAdjoint(...)",
    *mvIn.domain(), *mvTransOut->range()
    );
  THYRA_ASSERT_VEC_SPACES("doExplicitMultiVectorAdjoint(...)",
    *mvIn.range(), *mvTransOut->domain()
    );
#endif
  ConstDetachedMultiVectorView<Scalar> dMvIn(mvIn);
  DetachedMultiVectorView<Scalar> dMvTransOut(*mvTransOut);
  const int m = dMvIn.subDim();
  const int n = dMvIn.numSubCols();
  for ( int j = 0; j < n; ++j ) {
    for ( int i = 0; i < m; ++i ) {
      dMvTransOut(j,i) = ST::conjugate(dMvIn(i,j));
    }
  }
}

SolveStatus<Scalar>
BelosLinearOpWithSolve<Scalar>::solveImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &B,
  const Ptr<MultiVectorBase<Scalar> > &X,
  const Ptr<const SolveCriteria<Scalar> > solveCriteria
  ) const
{

  TEUCHOS_FUNC_TIME_MONITOR("BelosLOWS");

  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using Teuchos::rcpFromPtr;
  using Teuchos::FancyOStream;
  using Teuchos::OSTab;
  using Teuchos::describe;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType ScalarMag;
  Teuchos::Time totalTimer(""), timer("");
  totalTimer.start(true);

  assertSolveSupports(*this, M_trans, solveCriteria);
  // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
  // solve(...).

  const int numRhs = B.domain()->dim();
  const int numEquations = B.range()->dim();

  const RCP<FancyOStream> out = this->getOStream();
  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  OSTab tab = this->getOSTab();
  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
    *out << "\nStarting iterations with Belos:\n";
    OSTab tab2(out);
    *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
    *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
    *out << "With #Eqns="<<numEquations<<", #RHSs="<<numRhs<<" ...\n";
  }

  //
  // Set RHS and LHS
  //

  bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
  TEST_FOR_EXCEPTION(
    ret == false, CatastrophicSolveFailure
    ,"Error, the Belos::LinearProblem could not be set for the current solve!"
    );

  //
  // Set the solution criteria
  //

  const RCP<Teuchos::ParameterList> tmpPL = Teuchos::parameterList();

  SolveMeasureType solveMeasureType;
  RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
  if (nonnull(solveCriteria)) {
    solveMeasureType = solveCriteria->solveMeasureType;
    const ScalarMag requestedTol = solveCriteria->requestedTol;
    if (solveMeasureType.useDefault()) {
      tmpPL->set("Convergence Tolerance", defaultTol_);
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      tmpPL->set("Explicit Residual Scaling", "Norm of RHS");
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      tmpPL->set("Explicit Residual Scaling", "Norm of Initial Residual");
    }
    else {
      // Set the most generic (and inefficient) solve criteria
      generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
        *solveCriteria, convergenceTestFrequency_);
      // Set the verbosity level (one level down)
      generalSolveCriteriaBelosStatusTest->setOStream(out);
      generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
      // Set the default convergence tolerance to always converged to allow
      // the above status test to control things.
      tmpPL->set("Convergence Tolerance", 1.0);
    }
  }
  else {
    // No solveCriteria was even passed in!
    tmpPL->set("Convergence Tolerance", defaultTol_);
  }

  //
//.........这里部分代码省略.........

SolveStatus<Scalar>
BelosLinearOpWithSolve<Scalar>::solveImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &B,
  const Ptr<MultiVectorBase<Scalar> > &X,
  const Ptr<const SolveCriteria<Scalar> > solveCriteria
  ) const
{

  THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");

  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using Teuchos::rcpFromPtr;
  using Teuchos::FancyOStream;
  using Teuchos::OSTab;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::describe;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType ScalarMag;
  Teuchos::Time totalTimer(""), timer("");
  totalTimer.start(true);

  assertSolveSupports(*this, M_trans, solveCriteria);
  // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
  // solve(...).

  const RCP<FancyOStream> out = this->getOStream();
  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  OSTab tab = this->getOSTab();
  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
    *out << "\nStarting iterations with Belos:\n";
    OSTab tab2(out);
    *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
    *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
    *out << "With #Eqns="<<B.range()->dim()<<", #RHSs="<<B.domain()->dim()<<" ...\n";
  }

  //
  // Set RHS and LHS
  //

  bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
  TEUCHOS_TEST_FOR_EXCEPTION(
    ret == false, CatastrophicSolveFailure
    ,"Error, the Belos::LinearProblem could not be set for the current solve!"
    );

  //
  // Set the solution criteria
  //

  // Parameter list for the current solve.
  const RCP<ParameterList> tmpPL = Teuchos::parameterList();

  // The solver's valid parameter list.
  RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();

  SolveMeasureType solveMeasureType;
  RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
  if (nonnull(solveCriteria)) {
    solveMeasureType = solveCriteria->solveMeasureType;
    const ScalarMag requestedTol = solveCriteria->requestedTol;
    if (solveMeasureType.useDefault()) {
      tmpPL->set("Convergence Tolerance", defaultTol_);
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      setResidualScalingType (tmpPL, validPL, "Norm of RHS");
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
    }
    else {
      // Set the most generic (and inefficient) solve criteria
      generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
        *solveCriteria, convergenceTestFrequency_);
      // Set the verbosity level (one level down)
      generalSolveCriteriaBelosStatusTest->setOStream(out);
      generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
      // Set the default convergence tolerance to always converged to allow
      // the above status test to control things.
      tmpPL->set("Convergence Tolerance", 1.0);
    }
    // maximum iterations
    if (nonnull(solveCriteria->extraParameters)) {
      if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,"Maximum Iterations")) {
        tmpPL->set("Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,"Maximum Iterations"));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (nonnull(solveCriteria->extraParameters)) {
      maxIterations = solveCriteria->extraParameters->get("Maximum Iterations",maxIterations);
    }
  }

  //
  // Get Epetra_MultiVector views of B and X
  //

  RCP<const Epetra_MultiVector> epetra_B;
  RCP<Epetra_MultiVector> epetra_X;

  const EpetraOperatorWrapper* opWrapper =
    dynamic_cast<const EpetraOperatorWrapper*>(aztecOp);

  if (opWrapper == 0) {
    epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
    epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
  }

  //
  // Use AztecOO to solve each RHS one at a time (which is all that I can do anyway)
  //

  int totalIterations = 0;
  SolveStatus<double> solveStatus;
  solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
  solveStatus.achievedTol = -1.0;

  /* Get the number of columns in the multivector. We use Thyra
   * functions rather than Epetra functions to do this, as we
   * might not yet have created an Epetra multivector. - KL */
  //const int m = epetra_B->NumVectors();
  const int m = B.domain()->dim();

  for( int j = 0; j < m; ++j ) {

    THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: AztecOOLOWS:SingleSolve", SingleSolve);

    //
    // Get Epetra_Vector views of B(:,j) and X(:,j)
    // How this is done will depend on whether we have a true Epetra operator
    // or we are wrapping a general Thyra operator in an Epetra operator.
    //

    // We need to declare epetra_x_j as non-const because when we have a phony
    // Epetra operator we'll have to copy a thyra vector into it.
    RCP<Epetra_Vector> epetra_b_j;
    RCP<Epetra_Vector> epetra_x_j;

    if (opWrapper == 0) {
      epetra_b_j = rcpFromRef(*const_cast<Epetra_Vector*>((*epetra_B)(j)));
      epetra_x_j = rcpFromRef(*(*epetra_X)(j));
    }
    else {
      if (is_null(epetra_b_j)) {
        epetra_b_j = rcp(new Epetra_Vector(opRangeMap));
        epetra_x_j = rcp(new Epetra_Vector(opDomainMap));
      }
      opWrapper->copyThyraIntoEpetra(*B.col(j), *epetra_b_j);
      opWrapper->copyThyraIntoEpetra(*X->col(j), *epetra_x_j);
    }

    //
    // Set the RHS and LHS
    //

void EpetraLinearOp::applyImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<double> &X_in,
  const Ptr<MultiVectorBase<double> > &Y_inout,
  const double alpha,
  const double beta
  ) const
{

  THYRA_FUNC_TIME_MONITOR("Thyra::EpetraLinearOp::euclideanApply");

  const EOpTransp real_M_trans = real_trans(M_trans);

#ifdef TEUCHOS_DEBUG
  TEUCHOS_TEST_FOR_EXCEPT(!isFullyInitialized_);
  THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
    "EpetraLinearOp::euclideanApply(...)", *this, M_trans, X_in, Y_inout
    );
  TEUCHOS_TEST_FOR_EXCEPTION(
    real_M_trans==TRANS && adjointSupport_==EPETRA_OP_ADJOINT_UNSUPPORTED,
    Exceptions::OpNotSupported,
    "EpetraLinearOp::apply(...): *this was informed that adjoints "
    "are not supported when initialized." 
    );
#endif

  const RCP<const VectorSpaceBase<double> > XY_domain = X_in.domain();
  const int numCols = XY_domain->dim();
 
  //
  // Get Epetra_MultiVector objects for the arguments
  //
  // 2007/08/18: rabartl: Note: After profiling, I found that calling the more
  // general functions get_Epetra_MultiVector(...) was too slow. These
  // functions must ensure that memory is being remembered efficiently and the
  // use of extra data with the RCP and other things is slow.
  //
  RCP<const Epetra_MultiVector> X;
  RCP<Epetra_MultiVector> Y;
  {
    THYRA_FUNC_TIME_MONITOR_DIFF(
      "Thyra::EpetraLinearOp::euclideanApply: Convert MultiVectors", MultiVectors);
    // X
    X = get_Epetra_MultiVector(
      real_M_trans==NOTRANS ? getDomainMap() : getRangeMap(), X_in );
    // Y
    if( beta == 0 ) {
      Y = get_Epetra_MultiVector(
        real_M_trans==NOTRANS ? getRangeMap() : getDomainMap(), *Y_inout );
    }
  }

  //
  // Set the operator mode
  //

  /* We need to save the transpose state here, and then reset it after 
   * application. The reason for this is that if we later apply the 
   * operator outside Thyra (in Aztec, for instance), it will remember
   * the transpose flag set here. */
  bool oldState = op_->UseTranspose();
  op_->SetUseTranspose(
    real_trans(trans_trans(opTrans_,M_trans)) == NOTRANS ? false : true );

  //
  // Perform the apply operation
  //
  {
    THYRA_FUNC_TIME_MONITOR_DIFF(
      "Thyra::EpetraLinearOp::euclideanApply: Apply", Apply);
    if( beta == 0.0 ) {
      // Y = M * X
      if( applyAs_ == EPETRA_OP_APPLY_APPLY ) {
        THYRA_FUNC_TIME_MONITOR_DIFF(
          "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Apply",
          ApplyApply);
        op_->Apply( *X, *Y );
      }
      else if( applyAs_ == EPETRA_OP_APPLY_APPLY_INVERSE ) {
        THYRA_FUNC_TIME_MONITOR_DIFF(
          "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): ApplyInverse",
          ApplyApplyInverse);
        op_->ApplyInverse( *X, *Y );
      }
      else {
#ifdef TEUCHOS_DEBUG
        TEUCHOS_TEST_FOR_EXCEPT(true);
#endif
      }
      // Y = alpha * Y
      if( alpha != 1.0 ) {
        THYRA_FUNC_TIME_MONITOR_DIFF(
          "Thyra::EpetraLinearOp::euclideanApply: Apply(beta==0): Scale Y",
          Scale);
        Y->Scale(alpha);
      }
    }
    else {  // beta != 0.0
      // Y_inout = beta * Y_inout
      if(beta != 0.0) {
//.........这里部分代码省略.........