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

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