C++ multivector::DenseMatrix类代码示例

void 
LOCA::Extended::MultiVector::multiply(
			    double alpha, 
			    const LOCA::Extended::MultiVector& y,
			    NOX::Abstract::MultiVector::DenseMatrix& b) const
{
  // Verify dimensions are consistent
  if (y.numMultiVecRows != numMultiVecRows || y.numColumns != b.numRows() ||
      y.numScalarRows != numScalarRows || numColumns != b.numCols()) 
    globalData->locaErrorCheck->throwError(
  "LOCA::Extended::MultiVector::multiply()",
  "Size of supplied multivector/matrix is incompatible with this multivector");

  // Zero out b
  b.putScalar(0.0);

  // Create temporary matrix to hold product for each multivec
  NOX::Abstract::MultiVector::DenseMatrix tmp(b);

  // Compute and sum products for each multivec
  for (int i=0; i<numMultiVecRows; i++) {
    multiVectorPtrs[i]->multiply(alpha, *(y.multiVectorPtrs[i]), tmp);
    b += tmp;
  }

  // Compute and add in product for scalars
  if (numScalarRows > 0)
    b.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, alpha, *y.scalarsPtr,
	       *scalarsPtr, 1.0);

}

NOX::Abstract::Group::ReturnType
LOCA::MultiContinuation::CompositeConstraint::multiplyDX(
		      double alpha, 
		      const NOX::Abstract::MultiVector& input_x,
	              NOX::Abstract::MultiVector::DenseMatrix& result_p) const
{
  std::string callingFunction = 
    "LOCA::MultiContinuation::CompositeConstraint::multiplyDX()";
  NOX::Abstract::Group::ReturnType status;
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;

  // If dg/dx is zero for every constraint, result_p is zero
  if (isDXZero()) {
    result_p.putScalar(0.0);
    return finalStatus;
  }

  Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_p_sub;
  int num_rows;
  int num_cols = result_p.numCols();
  for (int i=0; i<numConstraintObjects; i++) {

    num_rows = constraintPtrs[i]->numConstraints();

    // if dg/dx is zero for this constraint, set corresponding entries of
    // result_p to zero
    if (constraintPtrs[i]->isDXZero()) {
      for (int j=0; j<num_rows; j++)
	for (int k=0; k<num_cols; k++)
	  result_p(indices[i][j],k) = 0.0;
    }
    else {

      // Create a sub view of rows indices[i][0] -- indices[i][end] 
      // of result_p 
      result_p_sub = 
	Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
								 result_p,
								 num_rows,
								 num_cols,
								 indices[i][0],
								 0));

      status = constraintPtrs[i]->multiplyDX(alpha, input_x, 
					     *result_p_sub);
      finalStatus = 
	globalData->locaErrorCheck->combineAndCheckReturnTypes(
							     status, 
							     finalStatus,
							     callingFunction);
    }

  }

  return finalStatus;
}

void
LOCA::BorderedSolver::HouseholderQR::applyCompactWY(
                const NOX::Abstract::MultiVector::DenseMatrix& Y1,
                const NOX::Abstract::MultiVector& Y2,
                const NOX::Abstract::MultiVector::DenseMatrix& T,
                NOX::Abstract::MultiVector::DenseMatrix& X1,
                NOX::Abstract::MultiVector& X2,
                bool isZeroX1, bool isZeroX2,
                bool useTranspose) const
{
  if (isZeroX1 && isZeroX2) {
    X1.putScalar(0.0);
    X2.init(0.0);
    return;
  }

  int m = Y2.numVectors();

  Teuchos::ETransp T_flag;
  if (useTranspose)
    T_flag = Teuchos::TRANS;
  else
    T_flag = Teuchos::NO_TRANS;

  NOX::Abstract::MultiVector::DenseMatrix tmp(m, X2.numVectors());

  // Compute Y1^T*X1 + Y2^T*X2
  if (!isZeroX2)
    X2.multiply(1.0, Y2, tmp);

  // Opportunity for optimization here since Y1 is a lower-triangular
  // matrix with unit diagonal
  if (!isZeroX2 && !isZeroX1)
    tmp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Y1, X1, 1.0);
  else if (!isZeroX1)
    tmp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Y1, X1, 0.0);

  // Compute op(T)*(Y1^T*X1 + Y2^T*X2)
  dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, T_flag,
         Teuchos::NON_UNIT_DIAG, tmp.numRows(), tmp.numCols(), 1.0,
         T.values(), T.numRows(), tmp.values(), tmp.numRows());

  // Compute X1 = X1 + Y1*op(T)*(Y1^T*X1 + Y2^T*X2)
  // Opportunity for optimization here since Y1 is a lower-triangular
  // matrix with unit diagonal
  if (isZeroX1)
    X1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Y1, tmp, 0.0);
  else
    X1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Y1, tmp, 1.0);

  // Compute X2 = X2 + Y1*op(T)*(Y1^T*X1 + Y2^T*X2)
  if (isZeroX2)
    X2.update(Teuchos::NO_TRANS, 1.0, Y2, tmp, 0.0);
  else
    X2.update(Teuchos::NO_TRANS, 1.0, Y2, tmp, 1.0);
}

