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

void
LOCA::MultiContinuation::ConstrainedGroup::fillB(
	                                  NOX::Abstract::MultiVector& B) const
{
  std::string callingFunction = 
    "LOCA::MultiContinuation::ConstrainedGroup::fillB";

  bool isZeroB = constraintsPtr->isDXZero();
  Teuchos::RCP<const NOX::Abstract::MultiVector> my_B;

  if (!isZeroB) {
    Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterfaceMVDX> constraints_mvdx = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(constraintsPtr);
    if (constraints_mvdx == Teuchos::null)
      globalData->locaErrorCheck->throwError(
				callingFunction,
				std::string("Constraints object must be of type") +
				std::string("ConstraintInterfaceMVDX"));

    my_B = Teuchos::rcp(constraints_mvdx->getDX(),false);
  }

  // If the underlying system isn't bordered, we're done
  if (!isBordered) {
    if (isZeroB)
      B.init(0.0);
    else
      B = *my_B;
    return;
  }

  // Create views for underlying group
  int w = bordered_grp->getBorderedWidth();
  std::vector<int> idx1(w);
  for (int i=0; i<w; i++)
    idx1[i] = i;
  Teuchos::RCP<NOX::Abstract::MultiVector> underlyingB = 
    B.subView(idx1);

  // Combine blocks in underlying group
  bordered_grp->fillB(*underlyingB);

  // Create views for my blocks
  std::vector<int> idx2(numParams);
  for (int i=0; i<numParams; i++)
    idx2[i] = w+i;
  Teuchos::RCP<NOX::Abstract::MultiVector> my_B_x = 
    B.subView(idx2);

  // Extract solution component from my_B and store in B
  if (isZeroB)
    my_B_x->init(0.0);
  else
    bordered_grp->extractSolutionComponent(*my_B, *my_B_x);
}

NOX::Abstract::Group::ReturnType
LOCA::Thyra::Group::computeDfDpMulti(const std::vector<int>& paramIDs,
                     NOX::Abstract::MultiVector& fdfdp,
                     bool isValidF)
{
  // Currently this does not work because the thyra modelevaluator is not
  // setting the parameter names correctly in the epetraext modelevalator,
  // so we are disabling this for now
  implement_dfdp = false;

  // Use default implementation if we don't want to use model evaluator, or
  // it doesn't support it
  if (!implement_dfdp ||
      !out_args_.supports(::Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,
              param_index).supports(::Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL)) {
    NOX::Abstract::Group::ReturnType res =
      LOCA::Abstract::Group::computeDfDpMulti(paramIDs, fdfdp, isValidF);
    return res;
  }

  // Split fdfdp into f and df/dp
  int num_vecs = fdfdp.numVectors()-1;
  std::vector<int> index_dfdp(num_vecs);
  for (int i=0; i<num_vecs; i++)
    index_dfdp[i] = i+1;
  NOX::Thyra::Vector& f = dynamic_cast<NOX::Thyra::Vector&>(fdfdp[0]);
  Teuchos::RCP<NOX::Abstract::MultiVector> dfdp =
    fdfdp.subView(index_dfdp);

  // Right now this isn't very efficient because we have to compute
  // derivatives with respect to all of the parameters, not just
  // paramIDs.  Will have to work out with Ross how to selectively get
  // parameter derivatives
  int np = params.length();
  Teuchos::RCP<NOX::Thyra::MultiVector> dfdp_full =
    Teuchos::rcp_dynamic_cast<NOX::Thyra::MultiVector>(dfdp->clone(np));

  ::Thyra::ModelEvaluatorBase::DerivativeMultiVector<double> dmv(dfdp_full->getThyraMultiVector(), ::Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL);
  ::Thyra::ModelEvaluatorBase::Derivative<double> deriv(dmv);

  in_args_.set_x(x_vec_->getThyraRCPVector().assert_not_null());
  if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_x_dot))
    in_args_.set_x_dot(x_dot_vec);
  in_args_.set_p(param_index, param_thyra_vec);
  if (!isValidF)
    out_args_.set_f(f.getThyraRCPVector().assert_not_null());
  out_args_.set_DfDp(param_index, deriv);

  // Evaluate model
  model_->evalModel(in_args_, out_args_);

  // Copy back dfdp
  for (int i=0; i<num_vecs; i++)
    (*dfdp)[i] = (*dfdp_full)[paramIDs[i]];

  // Reset inargs/outargs
  in_args_.set_x(Teuchos::null);
  in_args_.set_p(param_index, Teuchos::null);
  out_args_.set_f(Teuchos::null);
  out_args_.set_DfDp(param_index,
             ::Thyra::ModelEvaluatorBase::Derivative<double>());

  if (out_args_.isFailed())
    return NOX::Abstract::Group::Failed;

  return NOX::Abstract::Group::Ok;
}

