本文整理汇总了C++中nox::abstract::Group类的典型用法代码示例。如果您正苦于以下问题:C++ Group类的具体用法?C++ Group怎么用?C++ Group使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Group类的19个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1:
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;
}
开发者ID:00liujj,项目名称:trilinos,代码行数:33,代码来源:NOX_MeritFunction_SumOfSquares.C
示例2:
double
NOX::Solver::TensorBased::getNormModelResidual(
const NOX::Abstract::Vector& dir,
const NOX::Abstract::Group& soln,
bool isTensorModel) const
{
// Compute residual of Newton model...
Teuchos::RCP<NOX::Abstract::Vector> residualPtr =
soln.getF().clone(ShapeCopy);
soln.applyJacobian(dir, *residualPtr);
numJvMults++;
residualPtr->update(1.0, soln.getF(), 1.0);
// Compute residual of Tensor model, if requested...
if (isTensorModel)
{
double tmp = sVecPtr->innerProduct(dir);
if (utilsPtr->isPrintType(NOX::Utils::Details))
utilsPtr->out() << " sc'*dt = " << utilsPtr->sciformat(tmp, 6) << std::endl;
residualPtr->update(tmp*tmp, *aVecPtr, 1.0);
}
double modelNorm = residualPtr->norm();
return modelNorm;
}
开发者ID:KineticTheory,项目名称:Trilinos,代码行数:26,代码来源:NOX_Solver_TensorBased.C
示例3: switch
double NOX::StatusTest::NormF::computeNorm(const NOX::Abstract::Group& grp)
{
if (!grp.isF())
return -1.0;
double norm;
int n = grp.getX().length();
switch (normType)
{
case NOX::Abstract::Vector::TwoNorm:
norm = grp.getNormF();
if (scaleType == Scaled)
norm /= sqrt(1.0 * n);
break;
default:
norm = grp.getF().norm(normType);
if (scaleType == Scaled)
norm /= n;
break;
}
return norm;
}
开发者ID:,项目名称:,代码行数:27,代码来源:
示例4: if
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;
}
开发者ID:gitter-badger,项目名称:quinoa,代码行数:57,代码来源:NOX_Direction_ModifiedNewton.C
示例5: if
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;
}
开发者ID:haripandey,项目名称:trilinos,代码行数:51,代码来源:NOX_Direction_Newton.C
示例6: return
double NOX::MeritFunction::SumOfSquares::
computef(const NOX::Abstract::Group& grp) const
{
if ( !(grp.isF()) ) {
utils->err()
<< "ERROR: NOX::MeritFunction::SumOfSquares::computef() - "
<< "F has not been computed yet!. Please call "
<< "computeF() on the group passed into this function."
<< std::endl;
throw "NOX Error";
}
return (0.5 * grp.getNormF() * grp.getNormF());
}
开发者ID:00liujj,项目名称:trilinos,代码行数:14,代码来源:NOX_MeritFunction_SumOfSquares.C
示例7:
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;
}
开发者ID:00liujj,项目名称:trilinos,代码行数:14,代码来源:NOX_LineSearch_Polynomial.C
示例8:
bool NOX::Direction::ModifiedNewton::rescueBadNewtonSolve(const NOX::Abstract::Group& grp) const
{
//! Check if the "rescue" option has been selected
if (!doRescue)
return false;
//! See if the group has compute the accuracy
double accuracy;
NOX::Abstract::Group::ReturnType status = oldJacobianGrpPtr->getNormLastLinearSolveResidual(accuracy);
// If this functionality is not supported in the group, return false
/* NOTE FROM TAMMY: We could later modify this to acutally caluclate
the error itself if it's just a matter of the status being
NotDefined. */
if (status != NOX::Abstract::Group::Ok)
return false;
// Check if there is any improvement in the relative residual
double normF = grp.getNormF();
// If we can't reduce the relative norm at all, we're not happy
if (accuracy >= normF)
return false;
// Otherwise, we just print a warning and keep going
if (utils->isPrintType(NOX::Utils::Warning))
utils->out() << "WARNING: NOX::Direction::ModifiedNewton::compute - Unable to achieve desired linear solve accuracy." << std::endl;
return true;
}
开发者ID:gitter-badger,项目名称:quinoa,代码行数:30,代码来源:NOX_Direction_ModifiedNewton.C
示例9: getNormModelResidual
void
NOX::Solver::TensorBased::printDirectionInfo(std::string dirName,
const NOX::Abstract::Vector& dir,
const NOX::Abstract::Group& soln,
bool isTensorModel) const
{
double dirNorm = dir.norm();
double residual = getNormModelResidual(dir, soln, isTensorModel);
double residualRel = residual / soln.getNormF();
double fprime = getDirectionalDerivative(dir, soln);
double fprimeRel = fprime / dirNorm;
if (utilsPtr->isPrintType(NOX::Utils::Details))
{
utilsPtr->out() << " " << dirName << " norm of model residual = "
<< utilsPtr->sciformat(residual, 6) << " (abs) "
<< utilsPtr->sciformat(residualRel, 6) << " (rel)" << std::endl;
utilsPtr->out() << " " << dirName << " directional derivative = "
<< utilsPtr->sciformat(fprime, 6) << " (abs) "
<< utilsPtr->sciformat(fprimeRel, 6) << " (rel)" << std::endl;
utilsPtr->out() << " " << dirName << " norm = "
<< utilsPtr->sciformat(dirNorm, 6) << std::endl;
}
}
开发者ID:KineticTheory,项目名称:Trilinos,代码行数:26,代码来源:NOX_Solver_TensorBased.C
示例10: if
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()));
}
开发者ID:00liujj,项目名称:trilinos,代码行数:21,代码来源:NOX_MeritFunction_SumOfSquares.C
示例11: throwError
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;
}
开发者ID:,项目名称:,代码行数:59,代码来源:
示例12: computeNorm
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;
}
开发者ID:,项目名称:,代码行数:14,代码来源:
示例13: if
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;
}
开发者ID:KineticTheory,项目名称:Trilinos,代码行数:53,代码来源:NOX_Solver_TensorBased.C
示例14:
bool NOX::Direction::Broyden::doRestart(NOX::Abstract::Group& soln,
const NOX::Solver::LineSearchBased& solver)
{
// Test 1 - First iteration!
if (solver.getNumIterations() == 0)
return true;
// Test 2 - Frequency
if (cnt >= cntMax)
return true;
// Test 3 - Last step was zero!
if (solver.getStepSize() == 0.0)
return true;
// Test 4 - Check for convergence rate
convRate = soln.getNormF() / solver.getPreviousSolutionGroup().getNormF();
if (convRate > maxConvRate)
return true;
return false;
}
开发者ID:haripandey,项目名称:trilinos,代码行数:22,代码来源:NOX_Direction_Broyden.C
示例15: while
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);
}
开发者ID:00liujj,项目名称:trilinos,代码行数:83,代码来源:NOX_LineSearch_Backtrack.C
示例16: if
// **************************************************************************
// *** 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;
}
//.........这里部分代码省略.........
开发者ID:gitter-badger,项目名称:quinoa,代码行数:101,代码来源:NOX_Direction_Utils_InexactNewton.C
示例17: setGroupForComputeF
void MatrixFree::setGroupForComputeF(const NOX::Abstract::Group& group)
{
useGroupForComputeF = true;
groupPtr = group.clone();
return;
}
开发者ID:00liujj,项目名称:trilinos,代码行数:6,代码来源:NOX_Epetra_MatrixFree.C
示例18: if
NOX::Abstract::Group::ReturnType
LOCA::Eigensolver::DGGEVStrategy::computeEigenvalues(
NOX::Abstract::Group& group,
Teuchos::RCP< std::vector<double> >& evals_r,
Teuchos::RCP< std::vector<double> >& evals_i,
Teuchos::RCP< NOX::Abstract::MultiVector >& evecs_r,
Teuchos::RCP< NOX::Abstract::MultiVector >& evecs_i)
{
// Get LAPACK group
LOCA::LAPACK::Group* grp =
dynamic_cast<LOCA::LAPACK::Group*>(&group);
bool hasMassMatrix = true;
if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration)) {
globalData->locaUtils->out()
<< std::endl << globalData->locaUtils->fill(64,'=') << std::endl
<< "LAPACK ";
if (hasMassMatrix)
globalData->locaUtils->out() << "DGGEV ";
else
globalData->locaUtils->out() << "DGEEV ";
globalData->locaUtils->out() << "Eigensolver starting."
<< std::endl << std::endl;;
}
// Make sure Jacobian & mass matrices are fresh
grp->computeJacobian();
if (hasMassMatrix)
grp->computeShiftedMatrix(0.0, 1.0);
// Get data out of group
NOX::LAPACK::Matrix<double>& jacobianMatrix = grp->getJacobianMatrix();
NOX::LAPACK::Matrix<double>& massMatrix = grp->getShiftedMatrix();
// Size of matrix
int n = jacobianMatrix.numRows();
int lda = jacobianMatrix.numRows();
int ldb = massMatrix.numRows();
// Space to hold right eigenvectors
double *vr = new double[n*n];
// Space to hold real and imaginary eigenvalues
double *alphar = new double[n];
double *alphai = new double[n];
double *beta = new double[n];
// Size of work array, set to -1 to do a workspace query
int lwork = -1;
// Initial work "array"
double work0;
// Actual work array
double *work;
// Return code
int info;
// Copy Jacobian matrix since lapack routines overwrite it
NOX::LAPACK::Matrix<double> J(jacobianMatrix);
NOX::LAPACK::Matrix<double> M;
// First do a workspace query
if (hasMassMatrix) {
// Copy mass matrix since lapack routines overwrite it
M = massMatrix;
DGGEV_F77("N", "V", &n, &J(0,0), &lda, &M(0,0), &ldb, alphar, alphai, beta,
vr, &n, vr, &n, &work0, &lwork, &info);
}
else {
DGEEV_F77("N", "V", &n, &J(0,0), &lda, alphar, alphai,
vr, &n, vr, &n, &work0, &lwork, &info);
}
// Allocate work array
lwork = (int) work0;
work = new double[lwork];
// Calculate eigenvalues, eigenvectors
if (hasMassMatrix) {
DGGEV_F77("N", "V", &n, &J(0,0), &lda, &M(0,0), &ldb, alphar, alphai, beta,
vr, &n, vr, &n, work, &lwork, &info);
}
else {
DGEEV_F77("N", "V", &n, &J(0,0), &lda, alphar, alphai,
vr, &n, vr, &n, work, &lwork, &info);
}
// Check for success
if (info != 0)
return NOX::Abstract::Group::Failed;
// Compute all of the eigenvalues and eigenvectors before sorting
std::vector<double> evals_r_tmp(n);
//.........这里部分代码省略.........
开发者ID:00liujj,项目名称:trilinos,代码行数:101,代码来源:LOCA_Eigensolver_DGGEVStrategy.C
示例19: throwError
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>
//.........这里部分代码省略.........
开发者ID:haripandey,项目名称:trilinos,代码行数:101,代码来源:NOX_Direction_Broyden.C
注:本文中的nox::abstract::Group类示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论