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