NOX::Abstract::Group::ReturnType
LOCA::TurningPoint::MinimallyAugmented::Constraint::
computeDP(const std::vector<int>& paramIDs, 
	  NOX::Abstract::MultiVector::DenseMatrix& dgdp, 
	  bool isValidG)
{
  std::string callingFunction = 
    "LOCA::TurningPoint::MinimallyAugmented::Constraint::computeDP()";
  NOX::Abstract::Group::ReturnType status;
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;

  // Compute sigma, w and v if necessary
  if (!isValidConstraints) {
    status = computeConstraints();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Compute -(w^T*J*v)_p
  status = grpPtr->computeDwtJnDp(paramIDs, (*w_vector)[0], (*v_vector)[0], 
				  dgdp, false);
  finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  dgdp.scale(-1.0/sigma_scale);

  // Set the first column of dgdp
  dgdp(0,0) = constraints(0,0);

  return finalStatus;
}

int 
NOX::TestCompare::testMatrix(
		 const NOX::Abstract::MultiVector::DenseMatrix& mat, 
		 const NOX::Abstract::MultiVector::DenseMatrix& mat_expected, 
		 double rtol, double atol, 
		 const std::string& name)
{
  bool passed;

  NOX::Abstract::MultiVector::DenseMatrix tmp(mat_expected.numRows(),
					      mat_expected.numCols());

  for (int j=0; j<mat_expected.numCols(); j++)
    for (int i=0; i<mat_expected.numRows(); i++)
      tmp(i,j) = fabs(mat(i,j)-mat_expected(i,j)) / 
	(atol + rtol * fabs(mat_expected(i,j)));
 
  double inf_norm = tmp.normInf();

  if (inf_norm < 1)
    passed = true;
  else
    passed = false;

  if (utils.isPrintType(NOX::Utils::TestDetails)) {
    os << std::endl
	 << "\tChecking " << name << ":  ";
    if (passed)
      os << "Passed." << std::endl;
    else
      os << "Failed." << std::endl;
    os << "\t\tComputed norm:        " << utils.sciformat(inf_norm) 
       << std::endl
       << "\t\tRelative Tolerance:   " << utils.sciformat(rtol) 
       << std::endl
       << "\t\tAbsolute Tolerance:   " << utils.sciformat(rtol) 
       << std::endl;
  }

  if (passed)
    return 0;
  else
    return 1;
}

LOCA::MultiContinuation::ExtendedMultiVector::ExtendedMultiVector(
		     const Teuchos::RCP<LOCA::GlobalData>& global_data,
		     const NOX::Abstract::MultiVector& xVec,
		     const NOX::Abstract::MultiVector::DenseMatrix& params) :
  LOCA::Extended::MultiVector(global_data, xVec.numVectors(), 1, 
			      params.numRows())
{
  LOCA::Extended::MultiVector::setMultiVectorPtr(0, xVec.clone(NOX::DeepCopy));
  LOCA::Extended::MultiVector::getScalars()->assign(params);
}

void
LOCA::BorderedSolver::HouseholderQR::computeHouseholderVector(
              int col,
              const NOX::Abstract::MultiVector::DenseMatrix& A1,
              const NOX::Abstract::MultiVector& A2,
              NOX::Abstract::MultiVector::DenseMatrix& V1,
              NOX::Abstract::MultiVector& V2,
              double& beta)
{
  double houseP = A1(col,col);

  V1(0,0) = 1.0;
  V2[0] = A2[col];

  double sigma = A2[col].innerProduct(A2[col]);
  for (int i=col+1; i<A1.numRows(); i++)
    sigma += A1(i,col)*A1(i,col);

  if (sigma == 0.0)
    beta = 0.0;
  else {
    double mu = sqrt(houseP*houseP + sigma);
    if (houseP <= 0.0)
      houseP = houseP - mu;
    else
      houseP = -sigma / (houseP + mu);
    beta = 2.0*houseP*houseP/(sigma + houseP*houseP);

    V2.scale(1.0/houseP);
    for (int i=1; i<V1.numRows(); i++)
      V1(i,0) = A1(i+col,col) / houseP;
  }


  return;
}

NOX::Abstract::Group::ReturnType
LOCA::MultiContinuation::ConstraintInterfaceMVDX::multiplyDX(
               double alpha,
               const NOX::Abstract::MultiVector& input_x,
               NOX::Abstract::MultiVector::DenseMatrix& result_p) const
{


  if (!isDXZero()) {
    const NOX::Abstract::MultiVector* dgdx = getDX();
    input_x.multiply(alpha, *dgdx, result_p);
  }
  else
    result_p.putScalar(0.0);

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

// =============================================================================
// Compute result_p = alpha * dg/dx * input_x.
NOX::Abstract::Group::ReturnType
Ginla::FDM::Constraint::MinDist::
multiplyDX ( double                                    alpha,
             const NOX::Abstract::MultiVector        & input_x,
             NOX::Abstract::MultiVector::DenseMatrix & result_p
           ) const
{ 
  TEUCHOS_ASSERT( komplex_.is_valid_ptr() && !komplex_.is_null() );
  
  TEUCHOS_ASSERT_EQUALITY( result_p.numCols(), input_x.numVectors() );
  
  for ( int k=0; k<input_x.numVectors(); k++ )
  {
      const Epetra_Vector & xE =
              Teuchos::dyn_cast<const NOX::Epetra::Vector>( input_x[0] ).getEpetraVector();
      Teuchos::RCP<ComplexVector> xPsi = komplex_->real2complex( xE );

      result_p(0,k) = alpha * std::imag( psiRef_->dot(*xPsi) );
  }
  
  return NOX::Abstract::Group::Ok;
}

void
LOCA::BorderedSolver::HouseholderQR::applyCompactWY(
             const NOX::Abstract::MultiVector::DenseMatrix& Y1,
             const NOX::Abstract::MultiVector& Y2,
             const NOX::Abstract::MultiVector::DenseMatrix& T,
             const NOX::Abstract::MultiVector::DenseMatrix* input1,
             const NOX::Abstract::MultiVector* input2,
             NOX::Abstract::MultiVector::DenseMatrix& result1,
             NOX::Abstract::MultiVector& result2,
             bool useTranspose) const
{
  bool isZeroX1 = (input1 == NULL);
  bool isZeroX2 = (input2 == NULL);

  if (!isZeroX1)
    result1.assign(*input1);
  if (!isZeroX2)
    result2 = *input2;

  applyCompactWY(Y1, Y2, T, result1, result2, isZeroX1, isZeroX2,
         useTranspose);
}

void
LOCA::BorderedSolver::HouseholderQR::applyHouseholderVector(
               const NOX::Abstract::MultiVector::DenseMatrix& V1,
               const NOX::Abstract::MultiVector& V2,
               double beta,
               NOX::Abstract::MultiVector::DenseMatrix& A1,
               NOX::Abstract::MultiVector& A2)
{
  int nColsA = A2.numVectors();

  // Compute u = V2^T * A2
  NOX::Abstract::MultiVector::DenseMatrix u(1, nColsA);
  A2.multiply(1.0, V2, u);

  // Compute u = u + V1^T * A_P
  u.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, V1, A1, 1.0);

  // Compute A1 = A1 - b*V1*u
  A1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -beta, V1, u, 1.0);

  // Compute A2 = A2 - b*V2*u
  A2.update(Teuchos::NO_TRANS, -beta, V2, u, 1.0);
}

