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