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_Pitchfork_MinimallyAugmented_ExtendedGroup.H"
00043
00044 #include "Teuchos_ParameterList.hpp"
00045 #include "LOCA_Pitchfork_MinimallyAugmented_AbstractGroup.H"
00046 #include "LOCA_Pitchfork_MinimallyAugmented_Constraint.H"
00047 #include "LOCA_GlobalData.H"
00048 #include "LOCA_Factory.H"
00049 #include "LOCA_Parameter_SublistParser.H"
00050 #include "LOCA_BorderedSolver_AbstractStrategy.H"
00051 #include "LOCA_ErrorCheck.H"
00052 #include "NOX_Utils.H"
00053 #include "LOCA_Parameter_Vector.H"
00054 #include "LOCA_Abstract_TransposeSolveGroup.H"
00055 #include "LOCA_BorderedSolver_JacobianOperator.H"
00056
00057 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00058 ExtendedGroup(
00059 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00060 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00061 const Teuchos::RCP<Teuchos::ParameterList>& pfParams,
00062 const Teuchos::RCP<LOCA::Pitchfork::MinimallyAugmented::AbstractGroup>& g)
00063 : globalData(global_data),
00064 parsedParams(topParams),
00065 pitchforkParams(pfParams),
00066 grpPtr(g),
00067 bordered_grp(),
00068 constraintsPtr(),
00069 xMultiVec(globalData, g->getX(), 1, 2, NOX::DeepCopy),
00070 fMultiVec(globalData, g->getX(), 3, 2, NOX::ShapeCopy),
00071 newtonMultiVec(globalData, g->getX(), 1, 2, NOX::ShapeCopy),
00072 gradientMultiVec(globalData, g->getX(), 1, 2, NOX::ShapeCopy),
00073 xVec(),
00074 psiVec(),
00075 fVec(),
00076 ffMultiVec(),
00077 dfdpMultiVec(),
00078 fBifMultiVec(),
00079 newtonVec(),
00080 gradientVec(),
00081 jacOp(),
00082 borderedSolver(),
00083 index_f(1),
00084 index_dfdp(2),
00085 bifParamID(-1),
00086 isValidF(false),
00087 isValidJacobian(false),
00088 isValidNewton(false),
00089 isValidGradient(false),
00090 isBordered(false)
00091 {
00092 const char *func = "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup()";
00093
00094
00095 setupViews();
00096
00097
00098 if (!pitchforkParams->isParameter("Bifurcation Parameter")) {
00099 globalData->locaErrorCheck->throwError(func,
00100 "\"Bifurcation Parameter\" name is not set!");
00101 }
00102 string bifParamName = pitchforkParams->get("Bifurcation Parameter",
00103 "None");
00104 const ParameterVector& p = grpPtr->getParams();
00105 bifParamID = p.getIndex(bifParamName);
00106
00107
00108 if (!pitchforkParams->isParameter("Antisymmetric Vector")) {
00109 globalData->locaErrorCheck->throwError(func,
00110 "\"Antisymmetric Vector\" is not set!");
00111 }
00112 psiVec =
00113 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
00114 get< Teuchos::RCP<NOX::Abstract::Vector> >("Antisymmetric Vector");
00115
00116
00117 bool isSymmetric = pitchforkParams->get("Symmetric Jacobian", false);
00118
00119
00120 Teuchos::RCP<NOX::Abstract::Vector> aVecPtr;
00121 Teuchos::RCP<NOX::Abstract::Vector> bVecPtr;
00122 getInitialVectors(aVecPtr, bVecPtr, isSymmetric);
00123
00124
00125 constraintsPtr =
00126 Teuchos::rcp(new LOCA::Pitchfork::MinimallyAugmented::Constraint(
00127 globalData,
00128 parsedParams,
00129 pfParams,
00130 grpPtr,
00131 isSymmetric,
00132 *aVecPtr,
00133 bVecPtr.get(),
00134 psiVec,
00135 bifParamID));
00136
00137
00138 xVec->getScalar(0) = grpPtr->getParam(bifParamID);
00139 xVec->getScalar(1) = 0.0;
00140
00141
00142 constraintsPtr->setParam(bifParamID, xVec->getScalar(0));
00143 constraintsPtr->setX(*(xVec->getXVec()));
00144
00145
00146 borderedSolver = globalData->locaFactory->createBorderedSolverStrategy(
00147 parsedParams,
00148 pitchforkParams);
00149
00150
00151 bordered_grp =
00152 Teuchos::rcp_dynamic_cast<LOCA::BorderedSystem::AbstractGroup>(grpPtr);
00153 isBordered = (bordered_grp != Teuchos::null);
00154
00155
00156 jacOp = Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(grpPtr));
00157 }
00158
00159 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00160 ExtendedGroup(const LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup& source,
00161 NOX::CopyType type)
00162 : globalData(source.globalData),
00163 parsedParams(source.parsedParams),
00164 pitchforkParams(source.pitchforkParams),
00165 grpPtr(Teuchos::rcp_dynamic_cast<LOCA::Pitchfork::MinimallyAugmented::AbstractGroup>(source.grpPtr->clone(type))),
00166 bordered_grp(),
00167 constraintsPtr(Teuchos::rcp_dynamic_cast<LOCA::Pitchfork::MinimallyAugmented::Constraint>(source.constraintsPtr->clone(type))),
00168 xMultiVec(source.xMultiVec, type),
00169 fMultiVec(source.fMultiVec, type),
00170 newtonMultiVec(source.newtonMultiVec, type),
00171 gradientMultiVec(source.gradientMultiVec, type),
00172 xVec(),
00173 psiVec(source.psiVec),
00174 fVec(),
00175 ffMultiVec(),
00176 dfdpMultiVec(),
00177 fBifMultiVec(),
00178 newtonVec(),
00179 gradientVec(),
00180 jacOp(),
00181 borderedSolver(source.borderedSolver),
00182 index_f(1),
00183 index_dfdp(2),
00184 bifParamID(source.bifParamID),
00185 isValidF(source.isValidF),
00186 isValidJacobian(source.isValidJacobian),
00187 isValidNewton(source.isValidNewton),
00188 isValidGradient(source.isValidGradient),
00189 isBordered(false)
00190 {
00191
00192 setupViews();
00193
00194
00195 borderedSolver = globalData->locaFactory->createBorderedSolverStrategy(
00196 parsedParams,
00197 pitchforkParams);
00198
00199 if (type == NOX::ShapeCopy) {
00200 isValidF = false;
00201 isValidJacobian = false;
00202 isValidNewton = false;
00203 isValidGradient = false;
00204 }
00205
00206
00207 bordered_grp =
00208 Teuchos::rcp_dynamic_cast<LOCA::BorderedSystem::AbstractGroup>(grpPtr);
00209 isBordered = (bordered_grp != Teuchos::null);
00210
00211 constraintsPtr->setGroup(grpPtr);
00212
00213
00214 jacOp = Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(grpPtr));
00215
00216
00217 if (isValidJacobian) {
00218 borderedSolver->setMatrixBlocks(jacOp,
00219 dfdpMultiVec->getXMultiVec(),
00220 constraintsPtr,
00221 dfdpMultiVec->getScalars());
00222 NOX::Abstract::Group::ReturnType status = borderedSolver->initForSolve();
00223 globalData->locaErrorCheck->checkReturnType(status,
00224 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup()");
00225 }
00226 }
00227
00228
00229 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00230 ~ExtendedGroup()
00231 {
00232 }
00233
00234 double
00235 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00236 getBifParam() const
00237 {
00238 return grpPtr->getParam(bifParamID);
00239 }
00240
00241 NOX::Abstract::Group&
00242 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00243 operator=(const NOX::Abstract::Group& source)
00244 {
00245 copy(source);
00246 return *this;
00247 }
00248
00249 Teuchos::RCP<NOX::Abstract::Group>
00250 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00251 clone(NOX::CopyType type) const
00252 {
00253 return Teuchos::rcp(new ExtendedGroup(*this, type));
00254 }
00255
00256 void
00257 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00258 setX(const NOX::Abstract::Vector& y)
00259 {
00260 const LOCA::MultiContinuation::ExtendedVector& my =
00261 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(y);
00262
00263 grpPtr->setX( *(my.getXVec()) );
00264 grpPtr->setParam(bifParamID, my.getScalar(0));
00265 *xVec = my;
00266 constraintsPtr->setX( *(my.getXVec()) );
00267 constraintsPtr->setParam(bifParamID, my.getScalar(0));
00268
00269 resetIsValid();
00270 }
00271
00272 void
00273 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00274 computeX(const NOX::Abstract::Group& g,
00275 const NOX::Abstract::Vector& d,
00276 double step)
00277 {
00278 const LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup& mg =
00279 dynamic_cast<const LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup&>(g);
00280 const LOCA::MultiContinuation::ExtendedVector& md =
00281 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(d);
00282
00283 grpPtr->computeX(*(mg.grpPtr), *(md.getXVec()), step);
00284 xVec->update(1.0, mg.getX(), step, md, 0.0);
00285 grpPtr->setParam(bifParamID, xVec->getScalar(0));
00286 constraintsPtr->setX( *(xVec->getXVec()) );
00287 constraintsPtr->setParam(bifParamID, xVec->getScalar(0));
00288
00289 resetIsValid();
00290 }
00291
00292 NOX::Abstract::Group::ReturnType
00293 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00294 computeF()
00295 {
00296 if (isValidF)
00297 return NOX::Abstract::Group::Ok;
00298
00299 string callingFunction =
00300 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::computeF()";
00301 NOX::Abstract::Group::ReturnType status;
00302 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00303
00304
00305 if (!grpPtr->isF()) {
00306 status = grpPtr->computeF();
00307 finalStatus =
00308 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00309 finalStatus,
00310 callingFunction);
00311 }
00312 fVec->getXVec()->update(1.0, grpPtr->getF(), xVec->getScalar(1), *psiVec,
00313 0.0);
00314
00315
00316 if (!constraintsPtr->isConstraints()) {
00317 status = constraintsPtr->computeConstraints();
00318 }
00319 fVec->getScalars()->assign(constraintsPtr->getConstraints());
00320
00321 isValidF = true;
00322
00323 return finalStatus;
00324 }
00325
00326 NOX::Abstract::Group::ReturnType
00327 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00328 computeJacobian()
00329 {
00330 if (isValidJacobian)
00331 return NOX::Abstract::Group::Ok;
00332
00333 string callingFunction =
00334 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::computeJacobian()";
00335 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00336 NOX::Abstract::Group::ReturnType status;
00337
00338
00339
00340
00341 std::vector<int> paramIDs(1);
00342 paramIDs[0] = bifParamID;
00343 status = grpPtr->computeDfDpMulti(paramIDs,
00344 *(fBifMultiVec->getXMultiVec()),
00345 false);
00346 finalStatus =
00347 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00348 finalStatus,
00349 callingFunction);
00350
00351
00352 fVec->getXVec()->update(xVec->getScalar(1), *psiVec, 1.0);
00353
00354
00355 (*(dfdpMultiVec->getXMultiVec()))[1] = *psiVec;
00356
00357
00358
00359
00360
00361
00362 if (!constraintsPtr->isDX()) {
00363 status = constraintsPtr->computeDX();
00364 finalStatus =
00365 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00366 finalStatus,
00367 callingFunction);
00368 }
00369 status = constraintsPtr->computeDP(paramIDs,
00370 *(fBifMultiVec->getScalars()),
00371 isValidF);
00372 finalStatus =
00373 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00374 finalStatus,
00375 callingFunction);
00376
00377
00378 dfdpMultiVec->getScalar(0,1) = 0.0;
00379 dfdpMultiVec->getScalar(1,1) = 0.0;
00380
00381
00382 if (!grpPtr->isJacobian()) {
00383 status = grpPtr->computeJacobian();
00384 finalStatus =
00385 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00386 finalStatus,
00387 callingFunction);
00388 }
00389
00390
00391 borderedSolver->setMatrixBlocks(jacOp,
00392 dfdpMultiVec->getXMultiVec(),
00393 constraintsPtr,
00394 dfdpMultiVec->getScalars());
00395 status = borderedSolver->initForSolve();
00396 finalStatus =
00397 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00398 finalStatus,
00399 callingFunction);
00400
00401 isValidJacobian = true;
00402
00403 return finalStatus;
00404 }
00405
00406 NOX::Abstract::Group::ReturnType
00407 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00408 computeGradient()
00409 {
00410 if (isValidGradient)
00411 return NOX::Abstract::Group::Ok;
00412
00413 string callingFunction =
00414 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::computeGradient()";
00415 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00416 NOX::Abstract::Group::ReturnType status;
00417
00418
00419 if (!isF()) {
00420 status = computeF();
00421 finalStatus =
00422 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00423 finalStatus,
00424 callingFunction);
00425 }
00426
00427
00428 if (!isJacobian()) {
00429 status = computeJacobian();
00430 finalStatus =
00431 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00432 finalStatus,
00433 callingFunction);
00434 }
00435
00436
00437 status = applyJacobianTranspose(*fVec, *gradientVec);
00438 finalStatus =
00439 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00440 finalStatus,
00441 callingFunction);
00442
00443 isValidGradient = true;
00444
00445 return finalStatus;
00446 }
00447
00448 NOX::Abstract::Group::ReturnType
00449 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00450 computeNewton(Teuchos::ParameterList& params)
00451 {
00452 if (isValidNewton)
00453 return NOX::Abstract::Group::Ok;
00454
00455 string callingFunction =
00456 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::computeNewton()";
00457 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00458 NOX::Abstract::Group::ReturnType status;
00459
00460
00461 if (!isF()) {
00462 status = computeF();
00463 finalStatus =
00464 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00465 finalStatus,
00466 callingFunction);
00467 }
00468
00469
00470 if (!isJacobian()) {
00471 status = computeJacobian();
00472 finalStatus =
00473 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00474 finalStatus,
00475 callingFunction);
00476 }
00477
00478
00479 newtonMultiVec.init(0.0);
00480
00481 status = applyJacobianInverseMultiVector(params, *ffMultiVec,
00482 newtonMultiVec);
00483 finalStatus =
00484 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00485 finalStatus,
00486 callingFunction);
00487
00488 newtonMultiVec.scale(-1.0);
00489
00490 isValidNewton = true;
00491
00492 return finalStatus;
00493 }
00494
00495 NOX::Abstract::Group::ReturnType
00496 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00497 applyJacobian(const NOX::Abstract::Vector& input,
00498 NOX::Abstract::Vector& result) const
00499 {
00500
00501 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00502 input.createMultiVector(1, NOX::DeepCopy);
00503 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00504 result.createMultiVector(1, NOX::DeepCopy);
00505
00506
00507 NOX::Abstract::Group::ReturnType status =
00508 applyJacobianMultiVector(*mv_input, *mv_result);
00509
00510
00511 result = (*mv_result)[0];
00512
00513 return status;
00514 }
00515
00516 NOX::Abstract::Group::ReturnType
00517 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00518 applyJacobianTranspose(const NOX::Abstract::Vector& input,
00519 NOX::Abstract::Vector& result) const
00520 {
00521
00522 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00523 input.createMultiVector(1, NOX::DeepCopy);
00524 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00525 result.createMultiVector(1, NOX::DeepCopy);
00526
00527
00528 NOX::Abstract::Group::ReturnType status =
00529 applyJacobianTransposeMultiVector(*mv_input, *mv_result);
00530
00531
00532 result = (*mv_result)[0];
00533
00534 return status;
00535 }
00536
00537 NOX::Abstract::Group::ReturnType
00538 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00539 applyJacobianInverse(Teuchos::ParameterList& params,
00540 const NOX::Abstract::Vector& input,
00541 NOX::Abstract::Vector& result) const
00542 {
00543
00544 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00545 input.createMultiVector(1, NOX::DeepCopy);
00546 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00547 result.createMultiVector(1, NOX::DeepCopy);
00548
00549
00550 NOX::Abstract::Group::ReturnType status =
00551 applyJacobianInverseMultiVector(params, *mv_input, *mv_result);
00552
00553
00554 result = (*mv_result)[0];
00555
00556 return status;
00557 }
00558
00559 NOX::Abstract::Group::ReturnType
00560 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00561 applyJacobianMultiVector(const NOX::Abstract::MultiVector& input,
00562 NOX::Abstract::MultiVector& result) const
00563 {
00564 string callingFunction =
00565 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::applyJacobianMultiVector()";
00566
00567 if (!isJacobian()) {
00568 globalData->locaErrorCheck->throwError(callingFunction,
00569 "Called with invalid Jacobian!");
00570 }
00571
00572
00573 const LOCA::MultiContinuation::ExtendedMultiVector& c_input =
00574 dynamic_cast<const LOCA::MultiContinuation::ExtendedMultiVector&>(input);
00575 LOCA::MultiContinuation::ExtendedMultiVector& c_result =
00576 dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector&>(result);
00577
00578
00579 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00580 c_input.getXMultiVec();
00581 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = c_input.getScalars();
00582
00583
00584 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00585 c_result.getXMultiVec();
00586 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00587 c_result.getScalars();
00588
00589
00590 NOX::Abstract::Group::ReturnType status =
00591 borderedSolver->apply(*input_x, *input_param, *result_x, *result_param);
00592
00593 return status;
00594 }
00595
00596 NOX::Abstract::Group::ReturnType
00597 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00598 applyJacobianTransposeMultiVector(const NOX::Abstract::MultiVector& input,
00599 NOX::Abstract::MultiVector& result) const
00600 {
00601 string callingFunction =
00602 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::applyJacobianTransposeMultiVector()";
00603
00604 if (!isJacobian()) {
00605 globalData->locaErrorCheck->throwError(callingFunction,
00606 "Called with invalid Jacobian!");
00607 }
00608
00609
00610 const LOCA::MultiContinuation::ExtendedMultiVector& c_input =
00611 dynamic_cast<const LOCA::MultiContinuation::ExtendedMultiVector&>(input);
00612 LOCA::MultiContinuation::ExtendedMultiVector& c_result =
00613 dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector&>(result);
00614
00615
00616 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00617 c_input.getXMultiVec();
00618 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = c_input.getScalars();
00619
00620
00621 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00622 c_result.getXMultiVec();
00623 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00624 c_result.getScalars();
00625
00626
00627 NOX::Abstract::Group::ReturnType status =
00628 borderedSolver->applyTranspose(*input_x, *input_param, *result_x,
00629 *result_param);
00630
00631 return status;
00632 }
00633
00634 NOX::Abstract::Group::ReturnType
00635 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00636 applyJacobianInverseMultiVector(Teuchos::ParameterList& params,
00637 const NOX::Abstract::MultiVector& input,
00638 NOX::Abstract::MultiVector& result) const
00639 {
00640 string callingFunction =
00641 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::applyJacobianInverseMultiVector()";
00642
00643 if (!isJacobian()) {
00644 globalData->locaErrorCheck->throwError(callingFunction,
00645 "Called with invalid Jacobian!");
00646 }
00647
00648
00649 const LOCA::MultiContinuation::ExtendedMultiVector& c_input =
00650 dynamic_cast<const LOCA::MultiContinuation::ExtendedMultiVector&>(input);
00651 LOCA::MultiContinuation::ExtendedMultiVector& c_result =
00652 dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector&>(result);
00653
00654
00655 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00656 c_input.getXMultiVec();
00657 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = c_input.getScalars();
00658
00659
00660 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00661 c_result.getXMultiVec();
00662 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00663 c_result.getScalars();
00664
00665
00666 NOX::Abstract::Group::ReturnType status =
00667 borderedSolver->applyInverse(params, input_x.get(), input_param.get(),
00668 *result_x, *result_param);
00669
00670 return status;
00671 }
00672
00673 bool
00674 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00675 isF() const
00676 {
00677 return isValidF;
00678 }
00679
00680 bool
00681 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00682 isJacobian() const
00683 {
00684 return isValidJacobian;
00685 }
00686
00687 bool
00688 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00689 isGradient() const
00690 {
00691 return isValidGradient;
00692 }
00693
00694 bool
00695 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00696 isNewton() const
00697 {
00698 return isValidNewton;
00699 }
00700
00701 const NOX::Abstract::Vector&
00702 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00703 getX() const
00704 {
00705 return *xVec;
00706 }
00707
00708 const NOX::Abstract::Vector&
00709 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00710 getF() const
00711 {
00712 return *fVec;
00713 }
00714
00715 double
00716 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00717 getNormF() const
00718 {
00719 return fVec->norm();
00720 }
00721
00722 const NOX::Abstract::Vector&
00723 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00724 getGradient() const
00725 {
00726 return *gradientVec;
00727 }
00728
00729 const NOX::Abstract::Vector&
00730 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00731 getNewton() const
00732 {
00733 return *newtonVec;
00734 }
00735
00736 double
00737 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00738 getNormNewtonSolveResidual() const
00739 {
00740 string callingFunction =
00741 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::getNormNewtonSolveResidual()";
00742 NOX::Abstract::Group::ReturnType finalStatus;
00743 LOCA::MultiContinuation::ExtendedVector residual = *fVec;
00744
00745 finalStatus = applyJacobian(*newtonVec, residual);
00746 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00747
00748 residual = residual.update(1.0, *fVec, 1.0);
00749 return residual.norm();
00750 }
00751
00752 Teuchos::RCP<const LOCA::MultiContinuation::AbstractGroup>
00753 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00754 getUnderlyingGroup() const
00755 {
00756 return grpPtr;
00757 }
00758
00759 Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup>
00760 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00761 getUnderlyingGroup()
00762 {
00763 return grpPtr;
00764 }
00765
00766 void
00767 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00768 copy(const NOX::Abstract::Group& src)
00769 {
00770
00771 const LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup& source =
00772 dynamic_cast<const LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup&>(src);
00773
00774
00775 if (this != &source) {
00776 globalData = source.globalData;
00777 parsedParams = source.parsedParams;
00778 pitchforkParams = source.pitchforkParams;
00779 grpPtr->copy(*source.grpPtr);
00780 constraintsPtr->copy(*source.constraintsPtr);
00781 xMultiVec = source.xMultiVec;
00782 fMultiVec = source.fMultiVec;
00783 newtonMultiVec = source.newtonMultiVec;
00784 gradientMultiVec = source.gradientMultiVec;
00785 index_f = source.index_f;
00786 index_dfdp = source.index_dfdp;
00787 bifParamID = source.bifParamID;
00788 isValidF = source.isValidF;
00789 isValidJacobian = source.isValidJacobian;
00790 isValidNewton = source.isValidNewton;
00791 isValidGradient = source.isValidGradient;
00792
00793
00794 setupViews();
00795
00796
00797 borderedSolver =
00798 globalData->locaFactory->createBorderedSolverStrategy(parsedParams,
00799 pitchforkParams);
00800
00801
00802 if (isValidJacobian) {
00803 borderedSolver->setMatrixBlocks(jacOp,
00804 dfdpMultiVec->getXMultiVec(),
00805 constraintsPtr,
00806 dfdpMultiVec->getScalars());
00807 NOX::Abstract::Group::ReturnType status = borderedSolver->initForSolve();
00808 globalData->locaErrorCheck->checkReturnType(status,
00809 "LOCA::Pitchfork::MinimallyAugmented::copy()");
00810 }
00811 }
00812 }
00813
00814 void
00815 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00816 setParamsMulti(const vector<int>& paramIDs,
00817 const NOX::Abstract::MultiVector::DenseMatrix& vals)
00818 {
00819 grpPtr->setParamsMulti(paramIDs, vals);
00820 constraintsPtr->setParams(paramIDs, vals);
00821
00822 for (unsigned int i=0; i<paramIDs.size(); i++)
00823 if (paramIDs[i] == bifParamID)
00824 xVec->getScalar(0) = vals(i,0);
00825
00826 resetIsValid();
00827 }
00828
00829 void
00830 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00831 setParams(const LOCA::ParameterVector& p)
00832 {
00833 grpPtr->setParams(p);
00834 for (int i=0; i<p.length(); i++)
00835 constraintsPtr->setParam(i, p[i]);
00836 xVec->getScalar(0) = p[bifParamID];
00837
00838 resetIsValid();
00839 }
00840
00841 void
00842 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00843 setParam(int paramID, double val)
00844 {
00845 grpPtr->setParam(paramID, val);
00846 constraintsPtr->setParam(paramID, val);
00847
00848 if (paramID == bifParamID)
00849 xVec->getScalar(0) = val;
00850
00851 resetIsValid();
00852 }
00853
00854 void
00855 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00856 setParam(string paramID, double val)
00857 {
00858 const LOCA::ParameterVector& p = grpPtr->getParams();
00859 int id = p.getIndex(paramID);
00860 setParam(id, val);
00861 }
00862
00863 const LOCA::ParameterVector&
00864 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00865 getParams() const
00866 {
00867 return grpPtr->getParams();
00868 }
00869
00870 double
00871 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00872 getParam(int paramID) const
00873 {
00874 return grpPtr->getParam(paramID);
00875 }
00876
00877 double
00878 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00879 getParam(string paramID) const
00880 {
00881 return grpPtr->getParam(paramID);
00882 }
00883
00884 NOX::Abstract::Group::ReturnType
00885 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00886 computeDfDpMulti(const vector<int>& paramIDs,
00887 NOX::Abstract::MultiVector& dfdp,
00888 bool isValid_F)
00889 {
00890 string callingFunction =
00891 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::computeDfDpMulti()";
00892 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00893 NOX::Abstract::Group::ReturnType status;
00894
00895
00896 LOCA::MultiContinuation::ExtendedMultiVector& c_dfdp =
00897 dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector&>(dfdp);
00898
00899
00900 status = grpPtr->computeDfDpMulti(paramIDs, *c_dfdp.getXMultiVec(),
00901 isValid_F);
00902 finalStatus =
00903 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00904 finalStatus,
00905 callingFunction);
00906
00907
00908 status = constraintsPtr->computeDP(paramIDs,
00909 *c_dfdp.getScalars(),
00910 isValid_F);
00911 finalStatus =
00912 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00913 finalStatus,
00914 callingFunction);
00915
00916 return finalStatus;
00917 }
00918
00919 void
00920 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00921 preProcessContinuationStep(LOCA::Abstract::Iterator::StepStatus stepStatus)
00922 {
00923 grpPtr->preProcessContinuationStep(stepStatus);
00924 constraintsPtr->preProcessContinuationStep(stepStatus);
00925 }
00926
00927 void
00928 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00929 postProcessContinuationStep(LOCA::Abstract::Iterator::StepStatus stepStatus)
00930 {
00931 grpPtr->postProcessContinuationStep(stepStatus);
00932 constraintsPtr->postProcessContinuationStep(stepStatus);
00933 }
00934
00935 void
00936 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00937 projectToDraw(const NOX::Abstract::Vector& x,
00938 double *px) const
00939 {
00940 const LOCA::MultiContinuation::ExtendedVector& mx =
00941 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(x);
00942
00943 grpPtr->projectToDraw(*mx.getXVec(), px);
00944 for (int i=0; i<2; i++)
00945 px[grpPtr->projectToDrawDimension()+i] = mx.getScalar(i);
00946 }
00947
00948 int
00949 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00950 projectToDrawDimension() const
00951 {
00952 return grpPtr->projectToDrawDimension() + 2;
00953 }
00954
00955 double
00956 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00957 computeScaledDotProduct(const NOX::Abstract::Vector& a,
00958 const NOX::Abstract::Vector& b) const
00959 {
00960 const LOCA::MultiContinuation::ExtendedVector& ma =
00961 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(a);
00962 const LOCA::MultiContinuation::ExtendedVector& mb =
00963 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(b);
00964
00965 double val = grpPtr->computeScaledDotProduct(*ma.getXVec(), *mb.getXVec());
00966 for (int i=0; i<2; i++) {
00967 val += ma.getScalar(i) * mb.getScalar(i);
00968 }
00969
00970 return val;
00971 }
00972
00973 void
00974 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
00975 printSolution(const double conParam) const
00976 {
00977 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00978 globalData->locaUtils->out() <<
00979 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::printSolution\n";
00980
00981 globalData->locaUtils->out() << "Pitchfork located at: " <<
00982 globalData->locaUtils->sciformat(conParam) << " " <<
00983 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
00984
00985 globalData->locaUtils->out() <<
00986 "\tSlack variable = " <<
00987 globalData->locaUtils->sciformat(xVec->getScalar(1)) << std::endl;
00988
00989 globalData->locaUtils->out() <<
00990 "\tPrinting Solution Vector for conParam = " <<
00991 globalData->locaUtils->sciformat(conParam) << std::endl;
00992 }
00993 grpPtr->printSolution(conParam);
00994 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00995 globalData->locaUtils->out() <<
00996 "\tPrinting Right Null Vector for bif param = " <<
00997 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
00998 }
00999 grpPtr->printSolution(*(constraintsPtr->getRightNullVec()), getBifParam());
01000 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
01001 globalData->locaUtils->out() <<
01002 "\tPrinting Left Null Vector for sigma = " <<
01003 globalData->locaUtils->sciformat(constraintsPtr->getSigma()) <<
01004 std::endl;
01005 }
01006 grpPtr->printSolution(*(constraintsPtr->getLeftNullVec()),
01007 constraintsPtr->getSigma());
01008 }
01009
01010 void
01011 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01012 printSolution(const NOX::Abstract::Vector& x,
01013 const double conParam) const
01014 {
01015 const LOCA::MultiContinuation::ExtendedVector& mx =
01016 dynamic_cast<const LOCA::MultiContinuation::ExtendedVector&>(x);
01017
01018 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
01019 globalData->locaUtils->out() <<
01020 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::printSolution\n";
01021
01022 globalData->locaUtils->out() << "Pitchfork located at: " <<
01023 globalData->locaUtils->sciformat(conParam) << " " <<
01024 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
01025
01026 globalData->locaUtils->out() <<
01027 "\tSlack variable = " <<
01028 globalData->locaUtils->sciformat(mx.getScalar(1)) << std::endl;
01029
01030 globalData->locaUtils->out() <<
01031 "\tPrinting Solution Vector for conParam = " <<
01032 globalData->locaUtils->sciformat(conParam) << std::endl;
01033 }
01034 grpPtr->printSolution(*mx.getXVec(), conParam);
01035 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
01036 globalData->locaUtils->out() <<
01037 "\tPrinting Right Null Vector for bif param = " <<
01038 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
01039 }
01040 grpPtr->printSolution(*(constraintsPtr->getRightNullVec()), getBifParam());
01041 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
01042 globalData->locaUtils->out() <<
01043 "\tPrinting Left Null Vector for sigma = " <<
01044 globalData->locaUtils->sciformat(constraintsPtr->getSigma()) <<
01045 std::endl;
01046 }
01047 grpPtr->printSolution(*(constraintsPtr->getLeftNullVec()),
01048 constraintsPtr->getSigma());
01049 }
01050
01051 void
01052 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01053 scaleVector(NOX::Abstract::Vector& x) const
01054 {
01055 LOCA::MultiContinuation::ExtendedVector& mx =
01056 dynamic_cast<LOCA::MultiContinuation::ExtendedVector&>(x);
01057
01058 grpPtr->scaleVector(*mx.getXVec());
01059 }
01060
01061 int
01062 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01063 getBorderedWidth() const
01064 {
01065 int my_width = 2;
01066 if (isBordered)
01067 return my_width + bordered_grp->getBorderedWidth();
01068 else
01069 return my_width;
01070 }
01071
01072 Teuchos::RCP<const NOX::Abstract::Group>
01073 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01074 getUnborderedGroup() const
01075 {
01076 if (isBordered)
01077 return bordered_grp->getUnborderedGroup();
01078 else
01079 return grpPtr;
01080 }
01081
01082 bool
01083 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01084 isCombinedAZero() const
01085 {
01086 return false;
01087 }
01088
01089 bool
01090 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01091 isCombinedBZero() const
01092 {
01093 return false;
01094 }
01095
01096 bool
01097 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01098 isCombinedCZero() const
01099 {
01100 return false;
01101 }
01102
01103 void
01104 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01105 extractSolutionComponent(const NOX::Abstract::MultiVector& v,
01106 NOX::Abstract::MultiVector& v_x) const
01107 {
01108
01109 const LOCA::MultiContinuation::ExtendedMultiVector& mc_v =
01110 dynamic_cast<const LOCA::MultiContinuation::ExtendedMultiVector&>(v);
01111
01112
01113 Teuchos::RCP<const NOX::Abstract::MultiVector> mc_v_x =
01114 mc_v.getXMultiVec();
01115
01116
01117 if (!isBordered) {
01118 v_x = *mc_v_x;
01119 return;
01120 }
01121
01122
01123 bordered_grp->extractSolutionComponent(*mc_v_x, v_x);
01124 }
01125
01126 void
01127 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01128 extractParameterComponent(bool use_transpose,
01129 const NOX::Abstract::MultiVector& v,
01130 NOX::Abstract::MultiVector::DenseMatrix& v_p) const
01131 {
01132
01133 const LOCA::MultiContinuation::ExtendedMultiVector& mc_v =
01134 dynamic_cast<const LOCA::MultiContinuation::ExtendedMultiVector&>(v);
01135
01136
01137 Teuchos::RCP<const NOX::Abstract::MultiVector> mc_v_x =
01138 mc_v.getXMultiVec();
01139 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> mc_v_p =
01140 mc_v.getScalars();
01141
01142
01143 if (!isBordered) {
01144 if (!use_transpose)
01145 v_p.assign(*mc_v_p);
01146 else
01147 for (int j=0; j<v_p.numCols(); j++)
01148 for (int i=0; i<v_p.numRows(); i++)
01149 v_p(i,j) = (*mc_v_p)(j,i);
01150 return;
01151 }
01152
01153 int w = bordered_grp->getBorderedWidth();
01154 if (!use_transpose) {
01155
01156
01157 int num_cols = v_p.numCols();
01158 NOX::Abstract::MultiVector::DenseMatrix v_p_1(Teuchos::View, v_p,
01159 w, num_cols, 0, 0);
01160 NOX::Abstract::MultiVector::DenseMatrix v_p_2(Teuchos::View, v_p,
01161 2, num_cols, w, 0);
01162
01163
01164 bordered_grp->extractParameterComponent(use_transpose,*mc_v_x, v_p_1);
01165 v_p_2.assign(*mc_v_p);
01166 }
01167 else {
01168
01169
01170 int num_rows = v_p.numRows();
01171 NOX::Abstract::MultiVector::DenseMatrix v_p_1(Teuchos::View, v_p,
01172 num_rows, w, 0, 0);
01173 NOX::Abstract::MultiVector::DenseMatrix v_p_2(Teuchos::View, v_p,
01174 num_rows, 2, 0, w);
01175
01176
01177 bordered_grp->extractParameterComponent(use_transpose,*mc_v_x, v_p_1);
01178 for (int j=0; j<2; j++)
01179 for (int i=0; i<num_rows; i++)
01180 v_p_2(i,j) = (*mc_v_p)(j,i);
01181 }
01182 }
01183
01184 void
01185 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01186 loadNestedComponents(const NOX::Abstract::MultiVector& v_x,
01187 const NOX::Abstract::MultiVector::DenseMatrix& v_p,
01188 NOX::Abstract::MultiVector& v) const
01189 {
01190
01191 LOCA::MultiContinuation::ExtendedMultiVector& mc_v =
01192 dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector&>(v);
01193
01194
01195 Teuchos::RCP<NOX::Abstract::MultiVector> mc_v_x =
01196 mc_v.getXMultiVec();
01197 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> mc_v_p =
01198 mc_v.getScalars();
01199
01200
01201 if (!isBordered) {
01202 *mc_v_x = v_x;
01203 mc_v_p->assign(v_p);
01204 return;
01205 }
01206
01207
01208 int num_cols = v_p.numCols();
01209 int w = bordered_grp->getBorderedWidth();
01210 NOX::Abstract::MultiVector::DenseMatrix v_p_1(Teuchos::View, v_p,
01211 w, num_cols, 0, 0);
01212 NOX::Abstract::MultiVector::DenseMatrix v_p_2(Teuchos::View, v_p,
01213 2, num_cols, w, 0);
01214
01215
01216 bordered_grp->loadNestedComponents(v_x, v_p_1, *mc_v_x);
01217
01218
01219 mc_v_p->assign(v_p_2);
01220 }
01221
01222 void
01223 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01224 fillA(NOX::Abstract::MultiVector& A) const
01225 {
01226 string callingFunction =
01227 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::fillA";
01228
01229 Teuchos::RCP<const NOX::Abstract::MultiVector> my_A =
01230 dfdpMultiVec->getXMultiVec();
01231
01232
01233 if (!isBordered) {
01234 A = *my_A;
01235 return;
01236 }
01237
01238
01239 int w = bordered_grp->getBorderedWidth();
01240 std::vector<int> idx1(w);
01241 for (int i=0; i<w; i++)
01242 idx1[i] = i;
01243 Teuchos::RCP<NOX::Abstract::MultiVector> underlyingA =
01244 A.subView(idx1);
01245
01246
01247 bordered_grp->fillA(*underlyingA);
01248
01249
01250 std::vector<int> idx2(2);
01251 for (int i=0; i<2; i++)
01252 idx2[i] = w+i;
01253 Teuchos::RCP<NOX::Abstract::MultiVector> my_A_x =
01254 A.subView(idx2);
01255
01256
01257 bordered_grp->extractSolutionComponent(*my_A, *my_A_x);
01258 }
01259
01260 void
01261 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01262 fillB(NOX::Abstract::MultiVector& B) const
01263 {
01264 string callingFunction =
01265 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::fillB";
01266
01267 Teuchos::RCP<const NOX::Abstract::MultiVector> my_B =
01268 Teuchos::rcp(constraintsPtr->getDX(),false);
01269
01270
01271 if (!isBordered) {
01272 B = *my_B;
01273 return;
01274 }
01275
01276
01277 int w = bordered_grp->getBorderedWidth();
01278 std::vector<int> idx1(w);
01279 for (int i=0; i<w; i++)
01280 idx1[i] = i;
01281 Teuchos::RCP<NOX::Abstract::MultiVector> underlyingB =
01282 B.subView(idx1);
01283
01284
01285 bordered_grp->fillB(*underlyingB);
01286
01287
01288 std::vector<int> idx2(2);
01289 for (int i=0; i<2; i++)
01290 idx2[i] = w+i;
01291 Teuchos::RCP<NOX::Abstract::MultiVector> my_B_x =
01292 B.subView(idx2);
01293
01294
01295 bordered_grp->extractSolutionComponent(*my_B, *my_B_x);
01296 }
01297
01298 void
01299 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01300 fillC(NOX::Abstract::MultiVector::DenseMatrix& C) const
01301 {
01302 string callingFunction =
01303 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::fillC";
01304
01305 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> my_C =
01306 dfdpMultiVec->getScalars();
01307
01308
01309 if (!isBordered) {
01310 C.assign(*my_C);
01311 return;
01312 }
01313
01314 Teuchos::RCP<const NOX::Abstract::MultiVector> my_B =
01315 Teuchos::rcp(constraintsPtr->getDX(),false);
01316
01317 Teuchos::RCP<const NOX::Abstract::MultiVector> my_A =
01318 dfdpMultiVec->getXMultiVec();
01319
01320
01321 int w = bordered_grp->getBorderedWidth();
01322 NOX::Abstract::MultiVector::DenseMatrix underlyingC(Teuchos::View, C,
01323 w, w, 0, 0);
01324
01325
01326 bordered_grp->fillC(underlyingC);
01327
01328
01329 NOX::Abstract::MultiVector::DenseMatrix my_A_p(Teuchos::View, C,
01330 w, 2, 0, w);
01331 NOX::Abstract::MultiVector::DenseMatrix my_B_p(Teuchos::View, C,
01332 2, w, w, 0);
01333 NOX::Abstract::MultiVector::DenseMatrix my_CC(Teuchos::View, C,
01334 2, 2, w, w);
01335
01336
01337 bordered_grp->extractParameterComponent(false, *my_A, my_A_p);
01338
01339
01340 bordered_grp->extractParameterComponent(true, *my_B, my_B_p);
01341
01342
01343 my_CC.assign(*my_C);
01344 }
01345
01346 void
01347 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01348 resetIsValid()
01349 {
01350 isValidF = false;
01351 isValidJacobian = false;
01352 isValidNewton = false;
01353 isValidGradient = false;
01354 }
01355
01356 void
01357 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01358 setupViews()
01359 {
01360 index_f[0] = 0;
01361 for (int i=0; i<2; i++)
01362 index_dfdp[i] = i+1;
01363
01364 xVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedVector>(xMultiVec.getVector(0),true);
01365 fVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedVector>(fMultiVec.getVector(0),true);
01366 newtonVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedVector>(newtonMultiVec.getVector(0),true);
01367 gradientVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedVector>(gradientMultiVec.getVector(0),true);
01368
01369 ffMultiVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector>(fMultiVec.subView(index_f),true);
01370
01371 dfdpMultiVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector>(fMultiVec.subView(index_dfdp),true);
01372
01373 std::vector<int> index_fbif(2);
01374 index_fbif[0] = 0; index_fbif[1] = 1;
01375 fBifMultiVec = Teuchos::rcp_dynamic_cast<LOCA::MultiContinuation::ExtendedMultiVector>(fMultiVec.subView(index_fbif),true);
01376
01377 }
01378
01379 void
01380 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01381 setBifParam(double val)
01382 {
01383 grpPtr->setParam(bifParamID, val);
01384 xVec->getScalar(0) = val;
01385 constraintsPtr->setParam(bifParamID, val);
01386
01387 resetIsValid();
01388 }
01389
01390 void
01391 LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::
01392 getInitialVectors(Teuchos::RCP<NOX::Abstract::Vector>& aVecPtr,
01393 Teuchos::RCP<NOX::Abstract::Vector>& bVecPtr,
01394 bool isSymmetric)
01395 {
01396 string callingFunction =
01397 "LOCA::Pitchfork::MinimallyAugmented::ExtendedGroup::getIntitialVectors()";
01398
01399
01400 string method =
01401 pitchforkParams->get("Initial Null Vector Computation", "User Provided");
01402 if (method == "Solve df/dp") {
01403 NOX::Abstract::Group::ReturnType status;
01404 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
01405 std::vector<int> paramID(1);
01406 paramID[0] = bifParamID;
01407 Teuchos::RCP<NOX::Abstract::MultiVector> fdfdp =
01408 grpPtr->getX().createMultiVector(2);
01409 aVecPtr = grpPtr->getX().clone(NOX::ShapeCopy);
01410 bVecPtr = grpPtr->getX().clone(NOX::ShapeCopy);
01411 aVecPtr->init(0.0);
01412 bVecPtr->init(0.0);
01413
01414
01415 status = grpPtr->computeDfDpMulti(paramID, *fdfdp, false);
01416 finalStatus =
01417 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
01418 finalStatus,
01419 callingFunction);
01420
01421
01422 status = grpPtr->computeJacobian();
01423 finalStatus =
01424 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
01425 finalStatus,
01426 callingFunction);
01427
01428
01429 Teuchos::RCP<Teuchos::ParameterList> lsParams =
01430 parsedParams->getSublist("Linear Solver");
01431 status = grpPtr->applyJacobianInverse(*lsParams, (*fdfdp)[1], *bVecPtr);
01432 finalStatus =
01433 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
01434 finalStatus,
01435 callingFunction);
01436
01437
01438 if (!isSymmetric) {
01439
01440 Teuchos::RCP<LOCA::Abstract::TransposeSolveGroup> ts_grp =
01441 Teuchos::rcp_dynamic_cast<LOCA::Abstract::TransposeSolveGroup>(grpPtr);
01442 if (ts_grp == Teuchos::null)
01443 globalData->locaErrorCheck->throwError(
01444 callingFunction,
01445 string("Group must implement LOCA::Abstract::TransposeSolveGroup") +
01446 string(" to compute initial left null vector"));
01447
01448 Teuchos::RCP<Teuchos::ParameterList> lsParams =
01449 parsedParams->getSublist("Linear Solver");
01450 status =
01451 ts_grp->applyJacobianTransposeInverse(*lsParams, (*fdfdp)[1],
01452 *aVecPtr);
01453 finalStatus =
01454 globalData->locaErrorCheck->combineAndCheckReturnTypes(
01455 status,
01456 finalStatus,
01457 callingFunction);
01458 }
01459 else
01460 *aVecPtr = *bVecPtr;
01461
01462
01463 aVecPtr->scale(std::sqrt(static_cast<double>(aVecPtr->length())) / aVecPtr->norm());
01464 bVecPtr->scale(std::sqrt(static_cast<double>(bVecPtr->length())) / bVecPtr->norm());
01465 }
01466
01467 else {
01468
01469
01470 if (!pitchforkParams->isParameter("Initial A Vector")) {
01471 globalData->locaErrorCheck->throwError(callingFunction,
01472 "\"Initial A Vector\" is not set!");
01473 }
01474 aVecPtr =
01475 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
01476 get< Teuchos::RCP<NOX::Abstract::Vector> >("Initial A Vector");
01477
01478
01479 if (!isSymmetric) {
01480 if (!pitchforkParams->isParameter("Initial B Vector")) {
01481 globalData->locaErrorCheck->throwError(callingFunction,
01482 "\"Initial B Vector\" is not set!");
01483 }
01484 bVecPtr =
01485 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
01486 get< Teuchos::RCP<NOX::Abstract::Vector> >("Initial B Vector");
01487 }
01488 }
01489 }