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

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

bool NOX::Direction::ModifiedNewton::
compute(NOX::Abstract::Vector& dir, 
	NOX::Abstract::Group& soln, 
	const NOX::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");

  maxAgeOfJacobian = paramsPtr->sublist("Modified-Newton").get("Max Age of Jacobian", 10);

  if (Teuchos::is_null(oldJacobianGrpPtr)) {
    oldJacobianGrpPtr = soln.clone(DeepCopy);
  }
  NOX::Abstract::Group& oldJacobianGrp = *oldJacobianGrpPtr;
  
  status = NOX::Abstract::Group::Failed;
  while (status != NOX::Abstract::Group::Ok) {
    // Conditionally compute Jacobian at current solution.
    if ( (ageOfJacobian == -1) || (ageOfJacobian == maxAgeOfJacobian) ) {

      if (ageOfJacobian > 0) 
        oldJacobianGrp = soln;
      status = oldJacobianGrp.computeJacobian();
      if (status != NOX::Abstract::Group::Ok) 
        throwError("compute", "Unable to compute Jacobian");
      ageOfJacobian = 1;
    } 
    else 
      ageOfJacobian++;

    // Compute the Modified Newton direction
    status = oldJacobianGrp.applyJacobianInverse(paramsPtr->sublist("Modified-Newton").sublist("Linear Solver"), soln.getF(), dir);
    dir.scale(-1.0);

    // It didn't converge, but maybe we can recover.
    if ((status != NOX::Abstract::Group::Ok) &&
        (doRescue == false)) {
      throwError("compute", "Unable to solve Newton system");
    }
    else if ((status != NOX::Abstract::Group::Ok) &&
             (doRescue == true)) {
      if (utils->isPrintType(NOX::Utils::Warning))
        utils->out() << "WARNING: NOX::Direction::ModifiedNewton::compute() - "
             << "Linear solve failed to achieve convergence - "
             << "using the step anyway since \"Rescue Bad Newton Solve\" "
             << "is true. Also, flagging recompute of Jacobian." << std::endl;
      ageOfJacobian = maxAgeOfJacobian;
      status = NOX::Abstract::Group::Ok;
    }
  }

  return true;
}

bool
NOX::Solver::TensorBased::implementGlobalStrategy(NOX::Abstract::Group& newGrp,
                      double& in_stepSize,
                      const NOX::Solver::Generic& s)
{
  bool ok;
  counter.incrementNumLineSearches();
  isNewtonDirection = false;
  NOX::Abstract::Vector& searchDirection = *tensorVecPtr;

  if ((counter.getNumLineSearches() == 1)  ||  (lsType == Newton))
  {
    isNewtonDirection = true;
    searchDirection = *newtonVecPtr;
  }

  // Do line search and compute new soln.
  if ((lsType != Dual) || (isNewtonDirection))
    ok = performLinesearch(newGrp, in_stepSize, searchDirection, s);
  else if (lsType == Dual)
  {
    double fTensor = 0.0;
    double fNew = 0.0;
    double tensorStep = 1.0;
    bool isTensorDescent = false;

    const Abstract::Group& oldGrp = s.getPreviousSolutionGroup();
    double fprime = slopeObj.computeSlope(searchDirection, oldGrp);

    // Backtrack along tensor direction if it is descent direction.
    if (fprime < 0)
    {
      ok = performLinesearch(newGrp, in_stepSize, searchDirection, s);
      assert(ok);
      fTensor = 0.5 * newGrp.getNormF() * newGrp.getNormF();
      tensorStep = in_stepSize;
      isTensorDescent = true;
    }

    // Backtrack along the Newton direction.
    ok = performLinesearch(newGrp, in_stepSize, *newtonVecPtr, s);
    fNew = 0.5 * newGrp.getNormF() * newGrp.getNormF();

    // If backtracking on the tensor step produced a better step, then use it.
    if (isTensorDescent  &&  (fTensor <= fNew))
    {
      newGrp.computeX(oldGrp, *tensorVecPtr, tensorStep);
      newGrp.computeF();
    }
  }

  return ok;
}