void
NOX::Thyra::MultiVector::
multiply(double alpha,
     const NOX::Abstract::MultiVector& y,
     NOX::Abstract::MultiVector::DenseMatrix& b) const
{
  const NOX::Thyra::MultiVector& yy =
    dynamic_cast<const NOX::Thyra::MultiVector&>(y);

  int m = b.numRows();
  int n = b.numCols();
  Teuchos::RCP< ::Thyra::MultiVectorBase<double> > bb =
    ::Thyra::createMembersView(
      yy.thyraMultiVec->domain(),
      RTOpPack::SubMultiVectorView<double>(0, m, 0, n,
        Teuchos::arcp(b.values(), 0, b.stride()*b.numCols(), false),
        b.stride()
        )
      );

  ::Thyra::apply(*yy.thyraMultiVec, ::Thyra::CONJTRANS, *thyraMultiVec,
                   bb.ptr(), alpha, 0.0);
}

void
LOCA::MultiContinuation::ConstrainedGroup::fillC(
	                     NOX::Abstract::MultiVector::DenseMatrix& C) const
{
  std::string callingFunction = 
    "LOCA::MultiContinuation::ConstrainedGroup::fillC";

  Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> my_C = 
    dfdpMultiVec->getScalars();

  // If the underlying system isn't bordered, we're done
  if (!isBordered) {
    C.assign(*my_C);
    return;
  }

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

  Teuchos::RCP<const NOX::Abstract::MultiVector> my_A = 
    dfdpMultiVec->getXMultiVec();
  
  // Create views for underlying group
  int w = bordered_grp->getBorderedWidth();
  NOX::Abstract::MultiVector::DenseMatrix underlyingC(Teuchos::View, C,
						      w, w, 0, 0);

  // Combine blocks in underlying group
  bordered_grp->fillC(underlyingC);

  // Create views for my blocks
  NOX::Abstract::MultiVector::DenseMatrix my_A_p(Teuchos::View, C,
						 w, numParams, 0, w);
  NOX::Abstract::MultiVector::DenseMatrix my_B_p(Teuchos::View, C,
						 numParams, w, w, 0);
  NOX::Abstract::MultiVector::DenseMatrix my_CC(Teuchos::View, C,
						numParams, numParams, w, w);

  // Extract solution component from my_A and store in my_A_p
  bordered_grp->extractParameterComponent(false, *my_A, my_A_p);

  // Extract solution component from my_B and store in my_B_p
  if (isZeroB)
    my_B_p.putScalar(0.0);
  else
    bordered_grp->extractParameterComponent(true, *my_B, my_B_p);

  // Copy in my_C
  my_CC.assign(*my_C);
}