// Solves turning point equations via classic Salinger bordering
// The first m columns of input_x and input_null store the RHS while
// the last column stores df/dp, d(Jn)/dp respectively.  Note however
// input_param has only m columns (not m+1).  result_x, result_null,
// are result_param have the same dimensions as their input counterparts
NOX::Abstract::Group::ReturnType 
LOCA::TurningPoint::MooreSpence::PhippsBordering::solveTransposeContiguous(
		  Teuchos::ParameterList& params,
		  const NOX::Abstract::MultiVector& input_x,
		  const NOX::Abstract::MultiVector& input_null,
	          const NOX::Abstract::MultiVector::DenseMatrix& input_param,
		  NOX::Abstract::MultiVector& result_x,
		  NOX::Abstract::MultiVector& result_null,
	          NOX::Abstract::MultiVector::DenseMatrix& result_param) const
{
  std::string callingFunction = 
    "LOCA::TurningPoint::MooreSpence::PhippsBordering::solveTransposeContiguous()";
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
  NOX::Abstract::Group::ReturnType status;

  int m = input_x.numVectors()-2;
  std::vector<int> index_input(m);
  std::vector<int> index_input_dp(m+1);
  std::vector<int> index_null(1);
  std::vector<int> index_dp(1);
  for (int i=0; i<m; i++) {
    index_input[i] = i;
    index_input_dp[i] = i;
  }
  index_input_dp[m] = m;
  index_dp[0] = m;
  index_null[0] = m+1;

  NOX::Abstract::MultiVector::DenseMatrix tmp_mat_1(1, m+1);
  NOX::Abstract::MultiVector::DenseMatrix tmp_mat_2(1, m+2);

  // Create view of first m+1 columns of input_null, result_null
  Teuchos::RCP<NOX::Abstract::MultiVector> input_null_view = 
      input_null.subView(index_input_dp);
  Teuchos::RCP<NOX::Abstract::MultiVector> result_null_view = 
      result_null.subView(index_input_dp);

  // verify underlying Jacobian is valid
  if (!group->isJacobian()) {
    status = group->computeJacobian();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Solve  |J^T v||A B| = |G -phi|
  //        |u^T 0||a b|   |0   0 |
  status =
    transposeBorderedSolver->applyInverseTranspose(params, 
						   input_null_view.get(), 
						   NULL, 
						   *result_null_view, 
						   tmp_mat_1);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
							   callingFunction);
  Teuchos::RCP<NOX::Abstract::MultiVector> A = 
    result_null.subView(index_input);
  Teuchos::RCP<NOX::Abstract::MultiVector> B = 
    result_null.subView(index_dp);
  double b = tmp_mat_1(0,m);

  // compute (Jv)_x^T[A B u]
  result_null[m+1] = *uVector;
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp = 
    result_null.clone(NOX::ShapeCopy);
  status = group->computeDwtJnDxMulti(result_null, *nullVector, *tmp);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							   finalStatus,
							   callingFunction);

  // compute [F 0 0] - (Jv)_x^T[A B u]
  tmp->update(1.0, input_x, -1.0);

  // verify underlying Jacobian is valid
  if (!group->isJacobian()) {
    status = group->computeJacobian();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Solve  |J^T v||C D E| = |F - (Jv)_x^T A  -(Jv)_x^T B  -(Jv)_x^T u|
  //        |u^T 0||c d e|   |         0             0            0   |
  status = 
    transposeBorderedSolver->applyInverseTranspose(params, 
						   tmp.get(), 
						   NULL, 
						   result_x,
						   tmp_mat_2);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
//.........这里部分代码省略.........

// Solves Hopf equations via classic Salinger bordering
// The first m columns of input_x, input_y, input_z store the RHS, the
// next column stores df/dp, (Jy-wBz)_p and (Jz+wBy)_p respectively, the
// last column of input_y and input_z store Bz and -By respectively.  Note 
// input_x has m+1 columns, input_y and input_z have m+2, and input_w and
// input_p have m columns.  result_x, result_y, result_z, result_w and 
// result_param have the same dimensions as their input counterparts
NOX::Abstract::Group::ReturnType 
LOCA::Hopf::MooreSpence::SalingerBordering::solveContiguous(
		      Teuchos::ParameterList& params,
		      const NOX::Abstract::MultiVector& input_x,
		      const NOX::Abstract::MultiVector& input_y,
		      const NOX::Abstract::MultiVector& input_z,
		      const NOX::Abstract::MultiVector::DenseMatrix& input_w,
		      const NOX::Abstract::MultiVector::DenseMatrix& input_p,
		      NOX::Abstract::MultiVector& result_x,
		      NOX::Abstract::MultiVector& result_y,
		      NOX::Abstract::MultiVector& result_z,
		      NOX::Abstract::MultiVector::DenseMatrix& result_w,
	              NOX::Abstract::MultiVector::DenseMatrix& result_p) const
{
  std::string callingFunction = 
    "LOCA::Hopf::MooreSpence::SalingerBordering::solveContiguous()";
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
  NOX::Abstract::Group::ReturnType status;

  int m = input_x.numVectors()-1;
  std::vector<int> index_input(m);
  std::vector<int> index_dp(1);
  std::vector<int> index_B(1);
  std::vector<int> index_ip(m+1);
  for (int i=0; i<m; i++) {
    index_input[i] = i;
    index_ip[i] = i;
  }
  index_ip[m] = m;
  index_dp[0] = m;
  index_B[0] = m+1;

  // verify underlying Jacobian is valid
  if (!group->isJacobian()) {
    status = group->computeJacobian();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }
  
  // compute [A b] = J^-1 [F df/dp]
  status = group->applyJacobianInverseMultiVector(params, input_x, result_x);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
							   callingFunction);
  Teuchos::RCP<NOX::Abstract::MultiVector> A = 
    result_x.subView(index_input);
  Teuchos::RCP<NOX::Abstract::MultiVector> b = 
    result_x.subView(index_dp);

  // verify underlying complex matrix is valid
   if (!group->isComplex()) {
    status = group->computeComplex(w);
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // compute (J+iwB)(y+iz)_x [A b]
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp_real = 
    result_y.clone(NOX::ShapeCopy);
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp_real_sub =
    tmp_real->subView(index_ip);
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp_imag = 
    result_y.clone(NOX::ShapeCopy);
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp_imag_sub =
    tmp_imag->subView(index_ip);
  tmp_real->init(0.0);
  tmp_imag->init(0.0);
  status = group->computeDCeDxa(*yVector, *zVector, w, result_x,
				*CeRealVector, *CeImagVector, *tmp_real_sub,
				*tmp_imag_sub);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
							   callingFunction);

  // compute [G+iH d(J+iwB)(y+iz)/dp iB(y+iz)] - [(J+iwB)_x[A b] 0+i0]
  tmp_real->update(1.0, input_y, -1.0);
  tmp_imag->update(1.0, input_z, -1.0);

  // verify underlying complex matrix is valid
  if (!group->isComplex()) {
    status = group->computeComplex(w);
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // compute [C+iD e+if g+ih] = (J+iwB)^-1 (tmp_real + i tmp_imag)
  status = group->applyComplexInverseMultiVector(params, *tmp_real, *tmp_imag,
//.........这里部分代码省略.........

void
LOCA::BorderedSolver::HouseholderQR::computeQR(
                const NOX::Abstract::MultiVector::DenseMatrix& C,
                const NOX::Abstract::MultiVector& B,
                bool use_c_transpose,
                NOX::Abstract::MultiVector::DenseMatrix& Y1,
                NOX::Abstract::MultiVector& Y2,
                NOX::Abstract::MultiVector::DenseMatrix& T,
                NOX::Abstract::MultiVector::DenseMatrix& R)
{
  double beta;
  int m = B.numVectors();

  // Initialize
  Y1.putScalar(0.0);
  T.putScalar(0.0);
  Y2 = B;
  if (use_c_transpose) {
    for (int i=0; i<m; i++)
      for (int j=0; j<m; j++)
    R(i,j) = C(j,i);        // Copy transpose of C into R
  }
  else
    R.assign(C);

  // A temporary vector
  Teuchos::RCP<NOX::Abstract::MultiVector> v2 = Y2.clone(1);

  Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> v1;
  Teuchos::RCP<NOX::Abstract::MultiVector> h2;
  Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> h1;
  Teuchos::RCP<NOX::Abstract::MultiVector> y2;
  Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> y1;
  Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> z;
  std::vector<int> h_idx;
  std::vector<int> y_idx;
  y_idx.reserve(m);

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

    // Create view of column i of Y1 starting at row i
    v1 =
      Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
                                   Y1,
                                   m-i,
                                   1, i, i));

    // Create view of columns i through m-1 of Y2
    h_idx.resize(m-i);
    for (unsigned int j=0; j<h_idx.size(); j++)
      h_idx[j] = i+j;
    h2 = Y2.subView(h_idx);

    // Create view of columns i thru m-1 of R, starting at row i
    h1 =
      Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
                                   R,
                                   m-i,
                                   m-i,
                                   i, i));

    if (i > 0) {

      // Create view of columns 0 through i-1 of Y2
      y_idx.push_back(i-1);
      y2 = Y2.subView(y_idx);

      // Create view of columns 0 through i-1 of Y1, starting at row i
      y1 =
    Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
                                 Y1,
                                 m-i,
                                 i, i, 0));

      // Create view of column i, row 0 through i-1 of T
      z =
    Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
                                 T,
                                 i,
                                 1,
                                 0, i));
    }

    // Compute Householder Vector
    computeHouseholderVector(i, R, Y2, *v1, *v2, beta);

    // Apply Householder reflection
    applyHouseholderVector(*v1, *v2, beta, *h1, *h2);

    // Copy v2 into Y2
    Y2[i] = (*v2)[0];

    T(i,i) = -beta;

    if (i > 0) {

      // Compute z = y2^T * v2
      v2->multiply(1.0, *y2, *z);

      // Compute z = -beta * (z + y1^T * v1)
//.........这里部分代码省略.........