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

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

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

void MatrixFree::setGroupForComputeF(const NOX::Abstract::Group& group)
{
  useGroupForComputeF = true;
  groupPtr = group.clone();
  return;
}

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