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 "Teuchos_ParameterList.hpp"
00043
00044 #include "LOCA_Homotopy_Group.H"
00045 #include "LOCA_Homotopy_AbstractGroup.H"
00046 #include "LOCA_Parameter_Vector.H"
00047 #include "LOCA_GlobalData.H"
00048 #include "LOCA_ErrorCheck.H"
00049 #include "NOX_Utils.H"
00050
00051 LOCA::Homotopy::Group::Group(
00052 Teuchos::ParameterList& locaSublist,
00053 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00054 const Teuchos::RCP<LOCA::Homotopy::AbstractGroup>& g,
00055 double scalarRandom,
00056 double scalarInitialGuess) :
00057 globalData(global_data),
00058 grpPtr(g),
00059 gVecPtr(g->getX().clone(NOX::ShapeCopy)),
00060 randomVecPtr(gVecPtr->clone(NOX::ShapeCopy)),
00061 newtonVecPtr(),
00062 gradVecPtr(),
00063 paramVec(grpPtr->getParams()),
00064 conParam(0.0),
00065 conParamID(-1),
00066 conParamLabel("Homotopy Continuation Parameter"),
00067 augmentJacForHomotopyNotImplemented(false)
00068 {
00069
00070 randomVecPtr->random();
00071 randomVecPtr->abs(*randomVecPtr);
00072 randomVecPtr->update(scalarInitialGuess, grpPtr->getX(), scalarRandom);
00073
00074
00075 resetIsValidFlags();
00076
00077
00078
00079
00080 paramVec.addParameter(conParamLabel, conParam);
00081 grpPtr->setParams(paramVec);
00082
00083
00084 conParamID = paramVec.getIndex(conParamLabel);
00085
00086 setStepperParameters(locaSublist);
00087
00088 }
00089
00090 LOCA::Homotopy::Group::Group(
00091 Teuchos::ParameterList& locaSublist,
00092 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00093 const Teuchos::RCP<LOCA::Homotopy::AbstractGroup>& g,
00094 const NOX::Abstract::Vector& randomVector) :
00095 globalData(global_data),
00096 grpPtr(g),
00097 gVecPtr(g->getX().clone(NOX::ShapeCopy)),
00098 randomVecPtr(gVecPtr->clone(NOX::ShapeCopy)),
00099 newtonVecPtr(),
00100 gradVecPtr(),
00101 paramVec(grpPtr->getParams()),
00102 conParam(0.0),
00103 conParamID(-1),
00104 conParamLabel("Homotopy Continuation Parameter"),
00105 augmentJacForHomotopyNotImplemented(false)
00106 {
00107
00108 *randomVecPtr = randomVector;
00109
00110
00111 resetIsValidFlags();
00112
00113
00114
00115
00116 paramVec.addParameter(conParamLabel, conParam);
00117 grpPtr->setParams(paramVec);
00118
00119
00120 conParamID = paramVec.getIndex(conParamLabel);
00121
00122 setStepperParameters(locaSublist);
00123
00124 }
00125
00126 LOCA::Homotopy::Group::Group(const LOCA::Homotopy::Group& source,
00127 NOX::CopyType type) :
00128 globalData(source.globalData),
00129 grpPtr(Teuchos::rcp_dynamic_cast<LOCA::Homotopy::AbstractGroup>(source.grpPtr->clone(type))),
00130 gVecPtr(source.gVecPtr->clone(type)),
00131 randomVecPtr(source.randomVecPtr->clone(NOX::DeepCopy)),
00132 newtonVecPtr(),
00133 gradVecPtr(),
00134 paramVec(source.paramVec),
00135 conParam(source.conParam),
00136 conParamID(source.conParamID),
00137 conParamLabel(source.conParamLabel),
00138 augmentJacForHomotopyNotImplemented(source.augmentJacForHomotopyNotImplemented)
00139 {
00140 if (source.newtonVecPtr != Teuchos::null)
00141 newtonVecPtr = source.newtonVecPtr->clone(type);
00142
00143 if (source.gradVecPtr != Teuchos::null)
00144 newtonVecPtr = source.gradVecPtr->clone(type);
00145
00146 if (type == NOX::DeepCopy) {
00147 isValidF = source.isValidF;
00148 isValidJacobian = source.isValidJacobian;
00149 isValidNewton = source.isValidNewton;
00150 isValidGradient = source.isValidGradient;
00151 }
00152 else if (type == NOX::ShapeCopy) {
00153 resetIsValidFlags();
00154 }
00155 else {
00156 globalData->locaErrorCheck->throwError(
00157 "LOCA::Homotopy::Group::Group(copy ctor)",
00158 "CopyType is invalid!");
00159 }
00160
00161 }
00162
00163
00164 LOCA::Homotopy::Group::~Group()
00165 {
00166 }
00167
00168 NOX::Abstract::Group&
00169 LOCA::Homotopy::Group::operator=(const NOX::Abstract::Group& source)
00170 {
00171 copy(source);
00172 return *this;
00173 }
00174
00175 Teuchos::RCP<NOX::Abstract::Group>
00176 LOCA::Homotopy::Group::clone(NOX::CopyType type) const
00177 {
00178 return Teuchos::rcp(new LOCA::Homotopy::Group(*this, type));
00179 }
00180
00181 void
00182 LOCA::Homotopy::Group::setX(const NOX::Abstract::Vector& y)
00183 {
00184 resetIsValidFlags();
00185 grpPtr->setX(y);
00186 }
00187
00188 void
00189 LOCA::Homotopy::Group::computeX(const NOX::Abstract::Group& grp,
00190 const NOX::Abstract::Vector& d,
00191 double step)
00192 {
00193 resetIsValidFlags();
00194 const LOCA::Homotopy::Group& g =
00195 dynamic_cast<const LOCA::Homotopy::Group&>(grp);
00196 grpPtr->computeX(*(g.grpPtr), d, step);
00197 }
00198
00199 NOX::Abstract::Group::ReturnType
00200 LOCA::Homotopy::Group::computeF()
00201 {
00202 if (isValidF)
00203 return NOX::Abstract::Group::Ok;
00204
00205
00206 NOX::Abstract::Group::ReturnType status = grpPtr->computeF();
00207
00208
00209 *gVecPtr = grpPtr->getX();
00210 gVecPtr->update(-1.0, *randomVecPtr, 1.0);
00211 gVecPtr->scale( 1.0 - conParam);
00212 gVecPtr->update(conParam, grpPtr->getF(), 1.0);
00213
00214 isValidF = true;
00215
00216 return status;
00217 }
00218
00219 NOX::Abstract::Group::ReturnType
00220 LOCA::Homotopy::Group::computeJacobian()
00221 {
00222 if (isValidJacobian)
00223 return NOX::Abstract::Group::Ok;
00224
00225
00226 NOX::Abstract::Group::ReturnType status = grpPtr->computeJacobian();
00227
00228
00229 NOX::Abstract::Group::ReturnType augHomTest =
00230 grpPtr->augmentJacobianForHomotopy(conParam, 1.0-conParam);
00231
00232
00233 if (augHomTest == NOX::Abstract::Group::NotDefined)
00234 augmentJacForHomotopyNotImplemented = true;
00235
00236 isValidJacobian = true;
00237
00238 return status;
00239 }
00240
00241 NOX::Abstract::Group::ReturnType
00242 LOCA::Homotopy::Group::computeGradient()
00243 {
00244 if (isValidGradient)
00245 return NOX::Abstract::Group::Ok;
00246
00247 string callingFunction =
00248 "LOCA::Homotopy::Group::computeGradient()";
00249 NOX::Abstract::Group::ReturnType status, finalStatus;
00250
00251 if (gradVecPtr == Teuchos::null)
00252 gradVecPtr = gVecPtr->clone(NOX::ShapeCopy);
00253
00254 finalStatus = computeF();
00255 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00256
00257 status = computeJacobian();
00258 finalStatus =
00259 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00260 callingFunction);
00261
00262 status = applyJacobianTranspose(*gVecPtr, *gradVecPtr);
00263 finalStatus =
00264 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00265 callingFunction);
00266
00267 return finalStatus;
00268 }
00269
00270 NOX::Abstract::Group::ReturnType
00271 LOCA::Homotopy::Group::computeNewton(Teuchos::ParameterList& params)
00272 {
00273 if (isValidNewton)
00274 return NOX::Abstract::Group::Ok;
00275
00276 string callingFunction =
00277 "LOCA::Homotopy::Group::computeNewton()";
00278 NOX::Abstract::Group::ReturnType status, finalStatus;
00279
00280 if (newtonVecPtr == Teuchos::null)
00281 newtonVecPtr = gVecPtr->clone(NOX::ShapeCopy);
00282
00283 finalStatus = computeF();
00284 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00285
00286 status = computeJacobian();
00287 finalStatus =
00288 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00289 callingFunction);
00290
00291 status = applyJacobianInverse(params, *gVecPtr, *newtonVecPtr);
00292 finalStatus =
00293 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00294 callingFunction);
00295
00296 newtonVecPtr->scale(-1.0);
00297
00298 isValidNewton = true;
00299
00300 return finalStatus;
00301 }
00302
00303 NOX::Abstract::Group::ReturnType
00304 LOCA::Homotopy::Group::applyJacobian(const NOX::Abstract::Vector& input,
00305 NOX::Abstract::Vector& result) const
00306 {
00307 if (!isValidJacobian)
00308 return NOX::Abstract::Group::BadDependency;
00309
00310 NOX::Abstract::Group::ReturnType status =
00311 grpPtr->applyJacobian(input, result);
00312
00313
00314
00315 if (augmentJacForHomotopyNotImplemented)
00316 result.update(1.0-conParam, input, conParam);
00317
00318 return status;
00319 }
00320
00321 NOX::Abstract::Group::ReturnType
00322 LOCA::Homotopy::Group::applyJacobianTranspose(
00323 const NOX::Abstract::Vector& input,
00324 NOX::Abstract::Vector& result) const
00325 {
00326 if (!isValidJacobian)
00327 return NOX::Abstract::Group::BadDependency;
00328
00329 NOX::Abstract::Group::ReturnType status =
00330 grpPtr->applyJacobianTranspose(input, result);
00331
00332
00333
00334 if (augmentJacForHomotopyNotImplemented)
00335 result.update(1.0-conParam, input, conParam);
00336
00337 return status;
00338 }
00339
00340 NOX::Abstract::Group::ReturnType
00341 LOCA::Homotopy::Group::applyJacobianInverse(Teuchos::ParameterList& params,
00342 const NOX::Abstract::Vector& input,
00343 NOX::Abstract::Vector& result) const
00344 {
00345 return grpPtr->applyJacobianInverse(params, input, result);
00346 }
00347
00348 NOX::Abstract::Group::ReturnType
00349 LOCA::Homotopy::Group::applyJacobianMultiVector(
00350 const NOX::Abstract::MultiVector& input,
00351 NOX::Abstract::MultiVector& result) const
00352 {
00353 if (!isValidJacobian)
00354 return NOX::Abstract::Group::BadDependency;
00355
00356 NOX::Abstract::Group::ReturnType status =
00357 grpPtr->applyJacobianMultiVector(input, result);
00358
00359
00360
00361 if (augmentJacForHomotopyNotImplemented)
00362 result.update(1.0-conParam, input, conParam);
00363
00364 return status;
00365 }
00366
00367 NOX::Abstract::Group::ReturnType
00368 LOCA::Homotopy::Group::applyJacobianTransposeMultiVector(
00369 const NOX::Abstract::MultiVector& input,
00370 NOX::Abstract::MultiVector& result) const
00371 {
00372 if (!isValidJacobian)
00373 return NOX::Abstract::Group::BadDependency;
00374
00375 NOX::Abstract::Group::ReturnType status =
00376 grpPtr->applyJacobianTransposeMultiVector(input, result);
00377
00378
00379
00380 if (augmentJacForHomotopyNotImplemented)
00381 result.update(1.0-conParam, input, conParam);
00382
00383 return status;
00384 }
00385
00386 NOX::Abstract::Group::ReturnType
00387 LOCA::Homotopy::Group::applyJacobianInverseMultiVector(
00388 Teuchos::ParameterList& params,
00389 const NOX::Abstract::MultiVector& input,
00390 NOX::Abstract::MultiVector& result) const
00391 {
00392 return grpPtr->applyJacobianInverseMultiVector(params, input, result);
00393 }
00394
00395 bool
00396 LOCA::Homotopy::Group::isF() const
00397 {
00398 return isValidF;
00399 }
00400
00401 bool
00402 LOCA::Homotopy::Group::isJacobian() const
00403 {
00404 return isValidJacobian;
00405 }
00406
00407 bool
00408 LOCA::Homotopy::Group::isNewton() const
00409 {
00410 return isValidNewton;
00411 }
00412
00413 bool
00414 LOCA::Homotopy::Group::isGradient() const
00415 {
00416 return isValidGradient;
00417 }
00418
00419 const NOX::Abstract::Vector&
00420 LOCA::Homotopy::Group::getX() const
00421 {
00422 return grpPtr->getX();
00423 }
00424
00425 const NOX::Abstract::Vector&
00426 LOCA::Homotopy::Group::getF() const
00427 {
00428 return *gVecPtr;
00429 }
00430
00431 double
00432 LOCA::Homotopy::Group::getNormF() const
00433 {
00434 return gVecPtr->norm();
00435 }
00436
00437 const NOX::Abstract::Vector&
00438 LOCA::Homotopy::Group::getGradient() const
00439 {
00440 if (gradVecPtr == Teuchos::null) {
00441 globalData->locaErrorCheck->throwError(
00442 "LOCA::Homotopy::Group::getGradient",
00443 "gradVecPtr is NULL!");
00444 }
00445 return *gradVecPtr;
00446 }
00447
00448 const NOX::Abstract::Vector&
00449 LOCA::Homotopy::Group::getNewton() const
00450 {
00451 if (newtonVecPtr == Teuchos::null) {
00452 globalData->locaErrorCheck->throwError("LOCA::Homotopy::Group::getNewton",
00453 "newtonVecPtr is NULL!");
00454 }
00455 return *newtonVecPtr;
00456 }
00457
00458 Teuchos::RCP<const LOCA::MultiContinuation::AbstractGroup>
00459 LOCA::Homotopy::Group::getUnderlyingGroup() const
00460 {
00461 return grpPtr;
00462 }
00463
00464 Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup>
00465 LOCA::Homotopy::Group::getUnderlyingGroup()
00466 {
00467 return grpPtr;
00468 }
00469
00470
00471 void
00472 LOCA::Homotopy::Group::copy(const NOX::Abstract::Group& src)
00473 {
00474
00475 const LOCA::Homotopy::Group& source =
00476 dynamic_cast<const LOCA::Homotopy::Group&>(src);
00477
00478
00479 if (this != &source) {
00480
00481 globalData = source.globalData;
00482 grpPtr->copy(*(source.grpPtr));
00483 *gVecPtr = *(source.gVecPtr);
00484 *randomVecPtr = *(source.randomVecPtr);
00485 if (newtonVecPtr != Teuchos::null)
00486 *newtonVecPtr = *(source.newtonVecPtr);
00487 if (gradVecPtr != Teuchos::null)
00488 *gradVecPtr = *(source.gradVecPtr);
00489 isValidF = source.isValidF;
00490 isValidJacobian = source.isValidJacobian;
00491 isValidNewton = source.isValidNewton;
00492 isValidGradient = source.isValidGradient;
00493 paramVec = source.paramVec;
00494 conParam = source.conParam;
00495 conParamID = source.conParamID;
00496 augmentJacForHomotopyNotImplemented =
00497 source.augmentJacForHomotopyNotImplemented;
00498 }
00499
00500 }
00501
00502 void
00503 LOCA::Homotopy::Group::setParamsMulti(
00504 const vector<int>& paramIDs,
00505 const NOX::Abstract::MultiVector::DenseMatrix& vals)
00506 {
00507 resetIsValidFlags();
00508 grpPtr->setParamsMulti(paramIDs, vals);
00509 for (unsigned int i=0; i<paramIDs.size(); i++)
00510 if (paramIDs[i] == conParamID)
00511 conParam = vals(i,0);
00512 }
00513
00514 void
00515 LOCA::Homotopy::Group::setParams(const LOCA::ParameterVector& p)
00516 {
00517 resetIsValidFlags();
00518 grpPtr->setParams(p);
00519 conParam = p.getValue(conParamLabel);
00520 }
00521
00522 void
00523 LOCA::Homotopy::Group::setParam(int paramID, double val)
00524 {
00525 resetIsValidFlags();
00526 grpPtr->setParam(paramID, val);
00527 if (paramID == conParamID)
00528 conParam = val;
00529 }
00530
00531 void
00532 LOCA::Homotopy::Group::setParam(string paramID, double val)
00533 {
00534 resetIsValidFlags();
00535 grpPtr->setParam(paramID, val);
00536 if (paramID == conParamLabel)
00537 conParam = val;
00538 }
00539
00540 const LOCA::ParameterVector&
00541 LOCA::Homotopy::Group::getParams() const
00542 {
00543 return grpPtr->getParams();
00544 }
00545
00546 double
00547 LOCA::Homotopy::Group::getParam(int paramID) const
00548 {
00549 return grpPtr->getParam(paramID);
00550 }
00551
00552 double
00553 LOCA::Homotopy::Group::getParam(string paramID) const
00554 {
00555 return grpPtr->getParam(paramID);
00556 }
00557
00558 NOX::Abstract::Group::ReturnType
00559 LOCA::Homotopy::Group::computeDfDpMulti(const vector<int>& paramIDs,
00560 NOX::Abstract::MultiVector& dfdp,
00561 bool isValidF)
00562 {
00563
00564
00565
00566
00567
00568
00569
00570
00571 vector<int> pIDs;
00572 vector<int> idx(1);
00573 idx[0] = 0;
00574 for (unsigned int i=0; i<paramIDs.size(); i++)
00575 if (paramIDs[i] != conParamID) {
00576 pIDs.push_back(paramIDs[i]);
00577 idx.push_back(i+1);
00578 }
00579
00580
00581 Teuchos::RCP<NOX::Abstract::MultiVector> fp =
00582 dfdp.subView(idx);
00583
00584
00585
00586 NOX::Abstract::Group::ReturnType status =
00587 grpPtr->computeDfDpMulti(pIDs, *fp, false);
00588
00589
00590 fp->scale(conParam);
00591
00592
00593 double v = 1.0-conParam;
00594 dfdp[0].update(v, grpPtr->getX(), -v, *randomVecPtr, 1.0);
00595
00596
00597 grpPtr->computeF();
00598 for (unsigned int i=0; i<paramIDs.size(); i++)
00599 if (paramIDs[i] == conParamID) {
00600 dfdp[i+1] = grpPtr->getF();
00601 dfdp[i+1].update(-1.0, grpPtr->getX(), 1.0, *randomVecPtr, 1.0);
00602 }
00603
00604 return status;
00605 }
00606
00607 void
00608 LOCA::Homotopy::Group::preProcessContinuationStep(
00609 LOCA::Abstract::Iterator::StepStatus stepStatus)
00610 {
00611 grpPtr->preProcessContinuationStep(stepStatus);
00612 }
00613
00614 void
00615 LOCA::Homotopy::Group::postProcessContinuationStep(
00616 LOCA::Abstract::Iterator::StepStatus stepStatus)
00617 {
00618 grpPtr->postProcessContinuationStep(stepStatus);
00619 }
00620
00621 void
00622 LOCA::Homotopy::Group::projectToDraw(const NOX::Abstract::Vector& x,
00623 double *px) const
00624 {
00625 grpPtr->projectToDraw(x, px);
00626 px[this->projectToDrawDimension()] = conParam;
00627 }
00628
00629 int
00630 LOCA::Homotopy::Group::projectToDrawDimension() const
00631 {
00632 return grpPtr->projectToDrawDimension()+1;
00633 }
00634
00635 void
00636 LOCA::Homotopy::Group::printSolution(const double conParm) const
00637 {
00638 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00639 globalData->locaUtils->out() <<
00640 "\tPrinting Solution Vector for homotopy parameter = " <<
00641 globalData->locaUtils->sciformat(conParam) << std::endl;
00642 }
00643 grpPtr->printSolution(conParam);
00644 return;
00645 }
00646
00647 void
00648 LOCA::Homotopy::Group::printSolution(const NOX::Abstract::Vector& x_,
00649 const double conParm) const
00650 {
00651 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00652 globalData->locaUtils->out() <<
00653 "\tPrinting Solution Vector for homotopy parameter = " <<
00654 globalData->locaUtils->sciformat(conParam) << std::endl;
00655 }
00656 grpPtr->printSolution(conParam);
00657 return;
00658 }
00659
00660 void
00661 LOCA::Homotopy::Group::resetIsValidFlags()
00662 {
00663 isValidF = false;
00664 isValidJacobian = false;
00665 isValidNewton = false;
00666 isValidGradient = false;
00667 return;
00668 }
00669
00670 void
00671 LOCA::Homotopy::Group::setStepperParameters(Teuchos::ParameterList& params)
00672 {
00673
00674
00675 Teuchos::ParameterList& stepperList = params.sublist("Stepper");
00676 stepperList.set("Continuation Method", "Natural");
00677 stepperList.set("Continuation Parameter", conParamLabel);
00678 stepperList.set("Initial Value", 0.0);
00679 stepperList.set("Max Value", 1.0);
00680 stepperList.set("Min Value", -1.0);
00681 stepperList.set("Max Steps", 50);
00682
00683
00684 Teuchos::ParameterList& predictorList = params.sublist("Predictor");
00685 predictorList.set("Method", "Constant");
00686
00687
00688 Teuchos::ParameterList& stepSizeList = params.sublist("Step Size");
00689
00690 stepSizeList.set("Method", "Adaptive");
00691 stepSizeList.set("Initial Step Size", 0.1);
00692 stepSizeList.set("Min Step Size", 1.0e-2);
00693 stepSizeList.set("Max Step Size", 1.0);
00694 stepSizeList.set("Aggressiveness", 0.5);
00695 return;
00696 }