00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_TurningPoint_MinimallyAugmented_Constraint.H"
00043 #include "LOCA_TurningPoint_MinimallyAugmented_AbstractGroup.H"
00044 #include "LOCA_BorderedSolver_AbstractStrategy.H"
00045 #include "LOCA_Parameter_SublistParser.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_ErrorCheck.H"
00048 #include "LOCA_Factory.H"
00049 #include "NOX_Utils.H"
00050 #include "Teuchos_ParameterList.hpp"
00051 #include "LOCA_BorderedSolver_JacobianOperator.H"
00052
00053 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00054 Constraint(
00055 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00056 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00057 const Teuchos::RCP<Teuchos::ParameterList>& tpParams,
00058 const Teuchos::RCP<LOCA::TurningPoint::MinimallyAugmented::AbstractGroup>& g,
00059 bool is_symmetric,
00060 const NOX::Abstract::Vector& a,
00061 const NOX::Abstract::Vector* b,
00062 int bif_param) :
00063 globalData(global_data),
00064 parsedParams(topParams),
00065 turningPointParams(tpParams),
00066 grpPtr(g),
00067 a_vector(a.createMultiVector(1, NOX::DeepCopy)),
00068 b_vector(),
00069 w_vector(a.createMultiVector(1, NOX::ShapeCopy)),
00070 v_vector(a.createMultiVector(1, NOX::ShapeCopy)),
00071 Jv_vector(a.createMultiVector(1, NOX::ShapeCopy)),
00072 sigma_x(a.createMultiVector(1, NOX::ShapeCopy)),
00073 constraints(1, 1),
00074 borderedSolver(),
00075 dn(static_cast<double>(a_vector->length())),
00076 sigma_scale(1.0),
00077 isSymmetric(is_symmetric),
00078 isValidConstraints(false),
00079 isValidDX(false),
00080 bifParamID(1),
00081 updateVectorsEveryContinuationStep(true),
00082 updateVectorsEveryIteration(false)
00083 {
00084
00085 borderedSolver =
00086 globalData->locaFactory->createBorderedSolverStrategy(parsedParams,
00087 turningPointParams);
00088 if (!isSymmetric) {
00089 b_vector = b->createMultiVector(1, NOX::DeepCopy);
00090 }
00091 else {
00092 b_vector = a_vector->clone(NOX::DeepCopy);
00093 }
00094
00095
00096 updateVectorsEveryContinuationStep =
00097 turningPointParams->get(
00098 "Update Null Vectors Every Continuation Step",
00099 true);
00100 updateVectorsEveryIteration =
00101 turningPointParams->get(
00102 "Update Null Vectors Every Nonlinear Iteration",
00103 false);
00104 }
00105
00106 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00107 Constraint(const LOCA::TurningPoint::MinimallyAugmented::Constraint& source,
00108 NOX::CopyType type) :
00109 globalData(source.globalData),
00110 parsedParams(source.parsedParams),
00111 turningPointParams(source.turningPointParams),
00112 grpPtr(Teuchos::null),
00113 a_vector(source.a_vector->clone(type)),
00114 b_vector(source.b_vector->clone(type)),
00115 w_vector(source.w_vector->clone(type)),
00116 v_vector(source.v_vector->clone(type)),
00117 Jv_vector(source.Jv_vector->clone(type)),
00118 sigma_x(source.sigma_x->clone(type)),
00119 constraints(source.constraints),
00120 borderedSolver(),
00121 dn(source.dn),
00122 sigma_scale(source.sigma_scale),
00123 isSymmetric(source.isSymmetric),
00124 isValidConstraints(false),
00125 isValidDX(false),
00126 bifParamID(source.bifParamID),
00127 updateVectorsEveryContinuationStep(source.updateVectorsEveryContinuationStep),
00128 updateVectorsEveryIteration(source.updateVectorsEveryIteration)
00129 {
00130 if (source.isValidConstraints && type == NOX::DeepCopy)
00131 isValidConstraints = true;
00132 if (source.isValidDX && type == NOX::DeepCopy)
00133 isValidDX = true;
00134
00135
00136 borderedSolver =
00137 globalData->locaFactory->createBorderedSolverStrategy(parsedParams,
00138 turningPointParams);
00139
00140
00141
00142 }
00143
00144 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00145 ~Constraint()
00146 {
00147 }
00148
00149 void
00150 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00151 setGroup(const Teuchos::RCP<LOCA::TurningPoint::MinimallyAugmented::AbstractGroup>& g)
00152 {
00153 grpPtr = g;
00154 }
00155
00156 Teuchos::RCP<const NOX::Abstract::Vector>
00157 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00158 getLeftNullVec() const
00159 {
00160 return Teuchos::rcp(&(*w_vector)[0], false);
00161 }
00162
00163 Teuchos::RCP<const NOX::Abstract::Vector>
00164 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00165 getRightNullVec() const
00166 {
00167 return Teuchos::rcp(&(*v_vector)[0], false);
00168 }
00169
00170 double
00171 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00172 getSigma() const
00173 {
00174 return constraints(0,0);
00175 }
00176
00177 void
00178 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00179 copy(const LOCA::MultiContinuation::ConstraintInterface& src)
00180 {
00181 const LOCA::TurningPoint::MinimallyAugmented::Constraint& source =
00182 dynamic_cast<const LOCA::TurningPoint::MinimallyAugmented::Constraint&>(src);
00183
00184 if (this != &source) {
00185 globalData = source.globalData;
00186 parsedParams = source.parsedParams;
00187 turningPointParams = source.turningPointParams;
00188 *a_vector = *(source.a_vector);
00189 *b_vector = *(source.b_vector);
00190 *w_vector = *(source.w_vector);
00191 *v_vector = *(source.v_vector);
00192 *Jv_vector = *(source.Jv_vector);
00193 *sigma_x = *(source.sigma_x);
00194 constraints.assign(source.constraints);
00195 dn = source.dn;
00196 sigma_scale = source.sigma_scale;
00197 isSymmetric = source.isSymmetric;
00198 isValidConstraints = source.isValidConstraints;
00199 isValidDX = source.isValidDX;
00200 bifParamID = source.bifParamID;
00201 updateVectorsEveryContinuationStep =
00202 source.updateVectorsEveryContinuationStep;
00203 updateVectorsEveryIteration =
00204 source.updateVectorsEveryIteration;
00205
00206
00207 borderedSolver =
00208 globalData->locaFactory->createBorderedSolverStrategy(
00209 parsedParams,
00210 turningPointParams);
00211
00212
00213
00214 }
00215 }
00216
00217 Teuchos::RCP<LOCA::MultiContinuation::ConstraintInterface>
00218 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00219 clone(NOX::CopyType type) const
00220 {
00221 return Teuchos::rcp(new Constraint(*this, type));
00222 }
00223
00224 int
00225 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00226 numConstraints() const
00227 {
00228 return 1;
00229 }
00230
00231 void
00232 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00233 setX(const NOX::Abstract::Vector& y)
00234 {
00235 grpPtr->setX(y);
00236 isValidConstraints = false;
00237 isValidDX = false;
00238 }
00239
00240 void
00241 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00242 setParam(int paramID, double val)
00243 {
00244 grpPtr->setParam(paramID, val);
00245 isValidConstraints = false;
00246 isValidDX = false;
00247 }
00248
00249 void
00250 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00251 setParams(const vector<int>& paramIDs,
00252 const NOX::Abstract::MultiVector::DenseMatrix& vals)
00253 {
00254 grpPtr->setParamsMulti(paramIDs, vals);
00255 isValidConstraints = false;
00256 isValidDX = false;
00257 }
00258
00259 NOX::Abstract::Group::ReturnType
00260 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00261 computeConstraints()
00262 {
00263 if (isValidConstraints)
00264 return NOX::Abstract::Group::Ok;
00265
00266 string callingFunction =
00267 "LOCA::TurningPoint::MinimallyAugmented::Constraint::computeConstraints()";
00268 NOX::Abstract::Group::ReturnType status;
00269 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00270
00271
00272 status = grpPtr->computeJacobian();
00273 finalStatus =
00274 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00275 finalStatus,
00276 callingFunction);
00277
00278
00279 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> op =
00280 Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(grpPtr));
00281 borderedSolver->setMatrixBlocksMultiVecConstraint(op,
00282 a_vector,
00283 b_vector,
00284 Teuchos::null);
00285
00286
00287 NOX::Abstract::MultiVector::DenseMatrix one(1,1);
00288 one(0,0) = dn;
00289
00290
00291 Teuchos::RCP<Teuchos::ParameterList> linear_solver_params =
00292 parsedParams->getSublist("Linear Solver");
00293
00294
00295 NOX::Abstract::MultiVector::DenseMatrix s1(1,1);
00296 status = borderedSolver->initForSolve();
00297 finalStatus =
00298 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00299 finalStatus,
00300 callingFunction);
00301 status = borderedSolver->applyInverse(*linear_solver_params,
00302 NULL,
00303 &one,
00304 *v_vector,
00305 s1);
00306 finalStatus =
00307 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00308 finalStatus,
00309 callingFunction);
00310
00311
00312 NOX::Abstract::MultiVector::DenseMatrix s2(1,1);
00313 if (!isSymmetric) {
00314 status = borderedSolver->initForTransposeSolve();
00315 finalStatus =
00316 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00317 finalStatus,
00318 callingFunction);
00319 status = borderedSolver->applyInverseTranspose(*linear_solver_params,
00320 NULL,
00321 &one,
00322 *w_vector,
00323 s2);
00324 finalStatus =
00325 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00326 finalStatus,
00327 callingFunction);
00328
00329 }
00330 else {
00331 *w_vector = *v_vector;
00332 s2.assign(s1);
00333 }
00334
00335
00336 status = grpPtr->applyJacobianMultiVector(*v_vector, *Jv_vector);
00337 finalStatus =
00338 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00339 finalStatus,
00340 callingFunction);
00341 Jv_vector->multiply(-1.0, *w_vector, constraints);
00342
00343
00344
00345
00346
00347 sigma_scale = dn;
00348 constraints.scale(1.0/sigma_scale);
00349
00350 if (globalData->locaUtils->isPrintType(NOX::Utils::OuterIteration)) {
00351 globalData->locaUtils->out() <<
00352 "\n\tEstimate for singularity of Jacobian (sigma1) = " <<
00353 globalData->locaUtils->sciformat(s1(0,0));
00354 globalData->locaUtils->out() <<
00355 "\n\tEstimate for singularity of Jacobian (sigma2) = " <<
00356 globalData->locaUtils->sciformat(s2(0,0));
00357 globalData->locaUtils->out() <<
00358 "\n\tEstimate for singularity of Jacobian (sigma) = " <<
00359 globalData->locaUtils->sciformat(constraints(0,0)) << std::endl;
00360 }
00361
00362 isValidConstraints = true;
00363
00364
00365 if (updateVectorsEveryIteration) {
00366 if (globalData->locaUtils->isPrintType(NOX::Utils::OuterIteration)) {
00367 globalData->locaUtils->out() <<
00368 "\n\tUpdating null vectors for the next nonlinear iteration" <<
00369 std::endl;
00370 }
00371 *a_vector = *w_vector;
00372 *b_vector = *v_vector;
00373
00374 a_vector->scale(std::sqrt(dn) / (*a_vector)[0].norm());
00375 b_vector->scale(std::sqrt(dn) / (*b_vector)[0].norm());
00376 }
00377
00378 return finalStatus;
00379 }
00380
00381 NOX::Abstract::Group::ReturnType
00382 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00383 computeDX()
00384 {
00385 if (isValidDX)
00386 return NOX::Abstract::Group::Ok;
00387
00388 string callingFunction =
00389 "LOCA::TurningPoint::MinimallyAugmented::Constraint::computeDX()";
00390 NOX::Abstract::Group::ReturnType status;
00391 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00392
00393
00394 if (!isValidConstraints) {
00395 status = computeConstraints();
00396 finalStatus =
00397 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00398 finalStatus,
00399 callingFunction);
00400 }
00401
00402
00403 status = grpPtr->computeDwtJnDx((*w_vector)[0], (*v_vector)[0],
00404 (*sigma_x)[0]);
00405 finalStatus =
00406 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00407 finalStatus,
00408 callingFunction);
00409 sigma_x->scale(-1.0/sigma_scale);
00410
00411 isValidDX = true;
00412
00413 return finalStatus;
00414 }
00415
00416 NOX::Abstract::Group::ReturnType
00417 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00418 computeDP(const vector<int>& paramIDs,
00419 NOX::Abstract::MultiVector::DenseMatrix& dgdp,
00420 bool isValidG)
00421 {
00422 string callingFunction =
00423 "LOCA::TurningPoint::MinimallyAugmented::Constraint::computeDP()";
00424 NOX::Abstract::Group::ReturnType status;
00425 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00426
00427
00428 if (!isValidConstraints) {
00429 status = computeConstraints();
00430 finalStatus =
00431 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00432 finalStatus,
00433 callingFunction);
00434 }
00435
00436
00437 status = grpPtr->computeDwtJnDp(paramIDs, (*w_vector)[0], (*v_vector)[0],
00438 dgdp, false);
00439 finalStatus =
00440 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00441 finalStatus,
00442 callingFunction);
00443 dgdp.scale(-1.0/sigma_scale);
00444
00445
00446 dgdp(0,0) = constraints(0,0);
00447
00448 return finalStatus;
00449 }
00450
00451 bool
00452 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00453 isConstraints() const
00454 {
00455 return isValidConstraints;
00456 }
00457
00458 bool
00459 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00460 isDX() const
00461 {
00462 return isValidDX;
00463 }
00464
00465 const NOX::Abstract::MultiVector::DenseMatrix&
00466 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00467 getConstraints() const
00468 {
00469 return constraints;
00470 }
00471
00472 const NOX::Abstract::MultiVector*
00473 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00474 getDX() const
00475 {
00476 return sigma_x.get();
00477 }
00478
00479 bool
00480 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00481 isDXZero() const
00482 {
00483 return false;
00484 }
00485
00486 void
00487 LOCA::TurningPoint::MinimallyAugmented::Constraint::
00488 postProcessContinuationStep(LOCA::Abstract::Iterator::StepStatus stepStatus)
00489 {
00490 if (stepStatus == LOCA::Abstract::Iterator::Successful &&
00491 updateVectorsEveryContinuationStep) {
00492 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00493 globalData->locaUtils->out() <<
00494 "\n\tUpdating null vectors for the next continuation step" << std::endl;
00495 }
00496 *a_vector = *w_vector;
00497 *b_vector = *v_vector;
00498
00499 a_vector->scale(std::sqrt(dn) / (*a_vector)[0].norm());
00500 b_vector->scale(std::sqrt(dn) / (*b_vector)[0].norm());
00501 }
00502 }