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_EpetraAugmented.H"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_Vector.h"
00045 #include "NOX_Epetra_MultiVector.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_ErrorCheck.H"
00048 #include "LOCA_MultiContinuation_ConstraintInterfaceMVDX.H"
00049 #include "LOCA_Epetra_Group.H"
00050 #include "LOCA_Epetra_AugmentedOp.H"
00051 #include "LOCA_BorderedSolver_LowerTriangularBlockElimination.H"
00052 #include "LOCA_BorderedSolver_UpperTriangularBlockElimination.H"
00053 #include "LOCA_BorderedSolver_JacobianOperator.H"
00054
00055 LOCA::BorderedSolver::EpetraAugmented::EpetraAugmented(
00056 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00057 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00058 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams):
00059 globalData(global_data),
00060 solverParams(slvrParams),
00061 grp(),
00062 op(),
00063 A(),
00064 B(),
00065 C(),
00066 constraints(),
00067 numConstraints(0),
00068 isZeroA(true),
00069 isZeroB(true),
00070 isZeroC(true)
00071 {
00072 }
00073
00074 LOCA::BorderedSolver::EpetraAugmented::~EpetraAugmented()
00075 {
00076 }
00077
00078 void
00079 LOCA::BorderedSolver::EpetraAugmented::setMatrixBlocks(
00080 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& op_,
00081 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00082 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00083 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC)
00084 {
00085 string callingFunction =
00086 "LOCA::BorderedSolver::EpetraAugmented::setMatrixBlocks";
00087
00088 op = op_;
00089
00090
00091 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> jacOp =
00092 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::JacobianOperator>(op);
00093
00094 Teuchos::RCP<const NOX::Abstract::Group> group =
00095 jacOp->getGroup();
00096
00097
00098 Teuchos::RCP<NOX::Abstract::Group> non_const_group =
00099 Teuchos::rcp_const_cast<NOX::Abstract::Group>(group);
00100
00101
00102 grp = Teuchos::rcp_dynamic_cast<LOCA::Epetra::Group>(non_const_group);
00103 if (grp.get() == NULL)
00104 globalData->locaErrorCheck->throwError(
00105 callingFunction,
00106 "Group object must be an Epetra group");
00107
00108 A = blockA;
00109
00110
00111 constraints = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(blockB);
00112 if (constraints.get() == NULL)
00113 globalData->locaErrorCheck->throwError(
00114 callingFunction,
00115 "Constraints object must be of type ConstraintInterfaceMVDX");
00116 C = blockC;
00117
00118
00119 isZeroA = (A.get() == NULL);
00120 isZeroB = constraints->isDXZero();
00121 isZeroC = (C.get() == NULL);
00122
00123
00124 if (isZeroB)
00125 B = Teuchos::null;
00126 else
00127 B = Teuchos::rcp(constraints->getDX(), false);
00128
00129
00130 if (isZeroB && isZeroC)
00131 globalData->locaErrorCheck->throwError(
00132 callingFunction,
00133 "Blocks B and C cannot both be zero");
00134
00135
00136 if (isZeroA && isZeroC)
00137 globalData->locaErrorCheck->throwError(
00138 callingFunction,
00139 "Blocks A and C cannot both be zero");
00140
00141 if (isZeroB)
00142 numConstraints = C->numRows();
00143 else
00144 numConstraints = B->numVectors();
00145
00146
00147
00148
00149 if (isZeroC && !isZeroA && !isZeroB) {
00150
00151 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> tmpC =
00152 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(
00153 B->numVectors(),
00154 B->numVectors()));
00155 tmpC->putScalar(0.0);
00156 C = tmpC;
00157 isZeroC = false;
00158 }
00159
00160 }
00161
00162 NOX::Abstract::Group::ReturnType
00163 LOCA::BorderedSolver::EpetraAugmented::initForSolve()
00164 {
00165 return NOX::Abstract::Group::Ok;
00166 }
00167
00168 NOX::Abstract::Group::ReturnType
00169 LOCA::BorderedSolver::EpetraAugmented::initForTransposeSolve()
00170 {
00171 return NOX::Abstract::Group::Ok;
00172 }
00173
00174
00175 NOX::Abstract::Group::ReturnType
00176 LOCA::BorderedSolver::EpetraAugmented::apply(
00177 const NOX::Abstract::MultiVector& X,
00178 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00179 NOX::Abstract::MultiVector& U,
00180 NOX::Abstract::MultiVector::DenseMatrix& V) const
00181 {
00182
00183 NOX::Abstract::Group::ReturnType status =
00184 grp->applyJacobianMultiVector(X, U);
00185
00186
00187 if (!isZeroA)
00188 U.update(Teuchos::NO_TRANS, 1.0, *A, Y, 1.0);
00189
00190
00191 if (!isZeroB)
00192 constraints->multiplyDX(1.0, X, V);
00193
00194
00195 if (!isZeroC) {
00196 int e;
00197 if (isZeroB)
00198 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00199 else
00200 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00201 if (e < 0)
00202 status = NOX::Abstract::Group::Failed;
00203 }
00204
00205 return status;
00206 }
00207
00208 NOX::Abstract::Group::ReturnType
00209 LOCA::BorderedSolver::EpetraAugmented::applyTranspose(
00210 const NOX::Abstract::MultiVector& X,
00211 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00212 NOX::Abstract::MultiVector& U,
00213 NOX::Abstract::MultiVector::DenseMatrix& V) const
00214 {
00215
00216 NOX::Abstract::Group::ReturnType status =
00217 grp->applyJacobianTransposeMultiVector(X, U);
00218
00219
00220 if (!isZeroA)
00221 constraints->addDX(Teuchos::NO_TRANS, 1.0, Y, 1.0, U);
00222
00223
00224 if (!isZeroB)
00225 X.multiply(1.0, *A, V);
00226
00227
00228 if (!isZeroC) {
00229 int e;
00230 if (isZeroB)
00231 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00232 else
00233 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00234 if (e < 0)
00235 status = NOX::Abstract::Group::Failed;
00236 }
00237
00238 return status;
00239 }
00240
00241 NOX::Abstract::Group::ReturnType
00242 LOCA::BorderedSolver::EpetraAugmented::applyInverse(
00243 Teuchos::ParameterList& params,
00244 const NOX::Abstract::MultiVector* F,
00245 const NOX::Abstract::MultiVector::DenseMatrix* G,
00246 NOX::Abstract::MultiVector& X,
00247 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00248 {
00249 string callingFunction =
00250 "LOCA::BorderedSolver::EpetraAugmented::applyInverse()";
00251 NOX::Abstract::Group::ReturnType status;
00252 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00253
00254 bool isZeroF = (F == NULL);
00255 bool isZeroG = (G == NULL);
00256
00257
00258 if (isZeroF && isZeroG) {
00259 X.init(0.0);
00260 Y.putScalar(0.0);
00261 return NOX::Abstract::Group::Ok;
00262 }
00263
00264 if (isZeroA) {
00265 LOCA::BorderedSolver::LowerTriangularBlockElimination ltbe(globalData);
00266 return ltbe.solve(params, *op, *constraints, *C, F, G, X, Y);
00267 }
00268
00269 if (isZeroB) {
00270 LOCA::BorderedSolver::UpperTriangularBlockElimination utbe(globalData);
00271 return utbe.solve(params, *op, A.get(), *C, F, G, X, Y);
00272 }
00273
00274
00275 Teuchos::RCP<const NOX::Epetra::MultiVector> nox_epetra_a =
00276 Teuchos::rcp_dynamic_cast<const NOX::Epetra::MultiVector>(A);
00277 Teuchos::RCP<const NOX::Epetra::MultiVector> nox_epetra_b =
00278 Teuchos::rcp_dynamic_cast<const NOX::Epetra::MultiVector>(B);
00279 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_x =
00280 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(Teuchos::rcp(&X,
00281 false));
00282 Teuchos::RCP<const NOX::Epetra::MultiVector> nox_epetra_f =
00283 Teuchos::rcp_dynamic_cast<const NOX::Epetra::MultiVector>(Teuchos::rcp(F, false));
00284 Teuchos::RCP<const Epetra_MultiVector> epetra_a =
00285 Teuchos::rcp(&(nox_epetra_a->getEpetraMultiVector()), false);
00286 Teuchos::RCP<const Epetra_MultiVector> epetra_b =
00287 Teuchos::rcp(&(nox_epetra_b->getEpetraMultiVector()), false);
00288 Teuchos::RCP<Epetra_MultiVector> epetra_x =
00289 Teuchos::rcp(&(nox_epetra_x->getEpetraMultiVector()), false);
00290 Teuchos::RCP<const Epetra_MultiVector> epetra_f =
00291 Teuchos::rcp(&(nox_epetra_f->getEpetraMultiVector()), false);
00292
00293 const NOX::Epetra::Vector& solution_vec =
00294 dynamic_cast<const NOX::Epetra::Vector&>(grp->getX());
00295
00296
00297 Teuchos::RCP<NOX::Epetra::LinearSystem> linSys =
00298 grp->getLinearSystem();
00299
00300
00301 Teuchos::RCP<Epetra_Operator> jac =
00302 linSys->getJacobianOperator();
00303
00304
00305 linSys->setJacobianOperatorForSolve(jac);
00306
00307
00308 linSys->destroyPreconditioner();
00309 linSys->createPreconditioner(solution_vec, params, false);
00310
00311
00312 Teuchos::RCP<Epetra_Operator> prec =
00313 linSys->getGeneratedPrecOperator();
00314
00315
00316 LOCA::Epetra::AugmentedOp extended_jac(globalData, jac, epetra_a,
00317 epetra_b, C);
00318 LOCA::Epetra::AugmentedOp extended_prec(globalData, prec, epetra_a,
00319 epetra_b, C);
00320
00321
00322 linSys->setJacobianOperatorForSolve(Teuchos::rcp(&extended_jac,false));
00323 linSys->setPrecOperatorForSolve(Teuchos::rcp(&extended_prec,false));
00324
00325
00326 Teuchos::RCP<Epetra_MultiVector> epetra_augmented_x =
00327 extended_jac.buildEpetraAugmentedMultiVec(*epetra_x, &Y, false);
00328 Teuchos::RCP<Epetra_MultiVector> epetra_augmented_f =
00329 extended_jac.buildEpetraAugmentedMultiVec(*epetra_f, G, true);
00330
00331
00332 NOX::Epetra::MultiVector nox_epetra_augmented_x(
00333 epetra_augmented_x,
00334 NOX::DeepCopy,
00335 NOX::Epetra::MultiVector::CreateView);
00336 NOX::Epetra::MultiVector nox_epetra_augmented_f(
00337 epetra_augmented_f,
00338 NOX::DeepCopy,
00339 NOX::Epetra::MultiVector::CreateView);
00340
00341
00342 int m = nox_epetra_augmented_x.numVectors();
00343 for (int i=0; i<m; i++) {
00344 extended_jac.init(*((*epetra_f)(i)));
00345 extended_prec.init(*((*epetra_f)(i)));
00346 bool stat =
00347 linSys->applyJacobianInverse(
00348 params,
00349 dynamic_cast<NOX::Epetra::Vector&>(nox_epetra_augmented_f[i]),
00350 dynamic_cast<NOX::Epetra::Vector&>(nox_epetra_augmented_x[i]));
00351 if (stat == true)
00352 status = NOX::Abstract::Group::Ok;
00353 else
00354 status = NOX::Abstract::Group::NotConverged;
00355 finalStatus =
00356 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00357 finalStatus,
00358 callingFunction);
00359 }
00360
00361
00362 extended_jac.init(*epetra_x);
00363 extended_jac.setEpetraAugmentedMultiVec(*epetra_x, Y,
00364 *epetra_augmented_x);
00365
00366
00367 linSys->setJacobianOperatorForSolve(jac);
00368 linSys->destroyPreconditioner();
00369
00370
00371 return finalStatus;
00372 }