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

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 );

  }

}

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

}

//.........这里部分代码省略.........
  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
    //

    aztecSolver->SetRHS(&*epetra_b_j);
    aztecSolver->SetLHS(&*epetra_x_j);

    //
    // Solve the linear system
    //
    timer.start(true);
    {
      SetAztecSolveState
        setAztecSolveState(aztecSolver,out,verbLevel,solveMeasureType);
      aztecSolver->Iterate( maxIterations, tol );
      // NOTE: We ignore the returned status but get it below
    }
    timer.stop();

    //
    // Scale the solution
    // (Originally, this was at the end of the loop after all columns had been
    // processed. It's moved here because we need to do it before copying the
    // solution back into a Thyra vector. - KL
    //
    if (aztecSolverScalar_ != 1.0)
      epetra_x_j->Scale(1.0/aztecSolverScalar_);

    //