C++ Group::getX方法代码示例

double NOX::StatusTest::NormF::computeNorm(const NOX::Abstract::Group& grp)
{
  if (!grp.isF())
    return -1.0;

  double norm;
  int n = grp.getX().length();

  switch (normType) 
  {
    
  case NOX::Abstract::Vector::TwoNorm:
    norm = grp.getNormF();
    if (scaleType == Scaled)
      norm /= sqrt(1.0 * n);
    break;

  default:
    norm = grp.getF().norm(normType);
    if (scaleType == Scaled)
      norm /= n;
    break;

  }

  return norm;
}

bool NOX::Direction::QuasiNewton::compute(NOX::Abstract::Vector& dir, 
					  NOX::Abstract::Group& soln, 
					  const Solver::Generic& solver)
{
  NOX::Abstract::Group::ReturnType status;
  
  // Compute F at current solution
  status = soln.computeF();
  if (status != NOX::Abstract::Group::Ok) 
    throwError("compute", "Unable to compute F");

  // Compute Jacobian at current solution.
  status = soln.computeJacobian();
  if (status != NOX::Abstract::Group::Ok) 
    throwError("compute", "Unable to compute Jacobian");

  // Compute the gradient at the current solution
  status = soln.computeGradient();
  if (status != NOX::Abstract::Group::Ok) 
    throwError("compute", "Unable to compute gradient");

  // Push the old information onto the memory, but only after at least one previous iteration
  if (solver.getNumIterations() > 0) 
  {
    const NOX::Abstract::Group& oldSoln = solver.getPreviousSolutionGroup();
    if (oldSoln.isGradient())
      memory.add(soln.getX(), oldSoln.getX(), soln.getGradient(), oldSoln.getGradient());
  }

  // *** Calculate the QN direction ***
  
  // d = -g
  dir = soln.getGradient();
  dir.scale(-1.0);

  if (!memory.empty()) 
  {

    int m = memory.size();
    vector<double> alpha(m);
    double beta;
  
    for (int i = m-1; i >= 0; i --)
    {
      alpha[i] = memory[i].rho() * dir.innerProduct( memory[i].s() );
      dir.update(-1.0 * alpha[i], memory[i].y(), 1.0);
    }

    dir.scale( memory[m-1].sdoty() / memory[m-1].ydoty() );

    for (int i = 0; i < m; i ++)
    {
      beta = memory[i].rho() * dir.innerProduct( memory[i].y() );
      dir.update(alpha[i] - beta, memory[i].s(), 1.0);
    }
  }

  return true;
}

double NOX::MeritFunction::SumOfSquares::
computeSlopeWithoutJacobian(const NOX::Abstract::Vector& dir,
                const NOX::Abstract::Group& grp) const
{
  if (Teuchos::is_null(tmpVecPtr))
    tmpVecPtr = grp.getF().clone(NOX::ShapeCopy);

  if (Teuchos::is_null(tmpGrpPtr))
    tmpGrpPtr = grp.clone(NOX::ShapeCopy);


  // Compute the perturbation parameter
  double lambda = 1.0e-6;
  double denominator = dir.norm();

  // Don't divide by zero
  if (denominator == 0.0)
    denominator = 1.0;

  double eta = lambda * (lambda + grp.getX().norm() / denominator);

  // Don't divide by zero
  if (eta == 0.0)
    eta = 1.0e-6;

  // Perturb the solution vector
  tmpVecPtr->update(eta, dir, 1.0, grp.getX(), 0.0);

  // Compute the new F --> F(x + eta * dir)
  tmpGrpPtr->setX(*(tmpVecPtr.get()));
  tmpGrpPtr->computeF();

  // Compute Js = (F(x + eta * dir) - F(x))/eta
  tmpVecPtr->update(-1.0/eta, grp.getF(), 1.0/eta, tmpGrpPtr->getF(), 0.0);

  return(tmpVecPtr->innerProduct(grp.getF()));
}

