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

void NOX::MeritFunction::SumOfSquares::
computeGradient(const NOX::Abstract::Group& grp,
        NOX::Abstract::Vector& result) const
{
  if ( !(grp.isF()) ) {
    utils->err()
      << "ERROR: NOX::MeritFunction::SumOfSquares::computeGradient() - "
      << "F has not been computed yet!.  Please call "
      << "computeF() on the group passed into this function."
      << std::endl;
    throw "NOX Error";
  }

  if ( !(grp.isJacobian()) ) {
    utils->err()
      << "ERROR: NOX::MeritFunction::SumOfSquares::computeGradient() - "
      << "Jacobian has not been computed yet!.  Please call "
      << "computeJacobian() on the group passed into this function."
      << std::endl;
    throw "NOX Error";
  }

  NOX::Abstract::Group::ReturnType status =
    grp.applyJacobianTranspose(grp.getF(), result);

  if (status != NOX::Abstract::Group::Ok) {
    utils->err() << "ERROR: NOX::MeritFunction::SumOfSquares::compute"
         << "Gradient - applyJacobianTranspose failed!" << std::endl;
    throw "NOX Error";
  }

  return;
}

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

  // If the Jacobian is not computed, approximate it with
  // directional derivatives.  dir^T J^T F = F^T Jd
  if (!(grp.isJacobian()))
    return this->computeSlopeWithoutJacobian(dir, grp);
  // If the Jacobian is computed but doesn't support a gradient,
  // employ a different form for the inner product, eg
  // return <v, F> = F' * J * dir = <J'F, dir> = <g, dir>
  else if(!(grp.isGradient()))
    return this->computeSlopeWithoutJacobianTranspose(dir, grp);

  this->computeGradient(grp, *(tmpVecPtr.get()));

  return dir.innerProduct(*(tmpVecPtr.get()));
}

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

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