bool NOX::Direction::Newton::compute(NOX::Abstract::Vector& dir, 
				     NOX::Abstract::Group& soln, 
				     const NOX::Solver::Generic& solver)
{
  NOX::Abstract::Group::ReturnType status;

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

  // Reset the linear solver tolerance.
  if (useAdjustableForcingTerm) {
    resetForcingTerm(soln, solver.getPreviousSolutionGroup(), 
		     solver.getNumIterations(), solver);
  }
  else { 
    if (utils->isPrintType(Utils::Details)) {
      utils->out() << "       CALCULATING FORCING TERM" << endl;
      utils->out() << "       Method: Constant" << endl;
      utils->out() << "       Forcing Term: " << eta_k << endl;
    }
  }

  // Compute Jacobian at current solution.
  status = soln.computeJacobian();
  if (status != NOX::Abstract::Group::Ok) 
    NOX::Direction::Newton::throwError("compute", "Unable to compute Jacobian");
  
  // Compute the Newton direction
  status = soln.computeNewton(paramsPtr->sublist("Newton").sublist("Linear Solver"));
  
  // It didn't converge, but maybe we can recover. 
  if ((status != NOX::Abstract::Group::Ok) &&
      (doRescue == false)) {
    NOX::Direction::Newton::throwError("compute", 
				       "Unable to solve Newton system");
  }
  else if ((status != NOX::Abstract::Group::Ok) &&
	   (doRescue == true)) {
    if (utils->isPrintType(NOX::Utils::Warning))
      utils->out() << "WARNING: NOX::Direction::Newton::compute() - Linear solve "
	   << "failed to achieve convergence - using the step anyway " 
	   << "since \"Rescue Bad Newton Solve\" is true " << endl;
  }

  // Set search direction.
  dir = soln.getNewton();

  return true;
}

bool NOX::LineSearch::Polynomial::
updateGrp(NOX::Abstract::Group& newGrp,
      const NOX::Abstract::Group& oldGrp,
      const NOX::Abstract::Vector& dir,
      double step) const
{
  newGrp.computeX(oldGrp, dir, step);

  NOX::Abstract::Group::ReturnType status = newGrp.computeF();
  if (status != NOX::Abstract::Group::Ok)
    return false;

  return true;
}

void NOX::StatusTest::NormF::relativeSetup(NOX::Abstract::Group& initialGuess)
{
  NOX::Abstract::Group::ReturnType rtype;
  rtype = initialGuess.computeF();
  if (rtype != NOX::Abstract::Group::Ok) 
  {
    utils.err() << "NOX::StatusTest::NormF::NormF - Unable to compute F" 
		<< endl;
    throw "NOX Error";
  }
    
  initialTolerance = computeNorm(initialGuess); 
  trueTolerance = specifiedTolerance * initialTolerance;
}

bool NOX::LineSearch::Backtrack::
compute(NOX::Abstract::Group& grp, double& step,
    const NOX::Abstract::Vector& dir,
    const NOX::Solver::Generic& s)
{
  const Abstract::Group& oldGrp = s.getPreviousSolutionGroup();
  double oldF = meritFunctionPtr->computef(oldGrp);
  double newF;
  bool isFailed = false;

  step = defaultStep;
  grp.computeX(oldGrp, dir, step);

  NOX::Abstract::Group::ReturnType rtype;

  rtype = grp.computeF();
  if (rtype != NOX::Abstract::Group::Ok)
  {
    utils->err() << "NOX::LineSearch::BackTrack::compute - Unable to compute F"
        << std::endl;
    throw "NOX Error";
  }

  newF = meritFunctionPtr->computef(grp);
  int nIters = 1;

  if (utils->isPrintType(Utils::InnerIteration))
  {
   utils->out() << "\n" << Utils::fill(72) << "\n"
           << "-- Backtrack Line Search -- \n";
  }

  NOX::StatusTest::FiniteValue checkNAN;

  while ( ((newF >= oldF) || (checkNAN.finiteNumberTest(newF) !=0))
     && (!isFailed))
  {

    if (utils->isPrintType(Utils::InnerIteration))
    {
      utils->out() << std::setw(3) << nIters << ":";
      utils->out() << " step = " << utils->sciformat(step);
      utils->out() << " old f = " << utils->sciformat(oldF);
      utils->out() << " new f = " << utils->sciformat(newF);
      utils->out() << std::endl;
    }

    nIters ++;
    step = step * reductionFactor;

    if ((step < minStep) || (nIters > maxIters))
    {
      isFailed = true;
      step = recoveryStep;
    }

    grp.computeX(oldGrp, dir, step);

    rtype = grp.computeF();
    if (rtype != NOX::Abstract::Group::Ok)
    {
      utils->err() << "NOX::LineSearch::BackTrack::compute - Unable to compute F" << std::endl;
      throw "NOX Error";
    }

    newF = meritFunctionPtr->computef(grp);
  }

  if (utils->isPrintType(Utils::InnerIteration))
  {
    utils->out() << std::setw(3) << nIters << ":";
    utils->out() << " step = " << utils->sciformat(step);
    utils->out() << " old f = " << utils->sciformat(oldF);
    utils->out() << " new f = " << utils->sciformat(newF);
    if (isFailed)
      utils->out() << " (USING RECOVERY STEP!)" << std::endl;
    else
      utils->out() << " (STEP ACCEPTED!)" << std::endl;
    utils->out() << Utils::fill(72) << "\n" << std::endl;
  }

  return (!isFailed);
}

