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_BorderedSolver_Nested.H"
00043 #include "LOCA_BorderedSolver_BorderedOperator.H"
00044 #include "LOCA_BorderedSolver_JacobianOperator.H"
00045 #include "LOCA_GlobalData.H"
00046 #include "LOCA_ErrorCheck.H"
00047 #include "LOCA_Factory.H"
00048 #include "LOCA_MultiContinuation_ConstraintInterface.H"
00049 #include "LOCA_MultiContinuation_ConstraintInterfaceMVDX.H"
00050 #include "LOCA_BorderedSystem_AbstractGroup.H"
00051 #include "Teuchos_ParameterList.hpp"
00052
00053 LOCA::BorderedSolver::Nested::Nested(
00054 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00055 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00056 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams):
00057 globalData(global_data),
00058 solverParams(slvrParams),
00059 solver(),
00060 grp(),
00061 unbordered_grp(),
00062 myWidth(0),
00063 underlyingWidth(0),
00064 numConstraints(0)
00065 {
00066
00067 Teuchos::RCP<Teuchos::ParameterList> nestedSolverList =
00068 Teuchos::rcp(&(solverParams->sublist("Nested Bordered Solver")),false);
00069
00070
00071 solver =
00072 globalData->locaFactory->createBorderedSolverStrategy(topParams,
00073 nestedSolverList);
00074 }
00075
00076 LOCA::BorderedSolver::Nested::~Nested()
00077 {
00078 }
00079
00080 void
00081 LOCA::BorderedSolver::Nested::setMatrixBlocks(
00082 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& oper,
00083 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00084 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00085 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC)
00086 {
00087 string callingFunction =
00088 "LOCA::BorderedSolver::Nested::setMatrixBlocks()";
00089
00090
00091 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> op =
00092 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::JacobianOperator>(oper);
00093 if (op == Teuchos::null)
00094 globalData->locaErrorCheck->throwError(
00095 callingFunction,
00096 string("Operaror must be of type LOCA::BorderedSolver::JacobianOperator")
00097 + string(" in order to use nested bordered solver strategy."));
00098
00099
00100 grp = Teuchos::rcp_dynamic_cast<const LOCA::BorderedSystem::AbstractGroup>(op->getGroup());
00101 if (grp == Teuchos::null)
00102 globalData->locaErrorCheck->throwError(
00103 callingFunction,
00104 string("Group must be of type LOCA::BorderedSystem::AbstractGroup")
00105 + string(" in order to use nested bordered solver strategy."));
00106
00107 Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterfaceMVDX> con_mvdx = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(blockB);
00108 if (con_mvdx == Teuchos::null)
00109 globalData->locaErrorCheck->throwError(
00110 callingFunction,
00111 "Constraints object must be of type ConstraintInterfaceMVDX");
00112
00113 bool isZeroA = (blockA.get() == NULL);
00114 bool isZeroB = con_mvdx->isDXZero();
00115 bool isZeroC = (blockC.get() == NULL);
00116 Teuchos::RCP<const NOX::Abstract::MultiVector> blockB_dx;
00117 if (!isZeroB)
00118 blockB_dx = Teuchos::rcp(con_mvdx->getDX(), false);
00119
00120
00121 if (isZeroB && isZeroC)
00122 globalData->locaErrorCheck->throwError(
00123 callingFunction,
00124 "Blocks B and C cannot both be zero");
00125
00126
00127 if (isZeroA && isZeroC)
00128 globalData->locaErrorCheck->throwError(
00129 callingFunction,
00130 "Blocks A and C cannot both be zero");
00131
00132
00133 unbordered_grp = grp->getUnborderedGroup();
00134
00135
00136 if (isZeroB)
00137 numConstraints = blockC->numRows();
00138 else
00139 numConstraints = blockB_dx->numVectors();
00140
00141
00142 underlyingWidth = grp->getBorderedWidth();
00143 myWidth = underlyingWidth + numConstraints;
00144
00145
00146 bool isCombinedAZero = grp->isCombinedAZero();
00147 bool isCombinedBZero = grp->isCombinedBZero();
00148 bool isCombinedCZero = grp->isCombinedCZero();
00149 Teuchos::RCP<NOX::Abstract::MultiVector> A;
00150 Teuchos::RCP<NOX::Abstract::MultiVector> B;
00151 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> C;
00152
00153 if (!isCombinedAZero || !isZeroA) {
00154 A = unbordered_grp->getX().createMultiVector(myWidth);
00155 A->init(0.0);
00156 }
00157 if (!isCombinedBZero || !isZeroB) {
00158 B = unbordered_grp->getX().createMultiVector(myWidth);
00159 B->init(0.0);
00160
00161 }
00162 if (!isCombinedCZero || !isZeroC) {
00163 C = Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(myWidth,
00164 myWidth));
00165 C->putScalar(0.0);
00166 }
00167
00168 std::vector<int> idx1(underlyingWidth);
00169 for (int i=0; i<underlyingWidth; i++)
00170 idx1[i] = i;
00171 if (!isCombinedAZero) {
00172 Teuchos::RCP<NOX::Abstract::MultiVector> underlyingA =
00173 A->subView(idx1);
00174 grp->fillA(*underlyingA);
00175 }
00176 if (!isCombinedBZero) {
00177 Teuchos::RCP<NOX::Abstract::MultiVector> underlyingB =
00178 B->subView(idx1);
00179 grp->fillB(*underlyingB);
00180 }
00181 if (!isCombinedCZero) {
00182 NOX::Abstract::MultiVector::DenseMatrix underlyingC(Teuchos::View,
00183 *C,
00184 underlyingWidth,
00185 underlyingWidth,
00186 0, 0);
00187 grp->fillC(underlyingC);
00188 }
00189
00190 std::vector<int> idx2(numConstraints);
00191 for (int i=0; i<numConstraints; i++)
00192 idx2[i] = underlyingWidth+i;
00193 if (!isZeroA) {
00194 Teuchos::RCP<NOX::Abstract::MultiVector> my_A_x = A->subView(idx2);
00195 NOX::Abstract::MultiVector::DenseMatrix my_A_p(Teuchos::View, *C,
00196 underlyingWidth,
00197 numConstraints, 0,
00198 underlyingWidth);
00199 grp->extractSolutionComponent(*blockA, *my_A_x);
00200 grp->extractParameterComponent(false, *blockA, my_A_p);
00201 }
00202
00203 if (!isZeroB) {
00204 Teuchos::RCP<NOX::Abstract::MultiVector> my_B_x = B->subView(idx2);
00205 NOX::Abstract::MultiVector::DenseMatrix my_B_p(Teuchos::View, *C,
00206 numConstraints,
00207 underlyingWidth,
00208 underlyingWidth, 0);
00209 grp->extractSolutionComponent(*blockB_dx, *my_B_x);
00210 grp->extractParameterComponent(true, *blockB_dx, my_B_p);
00211 }
00212
00213 if (!isZeroC) {
00214 NOX::Abstract::MultiVector::DenseMatrix my_CC(Teuchos::View, *C,
00215 numConstraints,
00216 numConstraints,
00217 underlyingWidth,
00218 underlyingWidth);
00219 my_CC.assign(*blockC);
00220 }
00221
00222
00223 Teuchos::RCP<LOCA::BorderedSolver::AbstractOperator> unbordered_op =
00224 Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(unbordered_grp));
00225
00226
00227 solver->setMatrixBlocksMultiVecConstraint(unbordered_op, A, B, C);
00228 }
00229
00230 NOX::Abstract::Group::ReturnType
00231 LOCA::BorderedSolver::Nested::initForSolve()
00232 {
00233 return solver->initForSolve();
00234 }
00235
00236 NOX::Abstract::Group::ReturnType
00237 LOCA::BorderedSolver::Nested::initForTransposeSolve()
00238 {
00239 return solver->initForSolve();
00240 }
00241
00242 NOX::Abstract::Group::ReturnType
00243 LOCA::BorderedSolver::Nested::apply(
00244 const NOX::Abstract::MultiVector& X,
00245 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00246 NOX::Abstract::MultiVector& U,
00247 NOX::Abstract::MultiVector::DenseMatrix& V) const
00248 {
00249 int num_cols = X.numVectors();
00250 Teuchos::RCP<NOX::Abstract::MultiVector> XX =
00251 unbordered_grp->getX().createMultiVector(num_cols);
00252 Teuchos::RCP<NOX::Abstract::MultiVector> UU =
00253 unbordered_grp->getX().createMultiVector(num_cols);
00254 NOX::Abstract::MultiVector::DenseMatrix YY(myWidth, num_cols);
00255 NOX::Abstract::MultiVector::DenseMatrix VV(myWidth, num_cols);
00256 NOX::Abstract::MultiVector::DenseMatrix YY1(Teuchos::View, YY,
00257 underlyingWidth, num_cols,
00258 0, 0);
00259 NOX::Abstract::MultiVector::DenseMatrix YY2(Teuchos::View, YY,
00260 numConstraints, num_cols,
00261 underlyingWidth, 0);
00262 NOX::Abstract::MultiVector::DenseMatrix VV1(Teuchos::View, VV,
00263 underlyingWidth, num_cols,
00264 0, 0);
00265 NOX::Abstract::MultiVector::DenseMatrix VV2(Teuchos::View, VV,
00266 numConstraints, num_cols,
00267 underlyingWidth, 0);
00268
00269 grp->extractSolutionComponent(X, *XX);
00270 grp->extractParameterComponent(false, X, YY1);
00271 YY2.assign(Y);
00272
00273 NOX::Abstract::Group::ReturnType status =
00274 solver->apply(*XX, YY, *UU, VV);
00275
00276 V.assign(VV2);
00277 grp->loadNestedComponents(*UU, VV1, U);
00278
00279 return status;
00280 }
00281
00282 NOX::Abstract::Group::ReturnType
00283 LOCA::BorderedSolver::Nested::applyTranspose(
00284 const NOX::Abstract::MultiVector& X,
00285 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00286 NOX::Abstract::MultiVector& U,
00287 NOX::Abstract::MultiVector::DenseMatrix& V) const
00288 {
00289 int num_cols = X.numVectors();
00290 Teuchos::RCP<NOX::Abstract::MultiVector> XX =
00291 unbordered_grp->getX().createMultiVector(num_cols);
00292 Teuchos::RCP<NOX::Abstract::MultiVector> UU =
00293 unbordered_grp->getX().createMultiVector(num_cols);
00294 NOX::Abstract::MultiVector::DenseMatrix YY(myWidth, num_cols);
00295 NOX::Abstract::MultiVector::DenseMatrix VV(myWidth, num_cols);
00296 NOX::Abstract::MultiVector::DenseMatrix YY1(Teuchos::View, YY,
00297 underlyingWidth, num_cols,
00298 0, 0);
00299 NOX::Abstract::MultiVector::DenseMatrix YY2(Teuchos::View, YY,
00300 numConstraints, num_cols,
00301 underlyingWidth, 0);
00302 NOX::Abstract::MultiVector::DenseMatrix VV1(Teuchos::View, VV,
00303 underlyingWidth, num_cols,
00304 0, 0);
00305 NOX::Abstract::MultiVector::DenseMatrix VV2(Teuchos::View, VV,
00306 numConstraints, num_cols,
00307 underlyingWidth, 0);
00308
00309 grp->extractSolutionComponent(X, *XX);
00310 grp->extractParameterComponent(false, X, YY1);
00311 YY2.assign(Y);
00312
00313 NOX::Abstract::Group::ReturnType status =
00314 solver->applyTranspose(*XX, YY, *UU, VV);
00315
00316 V.assign(VV2);
00317 grp->loadNestedComponents(*UU, VV1, U);
00318
00319 return status;
00320 }
00321
00322 NOX::Abstract::Group::ReturnType
00323 LOCA::BorderedSolver::Nested::applyInverse(
00324 Teuchos::ParameterList& params,
00325 const NOX::Abstract::MultiVector* F,
00326 const NOX::Abstract::MultiVector::DenseMatrix* G,
00327 NOX::Abstract::MultiVector& X,
00328 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00329 {
00330 bool isZeroF = (F == NULL);
00331 bool isZeroG = (G == NULL);
00332
00333 if (isZeroF && isZeroG) {
00334 X.init(0.0);
00335 Y.putScalar(0.0);
00336 }
00337
00338 int num_cols = X.numVectors();
00339 Teuchos::RCP<NOX::Abstract::MultiVector> FF;
00340 if (!isZeroF)
00341 FF = unbordered_grp->getX().createMultiVector(num_cols);
00342 NOX::Abstract::MultiVector::DenseMatrix GG(myWidth, num_cols);
00343 GG.putScalar(0.0);
00344
00345 if (!isZeroF) {
00346 NOX::Abstract::MultiVector::DenseMatrix GG1(Teuchos::View, GG,
00347 underlyingWidth, num_cols,
00348 0, 0);
00349 grp->extractSolutionComponent(*F, *FF);
00350 grp->extractParameterComponent(false, *F, GG1);
00351 }
00352 if (!isZeroG) {
00353 NOX::Abstract::MultiVector::DenseMatrix GG2(Teuchos::View, GG,
00354 numConstraints, num_cols,
00355 underlyingWidth, 0);
00356 GG2.assign(*G);
00357 }
00358
00359 Teuchos::RCP<NOX::Abstract::MultiVector> XX =
00360 unbordered_grp->getX().createMultiVector(num_cols);
00361 NOX::Abstract::MultiVector::DenseMatrix YY(myWidth, num_cols);
00362 NOX::Abstract::MultiVector::DenseMatrix YY1(Teuchos::View, YY,
00363 underlyingWidth, num_cols,
00364 0, 0);
00365 NOX::Abstract::MultiVector::DenseMatrix YY2(Teuchos::View, YY,
00366 numConstraints, num_cols,
00367 underlyingWidth, 0);
00368
00369 NOX::Abstract::Group::ReturnType status =
00370 solver->applyInverse(params, FF.get(), &GG, *XX, YY);
00371
00372 Y.assign(YY2);
00373 grp->loadNestedComponents(*XX, YY1, X);
00374
00375 return status;
00376 }
00377
00378 NOX::Abstract::Group::ReturnType
00379 LOCA::BorderedSolver::Nested::applyInverseTranspose(
00380 Teuchos::ParameterList& params,
00381 const NOX::Abstract::MultiVector* F,
00382 const NOX::Abstract::MultiVector::DenseMatrix* G,
00383 NOX::Abstract::MultiVector& X,
00384 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00385 {
00386 bool isZeroF = (F == NULL);
00387 bool isZeroG = (G == NULL);
00388
00389 if (isZeroF && isZeroG) {
00390 X.init(0.0);
00391 Y.putScalar(0.0);
00392 }
00393
00394 int num_cols = X.numVectors();
00395 Teuchos::RCP<NOX::Abstract::MultiVector> FF;
00396 if (!isZeroF)
00397 FF = unbordered_grp->getX().createMultiVector(num_cols);
00398 NOX::Abstract::MultiVector::DenseMatrix GG(myWidth, num_cols);
00399 GG.putScalar(0.0);
00400
00401 if (!isZeroF) {
00402 NOX::Abstract::MultiVector::DenseMatrix GG1(Teuchos::View, GG,
00403 underlyingWidth, num_cols,
00404 0, 0);
00405 grp->extractSolutionComponent(*F, *FF);
00406 grp->extractParameterComponent(false, *F, GG1);
00407 }
00408 if (!isZeroG) {
00409 NOX::Abstract::MultiVector::DenseMatrix GG2(Teuchos::View, GG,
00410 numConstraints, num_cols,
00411 underlyingWidth, 0);
00412 GG2.assign(*G);
00413 }
00414
00415 Teuchos::RCP<NOX::Abstract::MultiVector> XX =
00416 unbordered_grp->getX().createMultiVector(num_cols);
00417 NOX::Abstract::MultiVector::DenseMatrix YY(myWidth, num_cols);
00418 NOX::Abstract::MultiVector::DenseMatrix YY1(Teuchos::View, YY,
00419 underlyingWidth, num_cols,
00420 0, 0);
00421 NOX::Abstract::MultiVector::DenseMatrix YY2(Teuchos::View, YY,
00422 numConstraints, num_cols,
00423 underlyingWidth, 0);
00424
00425 NOX::Abstract::Group::ReturnType status =
00426 solver->applyInverseTranspose(params, FF.get(), &GG, *XX, YY);
00427
00428 Y.assign(YY2);
00429 grp->loadNestedComponents(*XX, YY1, X);
00430
00431 return status;
00432 }
00433