NOX::Abstract::Group::ReturnType 
LOCA::BorderedSolver::LowerTriangularBlockElimination::
solve(Teuchos::ParameterList& params,
      const LOCA::BorderedSolver::AbstractOperator& op,
      const LOCA::MultiContinuation::ConstraintInterface& B,
      const NOX::Abstract::MultiVector::DenseMatrix& C,
      const NOX::Abstract::MultiVector* F,
      const NOX::Abstract::MultiVector::DenseMatrix* G,
      NOX::Abstract::MultiVector& X,
      NOX::Abstract::MultiVector::DenseMatrix& Y) const
{
  string callingFunction = 
    "LOCA::BorderedSolver::LowerTriangularBlockElimination::solve()";
  NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
  NOX::Abstract::Group::ReturnType status;

  // Determine if X or Y is zero
  bool isZeroF = (F == NULL);
  bool isZeroG = (G == NULL);
  bool isZeroB = B.isDXZero();
  bool isZeroX = isZeroF;
  bool isZeroY = isZeroG && (isZeroB  || isZeroX);

  // First compute X
  if (isZeroX)
    X.init(0.0);
  else {
    // Solve X = J^-1 F, note F must be nonzero
    status = op.applyInverse(params, *F, X);
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Now compute Y
  if (isZeroY)
    Y.putScalar(0.0);
  else {
    // Compute G - B^T*X and store in Y
    if (isZeroG) 
      B.multiplyDX(-1.0, X, Y);
    else {
      Y.assign(*G);
      if (!isZeroB && !isZeroX) {
	NOX::Abstract::MultiVector::DenseMatrix T(Y.numRows(),Y.numCols());
	B.multiplyDX(1.0, X, T);
	Y -= T;
      }
    }

    // Overwrite Y with Y = C^-1 * (G - B^T*X)
    NOX::Abstract::MultiVector::DenseMatrix M(C);
    int *ipiv = new int[M.numRows()];
    Teuchos::LAPACK<int,double> L;
    int info;
    L.GETRF(M.numRows(), M.numCols(), M.values(), M.stride(), ipiv, &info);
    if (info != 0) {
      status = NOX::Abstract::Group::Failed;
      finalStatus = 
	globalData->locaErrorCheck->combineAndCheckReturnTypes(
							      status, 
							      finalStatus,
							      callingFunction);
    }
    L.GETRS('N', M.numRows(), Y.numCols(), M.values(), M.stride(), ipiv, 
	    Y.values(), Y.stride(), &info);
    delete [] ipiv;
    if (info != 0) {
      status = NOX::Abstract::Group::Failed;
      finalStatus = 
	globalData->locaErrorCheck->combineAndCheckReturnTypes(
							     status, 
							     finalStatus,
							     callingFunction);
    }
  }

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