bool NOX::Direction::Broyden::compute(NOX::Abstract::Vector& dir, 
				      NOX::Abstract::Group& soln, 
				      const NOX::Solver::LineSearchBased& solver)
{
  // Return value for group operations (temp variable)
  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");

  // Check for restart
  if (doRestart(soln, solver))
  {
    // Reset memory
    memory.reset();

    // Update group
    if (Teuchos::is_null(oldJacobianGrpPtr))
      oldJacobianGrpPtr = soln.clone(NOX::DeepCopy);
    else
      // RPP - update the entire group (this grabs state vectors in xyce).
      // Otherwise, xyce is forced to recalculate F at each iteration.
      //oldJacobianGrpPtr->setX(soln.getX());
      *oldJacobianGrpPtr = soln;

    // Calcuate new Jacobian
    if (utils->isPrintType(NOX::Utils::Details))
      utils->out() << "       Recomputing Jacobian" << endl;
 
    status = oldJacobianGrpPtr->computeJacobian();
    if (status != NOX::Abstract::Group::Ok) 
      throwError("compute", "Unable to compute Jacobian");

    // Reset counter
    cnt = 0;
  }

  // If necesary, scale the s-vector from the last iteration
  if (!memory.empty()) 
  {
    double step = solver.getStepSize();
    memory[memory.size() - 1].setStep(step);
  }

  // --- Calculate the Broyden direction ---

  // Compute inexact forcing term if requested.
  inexactNewtonUtils.computeForcingTerm(soln, 
					solver.getPreviousSolutionGroup(),
					solver.getNumIterations(),
					solver);

  // dir = - J_old^{-1} * F
  cnt ++;
  status = oldJacobianGrpPtr->applyJacobianInverse(*lsParamsPtr, 
						   soln.getF(), 
						   dir);
  if (status != NOX::Abstract::Group::Ok) 
    throwError("compute", "Unable to apply Jacobian inverse");
  dir.scale(-1.0);

  // Apply the Broyden modifications to the old Jacobian (implicitly)
  if (!memory.empty()) 
  {
    // Number of elements in the memory
    int m = memory.size();

    // Information corresponding to index i
    double step;
    Teuchos::RCP<const NOX::Abstract::Vector> sPtr;

    // Information corresponding to index i + 1 
    // (initialized for i = -1)
    double stepNext = memory[0].step();
    Teuchos::RCP<const NOX::Abstract::Vector> sPtrNext = 
      memory[0].sPtr();

    // Intermediate storage
    double a, b, c, denom;

    for (int i = 0; i < m-1; i ++)
    {
      step = stepNext;
      sPtr = sPtrNext;
      stepNext = memory[i+1].step();
      sPtrNext = memory[i+1].sPtr();

      a = step / stepNext;
      b = step - 1;
      c = sPtr->innerProduct(dir) / memory[i].sNormSqr();

      dir.update(a * c, *sPtrNext, b * c, *sPtr, 1.0);
    }

    step = stepNext;
    sPtr = sPtrNext;

    a = sPtr->innerProduct(dir);		// <s,z>
//.........这里部分代码省略.........