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

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