// **************************************************************************
// *** computeForcingTerm
// **************************************************************************
double NOX::Direction::Utils::InexactNewton::
computeForcingTerm(const NOX::Abstract::Group& soln,
		   const NOX::Abstract::Group& oldsoln, 
		   int niter,
		   const NOX::Solver::Generic& solver,
		   double eta_last)
{
  const std::string indent = "       ";

  if (forcingTermMethod == Constant) {
    if (printing->isPrintType(NOX::Utils::Details)) {
      printing->out() << indent << "CALCULATING FORCING TERM" << std::endl;
      printing->out() << indent << "Method: Constant" << std::endl;
      printing->out() << indent << "Forcing Term: " << eta_k << std::endl;
    }
    if (setTolerance)
      paramsPtr->sublist(directionMethod).sublist("Linear Solver").
	set("Tolerance", eta_k);

    return eta_k;
  }

  // Get linear solver current tolerance. 
  // NOTE: These values are changing at each nonlinear iteration and 
  // must either be updated from the parameter list each time a compute 
  // is called or supplied during the function call!
  double eta_km1 = 0.0;
  if (eta_last < 0.0)
    eta_km1 = paramsPtr->sublist(directionMethod).
      sublist("Linear Solver").get("Tolerance", 0.0);
  else
    eta_km1 = eta_last;

  // Tolerance may have been adjusted in a line search algorithm so we 
  // have to account for this.
  const NOX::Solver::LineSearchBased* solverPtr = 0;
  solverPtr = dynamic_cast<const NOX::Solver::LineSearchBased*>(&solver);
  if (solverPtr != 0) {
    eta_km1 = 1.0 - solverPtr->getStepSize() * (1.0 - eta_km1);
  }

  if (printing->isPrintType(NOX::Utils::Details)) {
    printing->out() << indent << "CALCULATING FORCING TERM" << std::endl;
    printing->out() << indent << "Method: " << method << std::endl;
  }


  if (forcingTermMethod == Type1) {
    
    if (niter == 0) {
      
      eta_k = eta_initial;

    }
    else {

      // Return norm of predicted F

      // do NOT use the following lines!! This does NOT account for 
      // line search step length taken.
      // const double normpredf = 0.0;
      // oldsoln.getNormLastLinearSolveResidual(normpredf);
      
      // Create a new vector to be the predicted RHS
      if (Teuchos::is_null(predRhs)) {
	predRhs = oldsoln.getF().clone(ShapeCopy);
      }
      if (Teuchos::is_null(stepDir)) {
	stepDir = oldsoln.getF().clone(ShapeCopy);
      }
      
      // stepDir = X - oldX (i.e., the step times the direction)
      stepDir->update(1.0, soln.getX(), -1.0, oldsoln.getX(), 0);
      
      // Compute predRhs = Jacobian * step * dir
      if (!(oldsoln.isJacobian())) {
	if (printing->isPrintType(NOX::Utils::Details)) {
	  printing->out() << "WARNING: NOX::InexactNewtonUtils::resetForcingTerm() - "
	       << "Jacobian is out of date! Recomputing Jacobian." << std::endl;
	}
	const_cast<NOX::Abstract::Group&>(oldsoln).computeJacobian();
      }
      oldsoln.applyJacobian(*stepDir, *predRhs);

      // Compute predRhs = RHSVector + predRhs (this is the predicted RHS)
      predRhs->update(1.0, oldsoln.getF(), 1.0);
      
      // Compute the norms
      double normpredf = predRhs->norm();
      double normf = soln.getNormF();
      double normoldf = oldsoln.getNormF();

      if (printing->isPrintType(NOX::Utils::Details)) {
	printing->out() << indent << "Forcing Term Norm: Using L-2 Norm."
			<< std::endl;
      }

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

//.........这里部分代码省略.........
    // Copy mass matrix since lapack routines overwrite it
    M = massMatrix;

    DGGEV_F77("N", "V", &n, &J(0,0), &lda, &M(0,0), &ldb, alphar, alphai, beta,
          vr, &n, vr, &n, &work0, &lwork, &info);
  }
  else {
    DGEEV_F77("N", "V", &n, &J(0,0), &lda, alphar, alphai,
          vr, &n, vr, &n, &work0, &lwork, &info);
  }

  // Allocate work array
  lwork = (int) work0;
  work = new double[lwork];

  // Calculate eigenvalues, eigenvectors
  if (hasMassMatrix) {
    DGGEV_F77("N", "V", &n, &J(0,0), &lda, &M(0,0), &ldb, alphar, alphai, beta,
          vr, &n, vr, &n, work, &lwork, &info);
  }
  else {
    DGEEV_F77("N", "V", &n, &J(0,0), &lda, alphar, alphai,
          vr, &n, vr, &n, work, &lwork, &info);
  }

  // Check for success
  if (info != 0)
    return NOX::Abstract::Group::Failed;

  // Compute all of the eigenvalues and eigenvectors before sorting
  std::vector<double> evals_r_tmp(n);
  std::vector<double> evals_i_tmp(n);
  Teuchos::RCP<NOX::Abstract::MultiVector> evecs_r_tmp =
    group.getX().createMultiVector(n, NOX::ShapeCopy);
  Teuchos::RCP<NOX::Abstract::MultiVector>evecs_i_tmp =
    group.getX().createMultiVector(n, NOX::ShapeCopy);
  NOX::LAPACK::Vector* tmpr;
  NOX::LAPACK::Vector* tmpi;
  double rnext;
  double inext;
  bool isComplexEval = false;
  bool isPrevComplexEval = false;
  for (int j=0; j<n; j++) {

    // Compute eigenvalues
    if (hasMassMatrix) {
      evals_r_tmp[j] = alphar[j]/beta[j];
      evals_i_tmp[j] = alphai[j]/beta[j];
    }
    else {
      evals_r_tmp[j] = alphar[j];
      evals_i_tmp[j] = alphai[j];
    }

    // Compute next eigenvalue
    if (!isPrevComplexEval && j < n-1) {
      if (hasMassMatrix) {
    rnext = alphar[j+1]/beta[j+1];
    inext = alphai[j+1]/beta[j+1];
      }
      else {
    rnext = alphar[j+1];
    inext = alphai[j+1];
      }

      // Determine if this eigenvalue is a complex conjugate pair

bool
NOX::Solver::TensorBased::computeTensorDirection(NOX::Abstract::Group& soln,
                     const NOX::Solver::Generic& solver)
{
  NOX::Abstract::Group::ReturnType dir_status;

  Teuchos::ParameterList& linearParams = paramsPtr->sublist("Direction").
    sublist(paramsPtr->sublist("Direction").
        get("Method","Tensor")).
    sublist("Linear Solver");

  // Compute F at current solution.
  dir_status = soln.computeF();
  if (dir_status != NOX::Abstract::Group::Ok)
    throwError("computeTensorDirection", "Unable to compute F");

  // Compute Jacobian at current solution.
  dir_status = soln.computeJacobian();
  if (dir_status != NOX::Abstract::Group::Ok)
    throwError("computeTensorDirection", "Unable to compute Jacobian");

  // Begin processing for the tensor step, if necessary.
  double sDotS = 0.0;
  int tempVal1 = 0;
  if ((nIter > 0)  &&  (requestedBaseStep == TensorStep))
  {
    // Compute the tensor term s = x_{k-1} - x_k
    *sVecPtr = soln.getX();
    sVecPtr->update(1.0, solver.getPreviousSolutionGroup().getX(), -1.0);
    double normS = sVecPtr->norm();
    sDotS = normS * normS;

    // Form the tensor term a = (F_{k-1} - F_k - J*s) / (s^T s)^2
    soln.applyJacobian(*sVecPtr, *aVecPtr);
    numJvMults++;
    aVecPtr->update(1.0, solver.getPreviousSolutionGroup().getF(), -1.0);
    aVecPtr->update(-1.0, soln.getF(), 1.0);
    if (sDotS != 0)
      aVecPtr->scale(1.0 / (sDotS * sDotS));

    // Save old Newton step as initial guess to second system
    *tmpVecPtr = *newtonVecPtr;
    tmpVecPtr->scale(-1.0);   // Rewrite to avoid this?

    // Compute residual of linear system using initial guess...
    soln.applyJacobian(*tmpVecPtr, *residualVecPtr);
    numJvMults++;
    residualVecPtr->update(1.0, solver.getPreviousSolutionGroup().getF(),-1.0);
    double residualNorm = residualVecPtr->norm();

#if DEBUG_LEVEL > 0
    double tmpVecNorm = tmpVecPtr->norm();
    double residualNormRel = residualNorm /
      solver.getPreviousSolutionGroup().getNormF();
    if (utilsPtr->isPrintType(NOX::Utils::Details))
    {
      utilsPtr->out() << "  Norm of initial guess: " << utilsPtr->sciformat(tmpVecNorm, 6)
       << std::endl;
      utilsPtr->out() << "  initg norm of model residual =   "
       << utilsPtr->sciformat(residualNorm, 6) << " (abs)     "
       << utilsPtr->sciformat(residualNormRel, 6) << " (rel)" << std::endl;
    }
#endif

    // Save some parameters and use them later...
    double tol = linearParams.get("Tolerance", 1e-4);
    double relativeResidual = residualNorm /
      solver.getPreviousSolutionGroup().getNormF();

    // Decide whether to use initial guess...
    bool isInitialGuessGood = false;
#ifdef USE_INITIAL_GUESS_LOGIC
    if (relativeResidual < 1.0)
    {
      if (utilsPtr->isPrintType(NOX::Utils::Details))
    utilsPtr->out() << "  Initial guess is good..." << std::endl;
      isInitialGuessGood = true;
      // RPP - Brett please make sure the line below is correct.
      *tensorVecPtr = *tmpVecPtr;
      double newTol = tol / relativeResidual;
      if (newTol > 0.99)
    newTol = 0.99;  // force at least one iteration
      linearParams.set("Tolerance",  newTol);
      if (utilsPtr->isPrintType(NOX::Utils::Details))
    utilsPtr->out() << "  Setting tolerance to " << utilsPtr->sciformat(newTol,6) << std::endl;
    }
    else
#endif // USE_INITIAL_GUESS_LOGIC
    {
      //utilsPtr->out() << "  Initial guess is BAD... do not use!\n";
      isInitialGuessGood = false;
      *residualVecPtr = solver.getPreviousSolutionGroup().getF();
    }

    // Compute the term inv(J)*Fp....
    tmpVecPtr->init(0.0);
    dir_status = soln.applyJacobianInverse(linearParams, *residualVecPtr,
                       *tmpVecPtr);

    // If it didn't converge, maybe we can recover.
//.........这里部分代码省略.........