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_EpetraHouseholder.H"
00043 #include "Epetra_Vector.h"
00044 #include "Epetra_MultiVector.h"
00045 #include "Epetra_RowMatrix.h"
00046 #include "Epetra_CrsMatrix.h"
00047 #include "NOX_Epetra_MultiVector.H"
00048 #include "LOCA_GlobalData.H"
00049 #include "LOCA_ErrorCheck.H"
00050 #include "LOCA_MultiContinuation_ConstraintInterfaceMVDX.H"
00051 #include "LOCA_Epetra_Group.H"
00052 #include "LOCA_Epetra_CompactWYOp.H"
00053 #include "LOCA_Epetra_LowRankUpdateOp.H"
00054 #include "LOCA_Epetra_LowRankUpdateRowMatrix.H"
00055 #include "LOCA_Epetra_TransposeLinearSystem_AbstractStrategy.H"
00056 #include "LOCA_Epetra_TransposeLinearSystem_Factory.H"
00057 #include "Teuchos_ParameterList.hpp"
00058 #include "LOCA_BorderedSolver_LowerTriangularBlockElimination.H"
00059 #include "LOCA_BorderedSolver_UpperTriangularBlockElimination.H"
00060 #include "LOCA_Abstract_TransposeSolveGroup.H"
00061 #include "LOCA_BorderedSolver_JacobianOperator.H"
00062 #include "LOCA_BorderedSolver_ComplexOperator.H"
00063 #include "LOCA_Hopf_ComplexMultiVector.H"
00064
00065 #ifdef HAVE_NOX_EPETRAEXT
00066 #include "EpetraExt_BlockCrsMatrix.h"
00067 #include "EpetraExt_BlockVector.h"
00068 #include "EpetraExt_BlockMultiVector.h"
00069 #endif
00070
00071 LOCA::BorderedSolver::EpetraHouseholder::EpetraHouseholder(
00072 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00073 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00074 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams):
00075 globalData(global_data),
00076 solverParams(slvrParams),
00077 grp(),
00078 op(),
00079 A(),
00080 B(),
00081 C(),
00082 constraints(),
00083 qrFact(globalData),
00084 house_x(),
00085 house_p(),
00086 T(),
00087 R(),
00088 U(),
00089 V(),
00090 house_x_trans(),
00091 house_p_trans(),
00092 T_trans(),
00093 R_trans(),
00094 U_trans(),
00095 V_trans(),
00096 Ablock(),
00097 Bblock(),
00098 Ascaled(),
00099 Bscaled(),
00100 Cscaled(),
00101 linSys(),
00102 epetraOp(),
00103 baseMap(),
00104 globalMap(),
00105 numConstraints(0),
00106 isZeroA(true),
00107 isZeroB(true),
00108 isZeroC(true),
00109 isValidForSolve(false),
00110 isValidForTransposeSolve(false),
00111 dblas(),
00112 scale_rows(true),
00113 scale_vals(),
00114 precMethod(JACOBIAN),
00115 includeUV(false),
00116 use_P_For_Prec(false),
00117 isComplex(false),
00118 omega(0.0)
00119 {
00120 scale_rows = solverParams->get("Scale Augmented Rows", true);
00121 std::string prec_method =
00122 solverParams->get("Preconditioner Method", "Jacobian");
00123 if (prec_method == "Jacobian")
00124 precMethod = JACOBIAN;
00125 else if (prec_method == "SMW")
00126 precMethod = SMW;
00127 else
00128 globalData->locaErrorCheck->throwError(
00129 "LOCA::BorderedSolver::EpetraHouseholder::EpetraHouseholder()",
00130 "Unknown preconditioner method! Choices are Jacobian, SMW");
00131 includeUV =
00132 solverParams->get("Include UV In Preconditioner", false);
00133 use_P_For_Prec =
00134 solverParams->get("Use P For Preconditioner", false);
00135 }
00136
00137 LOCA::BorderedSolver::EpetraHouseholder::~EpetraHouseholder()
00138 {
00139 }
00140
00141 void
00142 LOCA::BorderedSolver::EpetraHouseholder::setMatrixBlocks(
00143 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& op_,
00144 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00145 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00146 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC)
00147 {
00148 string callingFunction =
00149 "LOCA::BorderedSolver::EpetraHouseholder::setMatrixBlocks";
00150
00151 op = op_;
00152 A = blockA;
00153
00154
00155 constraints = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(blockB);
00156 if (constraints.get() == NULL)
00157 globalData->locaErrorCheck->throwError(
00158 callingFunction,
00159 "Constraints object must be of type ConstraintInterfaceMVDX");
00160 C = blockC;
00161
00162
00163 isZeroA = (A.get() == NULL);
00164 isZeroB = constraints->isDXZero();
00165 isZeroC = (C.get() == NULL);
00166
00167
00168 if (isZeroB)
00169 B = Teuchos::null;
00170 else
00171 B = Teuchos::rcp(constraints->getDX(), false);
00172
00173
00174 if (isZeroB && isZeroC)
00175 globalData->locaErrorCheck->throwError(
00176 callingFunction,
00177 "Blocks B and C cannot both be zero");
00178
00179
00180 if (isZeroA && isZeroC)
00181 globalData->locaErrorCheck->throwError(
00182 callingFunction,
00183 "Blocks A and C cannot both be zero");
00184
00185 if (isZeroB)
00186 numConstraints = C->numRows();
00187 else
00188 numConstraints = B->numVectors();
00189
00190
00191
00192
00193 if (isZeroC && !isZeroA && !isZeroB) {
00194
00195 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> tmpC =
00196 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(
00197 B->numVectors(),
00198 B->numVectors()));
00199 tmpC->putScalar(0.0);
00200 C = tmpC;
00201 isZeroC = false;
00202 }
00203
00204
00205 isValidForSolve = false;
00206 isValidForTransposeSolve = false;
00207
00208 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> jacOp =
00209 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::JacobianOperator>(op);
00210 Teuchos::RCP<const LOCA::BorderedSolver::ComplexOperator> complexOp =
00211 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::ComplexOperator>(op);
00212
00213 if (jacOp != Teuchos::null) {
00214
00215 isComplex = false;
00216
00217
00218 Teuchos::RCP<const LOCA::Epetra::Group> constGrp =
00219 Teuchos::rcp_dynamic_cast<const LOCA::Epetra::Group>(jacOp->getGroup());
00220 if (constGrp.get() == NULL)
00221 globalData->locaErrorCheck->throwError(
00222 callingFunction,
00223 "Group object must be an Epetra group");
00224
00225
00226 Ablock = A;
00227
00228
00229 Bblock = B;
00230
00231
00232 grp = Teuchos::rcp_const_cast<LOCA::Epetra::Group>(constGrp);
00233
00234
00235 linSys = grp->getLinearSystem();
00236 epetraOp = linSys->getJacobianOperator();
00237 }
00238 else if (complexOp != Teuchos::null) {
00239
00240 isComplex = true;
00241
00242 omega = complexOp->getFrequency();
00243
00244
00245 Teuchos::RCP<const LOCA::Epetra::Group> constGrp =
00246 Teuchos::rcp_dynamic_cast<const LOCA::Epetra::Group>(complexOp->getGroup());
00247 if (constGrp.get() == NULL)
00248 globalData->locaErrorCheck->throwError(
00249 callingFunction,
00250 "Group object must be an Epetra group");
00251
00252
00253 grp = Teuchos::rcp_const_cast<LOCA::Epetra::Group>(constGrp);
00254
00255
00256 linSys = grp->getComplexLinearSystem();
00257 epetraOp = linSys->getJacobianOperator();
00258
00259
00260 grp->getComplexMaps(baseMap, globalMap);
00261
00262
00263 if (!isZeroA)
00264 Ablock = createBlockMV(*A);
00265
00266
00267 if (!isZeroB)
00268 Bblock = createBlockMV(*B);
00269 }
00270 else {
00271 globalData->locaErrorCheck->throwError(
00272 callingFunction,
00273 string("Op argument must be of type !\n") +
00274 string("LOCA::BorderedSolver::JacobianOperator or \n") +
00275 string("LOCA::BorderedSolver::ComplexOperator."));
00276 }
00277
00278 Ascaled = Teuchos::null;
00279 Bscaled = Teuchos::null;
00280 Cscaled = Teuchos::null;
00281
00282 if (Ablock != Teuchos::null)
00283 Ascaled = Ablock->clone();
00284 if (Bblock != Teuchos::null)
00285 Bscaled = Bblock->clone();
00286 if (C != Teuchos::null) {
00287 Cscaled =
00288 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(
00289 C->numRows(),
00290 C->numCols()));
00291 Cscaled->assign(*C);
00292 }
00293
00294 }
00295
00296 NOX::Abstract::Group::ReturnType
00297 LOCA::BorderedSolver::EpetraHouseholder::initForSolve()
00298 {
00299 NOX::Abstract::Group::ReturnType res = NOX::Abstract::Group::Ok;
00300
00301
00302 if (isValidForSolve)
00303 return res;
00304
00305
00306 if (!isZeroA && !isZeroB) {
00307
00308
00309 if (scale_rows) {
00310 scale_vals.resize(numConstraints);
00311 double sn = std::sqrt( static_cast<double>(Bblock->length()) );
00312 for (int i=0; i<numConstraints; i++) {
00313 scale_vals[i] = (*Bblock)[i].norm() / sn;
00314 double t = 0.0;
00315 for (int j=0; j<numConstraints; j++)
00316 t += (*C)(i,j) * (*C)(i,j);
00317 scale_vals[i] += std::sqrt(t);
00318 scale_vals[i] = 1.0 / scale_vals[i];
00319 (*Bscaled)[i].scale(scale_vals[i]);
00320 for (int j=0; j<numConstraints; j++)
00321 (*Cscaled)(i,j) *= scale_vals[i];
00322 }
00323 }
00324
00325
00326 if (house_x.get() == NULL || house_x->numVectors() != numConstraints) {
00327 house_x = Bblock->clone(NOX::ShapeCopy);
00328 U = house_x->clone(NOX::ShapeCopy);
00329 V = house_x->clone(NOX::ShapeCopy);
00330 house_p.reshape(numConstraints, numConstraints);
00331 T.reshape(numConstraints, numConstraints);
00332 R.reshape(numConstraints, numConstraints);
00333 }
00334
00335
00336 qrFact.computeQR(*Cscaled, *Bscaled, true, house_p, *house_x, T, R);
00337
00338
00339 res = computeUV(house_p, *house_x, T, *Ablock, *U, *V, false);
00340 globalData->locaErrorCheck->checkReturnType(res,
00341 "LOCA::BorderedSolver::Epetra_Householder::initForSolve()");
00342 }
00343
00344 isValidForSolve = true;
00345
00346 return res;
00347 }
00348
00349 NOX::Abstract::Group::ReturnType
00350 LOCA::BorderedSolver::EpetraHouseholder::initForTransposeSolve()
00351 {
00352 NOX::Abstract::Group::ReturnType res = NOX::Abstract::Group::Ok;
00353
00354
00355 if (isValidForTransposeSolve)
00356 return res;
00357
00358
00359 if (!isZeroA && !isZeroB) {
00360
00361
00362 if (scale_rows) {
00363 scale_vals.resize(numConstraints);
00364 double sn = std::sqrt( static_cast<double>(Ablock->length()) );
00365 for (int i=0; i<numConstraints; i++) {
00366 scale_vals[i] = (*Ablock)[i].norm() / sn;
00367 double t = 0.0;
00368 for (int j=0; j<numConstraints; j++)
00369 t += (*C)(j,i) * (*C)(j,i);
00370 scale_vals[i] += std::sqrt(t);
00371 scale_vals[i] = 1.0 / scale_vals[i];
00372 (*Ascaled)[i].scale(scale_vals[i]);
00373 for (int j=0; j<numConstraints; j++)
00374 (*Cscaled)(j,i) *= scale_vals[i];
00375 }
00376 }
00377
00378
00379 if (house_x_trans.get() == NULL ||
00380 house_x_trans->numVectors() != numConstraints) {
00381 house_x_trans = Ablock->clone(NOX::ShapeCopy);
00382 U_trans = house_x_trans->clone(NOX::ShapeCopy);
00383 V_trans = house_x_trans->clone(NOX::ShapeCopy);
00384 house_p_trans.reshape(numConstraints, numConstraints);
00385 T_trans.reshape(numConstraints, numConstraints);
00386 R_trans.reshape(numConstraints, numConstraints);
00387 }
00388
00389
00390 qrFact.computeQR(*Cscaled, *Ascaled, false, house_p_trans, *house_x_trans,
00391 T_trans, R_trans);
00392
00393
00394 res = computeUV(house_p_trans, *house_x_trans, T_trans, *Bblock, *U_trans,
00395 *V_trans, true);
00396 globalData->locaErrorCheck->checkReturnType(res,
00397 "LOCA::BorderedSolver::Epetra_Householder::initForTransposeSolve()");
00398 }
00399
00400 isValidForTransposeSolve = true;
00401
00402 return res;
00403 }
00404
00405 NOX::Abstract::Group::ReturnType
00406 LOCA::BorderedSolver::EpetraHouseholder::apply(
00407 const NOX::Abstract::MultiVector& X,
00408 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00409 NOX::Abstract::MultiVector& U,
00410 NOX::Abstract::MultiVector::DenseMatrix& V) const
00411 {
00412
00413 NOX::Abstract::Group::ReturnType status = op->apply(X, U);
00414
00415
00416 if (!isZeroA)
00417 U.update(Teuchos::NO_TRANS, 1.0, *A, Y, 1.0);
00418
00419
00420 if (!isZeroB)
00421 constraints->multiplyDX(1.0, X, V);
00422
00423
00424 if (!isZeroC) {
00425 int e;
00426 if (isZeroB)
00427 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00428 else
00429 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00430 if (e < 0)
00431 status = NOX::Abstract::Group::Failed;
00432 }
00433
00434 return status;
00435 }
00436
00437 NOX::Abstract::Group::ReturnType
00438 LOCA::BorderedSolver::EpetraHouseholder::applyTranspose(
00439 const NOX::Abstract::MultiVector& X,
00440 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00441 NOX::Abstract::MultiVector& U,
00442 NOX::Abstract::MultiVector::DenseMatrix& V) const
00443 {
00444
00445 NOX::Abstract::Group::ReturnType status = op->applyTranspose(X, U);
00446
00447
00448 if (!isZeroA)
00449 constraints->addDX(Teuchos::NO_TRANS, 1.0, Y, 1.0, U);
00450
00451
00452 if (!isZeroB)
00453 X.multiply(1.0, *A, V);
00454
00455
00456 if (!isZeroC) {
00457 int e;
00458 if (isZeroB)
00459 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00460 else
00461 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00462 if (e < 0)
00463 status = NOX::Abstract::Group::Failed;
00464 }
00465
00466 return status;
00467 }
00468
00469 NOX::Abstract::Group::ReturnType
00470 LOCA::BorderedSolver::EpetraHouseholder::applyInverse(
00471 Teuchos::ParameterList& params,
00472 const NOX::Abstract::MultiVector* F,
00473 const NOX::Abstract::MultiVector::DenseMatrix* G,
00474 NOX::Abstract::MultiVector& X,
00475 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00476 {
00477 bool isZeroF = (F == NULL);
00478 bool isZeroG = (G == NULL);
00479
00480
00481 if (isZeroF && isZeroG) {
00482 X.init(0.0);
00483 Y.putScalar(0.0);
00484 return NOX::Abstract::Group::Ok;
00485 }
00486
00487 if (isZeroA) {
00488 LOCA::BorderedSolver::LowerTriangularBlockElimination ltbe(globalData);
00489 return ltbe.solve(params, *op, *constraints, *C, F, G, X, Y);
00490 }
00491
00492 if (isZeroB) {
00493 LOCA::BorderedSolver::UpperTriangularBlockElimination utbe(globalData);
00494 return utbe.solve(params, *op, A.get(), *C, F, G, X, Y);
00495 }
00496
00497
00498 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> Gscaled;
00499 if (G != NULL) {
00500 Gscaled =
00501 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(G->numRows(),
00502 G->numCols()));
00503 if (scale_rows)
00504 for (int i=0; i<G->numRows(); i++)
00505 for (int j=0; j<G->numCols(); j++)
00506 (*Gscaled)(i,j) = scale_vals[i]*(*G)(i,j);
00507 else
00508 Gscaled->assign(*G);
00509 }
00510
00511 NOX::Abstract::Group::ReturnType res;
00512 if (!isComplex)
00513 res = solve(params, F, Gscaled.get(), X, Y);
00514 else {
00515 Teuchos::RCP<NOX::Abstract::MultiVector> blockF;
00516 if (!isZeroF)
00517 blockF = createBlockMV(*F);
00518 Teuchos::RCP<NOX::Abstract::MultiVector> blockX = createBlockMV(X);
00519 res = solve(params, blockF.get(), Gscaled.get(), *blockX, Y);
00520 setBlockMV(*blockX, X);
00521 }
00522
00523 return res;
00524 }
00525
00526 NOX::Abstract::Group::ReturnType
00527 LOCA::BorderedSolver::EpetraHouseholder::applyInverseTranspose(
00528 Teuchos::ParameterList& params,
00529 const NOX::Abstract::MultiVector* F,
00530 const NOX::Abstract::MultiVector::DenseMatrix* G,
00531 NOX::Abstract::MultiVector& X,
00532 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00533 {
00534 bool isZeroF = (F == NULL);
00535 bool isZeroG = (G == NULL);
00536
00537
00538 if (isZeroF && isZeroG) {
00539 X.init(0.0);
00540 Y.putScalar(0.0);
00541 return NOX::Abstract::Group::Ok;
00542 }
00543
00544
00545 Teuchos::RCP<const LOCA::Abstract::TransposeSolveGroup> ts_grp;
00546 if (isZeroA || isZeroB) {
00547 ts_grp =
00548 Teuchos::rcp_dynamic_cast<const LOCA::Abstract::TransposeSolveGroup>(grp);
00549 if (ts_grp == Teuchos::null)
00550 globalData->locaErrorCheck->throwError(
00551 "LOCA::BorderedSolver::EpetraHouseholder::applyInverseTranspose()",
00552 string("Group must implement the LOCA::Abstract::TransposeSolveGroup")
00553 + string(" interface in order to solve the transpose of the bordered")
00554 + string(" system via bordering."));
00555 }
00556
00557 if (isZeroA) {
00558 LOCA::BorderedSolver::UpperTriangularBlockElimination utbe(globalData);
00559 return utbe.solveTranspose(params, *op, B.get(), *C, F, G, X, Y);
00560 }
00561
00562 else if (isZeroB) {
00563 LOCA::BorderedSolver::LowerTriangularBlockElimination ltbe(globalData);
00564 return ltbe.solveTranspose(params, *op, *A, *C, F, G, X, Y);
00565 }
00566
00567
00568 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> Gscaled;
00569 if (G != NULL) {
00570 Gscaled =
00571 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(G->numRows(),
00572 G->numCols()));
00573 if (scale_rows)
00574 for (int i=0; i<G->numRows(); i++)
00575 for (int j=0; j<G->numCols(); j++)
00576 (*Gscaled)(i,j) = scale_vals[i]*(*G)(i,j);
00577 else
00578 Gscaled->assign(*G);
00579 }
00580
00581 NOX::Abstract::Group::ReturnType res;
00582 if (!isComplex)
00583 res = solveTranspose(params, F, Gscaled.get(), X, Y);
00584 else {
00585 Teuchos::RCP<NOX::Abstract::MultiVector> blockF;
00586 if (!isZeroF)
00587 blockF = createBlockMV(*F);
00588 Teuchos::RCP<NOX::Abstract::MultiVector> blockX = createBlockMV(X);
00589 res = solveTranspose(params, blockF.get(), Gscaled.get(), *blockX, Y);
00590 setBlockMV(*blockX, X);
00591 }
00592
00593 return res;
00594 }
00595
00596 NOX::Abstract::Group::ReturnType
00597 LOCA::BorderedSolver::EpetraHouseholder::solve(
00598 Teuchos::ParameterList& params,
00599 const NOX::Abstract::MultiVector* F,
00600 const NOX::Abstract::MultiVector::DenseMatrix* G,
00601 NOX::Abstract::MultiVector& X,
00602 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00603 {
00604 string callingFunction =
00605 "LOCA::BorderedSolver::EpetraHouseholder::applyInverse()";
00606 NOX::Abstract::Group::ReturnType status;
00607 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00608
00609 bool isZeroF = (F == NULL);
00610 bool isZeroG = (G == NULL);
00611
00612
00613 if (!isValidForSolve)
00614 globalData->locaErrorCheck->throwError(
00615 callingFunction,
00616 string("applyInverse() called with invalid constraint") +
00617 string(" factorizations. Call initForSolve() first."));
00618
00619 Teuchos::RCP<const NOX::Abstract::MultiVector> cRHS;
00620
00621 if (!isZeroG) {
00622
00623 Teuchos::RCP<NOX::Epetra::MultiVector> tmp_x =
00624 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(X.clone(NOX::ShapeCopy));
00625 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> tmp_y =
00626 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(G->numRows(),
00627 G->numCols()));
00628 Teuchos::RCP<NOX::Epetra::MultiVector> RHS =
00629 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(X.clone(NOX::ShapeCopy));
00630
00631
00632 Y.assign(*G);
00633 dblas.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::TRANS,
00634 Teuchos::NON_UNIT_DIAG,
00635 G->numRows(), G->numCols(), 1.0, R.values(),
00636 R.numRows(), Y.values(), Y.numRows());
00637
00638
00639 qrFact.applyCompactWY(house_p, *house_x, T, &Y, NULL, *tmp_y, *tmp_x,
00640 false);
00641
00642
00643 int res = epetraOp->Apply(tmp_x->getEpetraMultiVector(),
00644 RHS->getEpetraMultiVector());
00645 if (res == 0)
00646 status = NOX::Abstract::Group::Ok;
00647 else
00648 status = NOX::Abstract::Group::Failed;
00649 finalStatus =
00650 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00651 finalStatus,
00652 callingFunction);
00653 RHS->update(Teuchos::NO_TRANS, -1.0, *Ablock, *tmp_y, -1.0);
00654
00655
00656 if (!isZeroF)
00657 RHS->update(1.0, *F, 1.0);
00658
00659 cRHS = RHS;
00660 }
00661 else
00662 cRHS = Teuchos::rcp(F, false);
00663
00664
00665 const NOX::Epetra::Vector& solution_vec =
00666 dynamic_cast<const NOX::Epetra::Vector&>(grp->getX());
00667
00668
00669 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_U =
00670 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(U);
00671 Teuchos::RCP<Epetra_MultiVector> epetra_U =
00672 Teuchos::rcp(&(nox_epetra_U->getEpetraMultiVector()), false);
00673 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_V =
00674 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(V);
00675 Teuchos::RCP<Epetra_MultiVector> epetra_V =
00676 Teuchos::rcp(&(nox_epetra_V->getEpetraMultiVector()), false);
00677
00678
00679 Teuchos::RCP<Epetra_RowMatrix> jac_rowmatrix =
00680 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(epetraOp);
00681 Teuchos::RCP<Epetra_Operator> op;
00682 if (jac_rowmatrix != Teuchos::null)
00683 op = Teuchos::rcp(new LOCA::Epetra::LowRankUpdateRowMatrix(globalData,
00684 jac_rowmatrix,
00685 epetra_U,
00686 epetra_V,
00687 false,
00688 includeUV));
00689 else
00690 op = Teuchos::rcp(new LOCA::Epetra::LowRankUpdateOp(globalData,
00691 epetraOp,
00692 epetra_U,
00693 epetra_V,
00694 false));
00695
00696
00697
00698 Teuchos::RCP<Epetra_CrsMatrix> jac_crs;
00699 if (precMethod == JACOBIAN && includeUV && !use_P_For_Prec) {
00700 jac_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(epetraOp);
00701 if (jac_crs != Teuchos::null) {
00702 updateJacobianForPreconditioner(*U, *V, *jac_crs);
00703 }
00704 }
00705
00706
00707 if (precMethod == JACOBIAN && use_P_For_Prec)
00708 linSys->setJacobianOperatorForSolve(op);
00709 else
00710 linSys->setJacobianOperatorForSolve(epetraOp);
00711
00712
00713 linSys->destroyPreconditioner();
00714 linSys->createPreconditioner(solution_vec, params, false);
00715
00716
00717 if (precMethod == JACOBIAN && includeUV && !use_P_For_Prec &&
00718 jac_crs != Teuchos::null) {
00719 grp->setX(solution_vec);
00720 if (isComplex)
00721 grp->computeComplex(omega);
00722 else
00723 grp->computeJacobian();
00724 }
00725
00726
00727 linSys->setJacobianOperatorForSolve(op);
00728
00729
00730 Teuchos::RCP<Epetra_Operator> prec_op;
00731 Teuchos::RCP<Epetra_Operator> epetraPrecOp;
00732 if (precMethod == SMW) {
00733 epetraPrecOp = linSys->getGeneratedPrecOperator();
00734 prec_op = Teuchos::rcp(new LOCA::Epetra::LowRankUpdateOp(globalData,
00735 epetraPrecOp,
00736 epetra_U,
00737 epetra_V,
00738 true));
00739 linSys->setPrecOperatorForSolve(prec_op);
00740 }
00741
00742
00743 int m = X.numVectors();
00744 X.init(0.0);
00745 for (int i=0; i<m; i++) {
00746 bool stat =
00747 linSys->applyJacobianInverse(
00748 params,
00749 dynamic_cast<const NOX::Epetra::Vector&>((*cRHS)[i]),
00750 dynamic_cast<NOX::Epetra::Vector&>(X[i]));
00751
00752 if (stat == true)
00753 status = NOX::Abstract::Group::Ok;
00754 else
00755 status = NOX::Abstract::Group::NotConverged;
00756 finalStatus =
00757 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00758 finalStatus,
00759 callingFunction);
00760 }
00761
00762 qrFact.applyCompactWY(house_p, *house_x, T, Y, X, isZeroG, false, false);
00763
00764
00765 linSys->setJacobianOperatorForSolve(epetraOp);
00766 if (precMethod == SMW)
00767 linSys->setPrecOperatorForSolve(epetraPrecOp);
00768 linSys->destroyPreconditioner();
00769
00770 return finalStatus;
00771 }
00772
00773 NOX::Abstract::Group::ReturnType
00774 LOCA::BorderedSolver::EpetraHouseholder::solveTranspose(
00775 Teuchos::ParameterList& params,
00776 const NOX::Abstract::MultiVector* F,
00777 const NOX::Abstract::MultiVector::DenseMatrix* G,
00778 NOX::Abstract::MultiVector& X,
00779 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00780 {
00781 string callingFunction =
00782 "LOCA::BorderedSolver::EpetraHouseholder::applyInverseTranspose()";
00783 NOX::Abstract::Group::ReturnType status;
00784 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00785
00786 bool isZeroF = (F == NULL);
00787 bool isZeroG = (G == NULL);
00788
00789
00790 if (!isValidForTransposeSolve)
00791 globalData->locaErrorCheck->throwError(
00792 callingFunction,
00793 string("applyInverseTranspose() called with invalid constraint") +
00794 string(" factorizations. Call initForSolve() first."));
00795
00796 Teuchos::RCP<const NOX::Abstract::MultiVector> cRHS;
00797
00798 if (!isZeroG) {
00799
00800 Teuchos::RCP<NOX::Epetra::MultiVector> tmp_x =
00801 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(X.clone(NOX::ShapeCopy));
00802 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> tmp_y =
00803 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(G->numRows(),
00804 G->numCols()));
00805 Teuchos::RCP<NOX::Epetra::MultiVector> RHS =
00806 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(X.clone(NOX::ShapeCopy));
00807
00808
00809 Y.assign(*G);
00810 dblas.TRSM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::TRANS,
00811 Teuchos::NON_UNIT_DIAG,
00812 G->numRows(), G->numCols(), 1.0, R_trans.values(),
00813 R_trans.numRows(), Y.values(), Y.numRows());
00814
00815
00816 qrFact.applyCompactWY(house_p_trans, *house_x_trans, T_trans, &Y, NULL,
00817 *tmp_y, *tmp_x, false);
00818
00819
00820 epetraOp->SetUseTranspose(true);
00821 int res = epetraOp->Apply(tmp_x->getEpetraMultiVector(),
00822 RHS->getEpetraMultiVector());
00823 epetraOp->SetUseTranspose(false);
00824 if (res == 0)
00825 status = NOX::Abstract::Group::Ok;
00826 else
00827 status = NOX::Abstract::Group::Failed;
00828 finalStatus =
00829 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00830 finalStatus,
00831 callingFunction);
00832 RHS->update(Teuchos::NO_TRANS, -1.0, *Bblock, *tmp_y, -1.0);
00833
00834
00835 if (!isZeroF)
00836 RHS->update(1.0, *F, 1.0);
00837
00838 cRHS = RHS;
00839 }
00840 else
00841 cRHS = Teuchos::rcp(F, false);
00842
00843
00844 const NOX::Epetra::Vector& solution_vec =
00845 dynamic_cast<const NOX::Epetra::Vector&>(grp->getX());
00846
00847
00848 LOCA::Epetra::TransposeLinearSystem::Factory tls_factory(globalData);
00849
00850 Teuchos::RCP<LOCA::Epetra::TransposeLinearSystem::AbstractStrategy> tls_strategy = tls_factory.create(solverParams, linSys);
00851
00852
00853 tls_strategy->createJacobianTranspose();
00854 Teuchos::RCP<Epetra_Operator> jac_trans =
00855 tls_strategy->getJacobianTransposeOperator();
00856
00857
00858 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_U =
00859 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(U_trans);
00860 Teuchos::RCP<Epetra_MultiVector> epetra_U =
00861 Teuchos::rcp(&(nox_epetra_U->getEpetraMultiVector()), false);
00862 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_V =
00863 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(V_trans);
00864 Teuchos::RCP<Epetra_MultiVector> epetra_V =
00865 Teuchos::rcp(&(nox_epetra_V->getEpetraMultiVector()), false);
00866
00867
00868 Teuchos::RCP<Epetra_RowMatrix> jac_trans_rowmatrix =
00869 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(jac_trans);
00870 Teuchos::RCP<Epetra_Operator> op;
00871 if (jac_trans_rowmatrix != Teuchos::null)
00872 op = Teuchos::rcp(new LOCA::Epetra::LowRankUpdateRowMatrix(
00873 globalData,
00874 jac_trans_rowmatrix,
00875 epetra_U,
00876 epetra_V,
00877 false,
00878 includeUV));
00879 else
00880 op = Teuchos::rcp(new LOCA::Epetra::LowRankUpdateOp(globalData,
00881 jac_trans,
00882 epetra_U,
00883 epetra_V,
00884 false));
00885
00886
00887
00888 Teuchos::RCP<Epetra_CrsMatrix> jac_trans_crs;
00889 if (includeUV && !use_P_For_Prec) {
00890 jac_trans_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(jac_trans);
00891 if (jac_trans_crs != Teuchos::null) {
00892 updateJacobianForPreconditioner(*U_trans, *V_trans, *jac_trans_crs);
00893 }
00894 }
00895
00896
00897 if (use_P_For_Prec)
00898 tls_strategy->setJacobianTransposeOperator(op);
00899
00900
00901 tls_strategy->createTransposePreconditioner(solution_vec, params);
00902
00903
00904 if (includeUV && !use_P_For_Prec && jac_trans_crs != Teuchos::null) {
00905
00906 grp->setX(solution_vec);
00907 if (isComplex)
00908 grp->computeComplex(omega);
00909 else
00910 grp->computeJacobian();
00911 tls_strategy->createJacobianTranspose();
00912 }
00913
00914
00915 if (!use_P_For_Prec)
00916 tls_strategy->setJacobianTransposeOperator(op);
00917
00918
00919 int m = X.numVectors();
00920 X.init(0.0);
00921 for (int i=0; i<m; i++) {
00922 bool stat =
00923 tls_strategy->applyJacobianTransposeInverse(
00924 params,
00925 dynamic_cast<const NOX::Epetra::Vector&>((*cRHS)[i]),
00926 dynamic_cast<NOX::Epetra::Vector&>(X[i]));
00927 if (stat == true)
00928 status = NOX::Abstract::Group::Ok;
00929 else
00930 status = NOX::Abstract::Group::NotConverged;
00931 finalStatus =
00932 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00933 finalStatus,
00934 callingFunction);
00935
00936
00937 }
00938
00939 qrFact.applyCompactWY(house_p_trans, *house_x_trans, T_trans, Y, X, isZeroG,
00940 false, false);
00941
00942 epetraOp->SetUseTranspose(false);
00943
00944
00945 linSys->setJacobianOperatorForSolve(epetraOp);
00946 linSys->destroyPreconditioner();
00947
00948 return finalStatus;
00949 }
00950
00951 NOX::Abstract::Group::ReturnType
00952 LOCA::BorderedSolver::EpetraHouseholder::computeUV(
00953 const NOX::Abstract::MultiVector::DenseMatrix& Y1,
00954 const NOX::Abstract::MultiVector& Y2,
00955 const NOX::Abstract::MultiVector::DenseMatrix& T,
00956 const NOX::Abstract::MultiVector& A,
00957 NOX::Abstract::MultiVector& UU,
00958 NOX::Abstract::MultiVector& VV,
00959 bool use_jac_transpose)
00960 {
00961
00962 const NOX::Epetra::MultiVector& Y2_epetra =
00963 dynamic_cast<const NOX::Epetra::MultiVector&>(Y2);
00964 NOX::Epetra::MultiVector& UU_epetra =
00965 dynamic_cast<NOX::Epetra::MultiVector&>(UU);
00966 bool use_trans = epetraOp->UseTranspose();
00967 epetraOp->SetUseTranspose(use_jac_transpose);
00968 int status = epetraOp->Apply(Y2_epetra.getEpetraMultiVector(),
00969 UU_epetra.getEpetraMultiVector());
00970 epetraOp->SetUseTranspose(use_trans);
00971
00972 UU.update(Teuchos::NO_TRANS, 1.0, A, Y1, 1.0);
00973 VV.update(Teuchos::TRANS, 1.0, Y2, T, 0.0);
00974
00975 if (status == 0)
00976 return NOX::Abstract::Group::Ok;
00977 else
00978 return NOX::Abstract::Group::Failed;
00979 }
00980
00981 void
00982 LOCA::BorderedSolver::EpetraHouseholder::updateJacobianForPreconditioner(
00983 const NOX::Abstract::MultiVector& UU,
00984 const NOX::Abstract::MultiVector& VV,
00985 Epetra_CrsMatrix& jac) const
00986 {
00987
00988 int numMyRows = jac.NumMyRows();
00989
00990 int numEntries;
00991 int U_LDA;
00992 int V_LDA;
00993 int *row_indices;
00994 double *row_values;
00995 double *U_values;
00996 double *V_values;
00997 double val;
00998
00999
01000 const NOX::Epetra::MultiVector nox_epetra_U =
01001 dynamic_cast<const NOX::Epetra::MultiVector&>(UU);
01002 const Epetra_MultiVector& epetra_U = nox_epetra_U.getEpetraMultiVector();
01003 const NOX::Epetra::MultiVector& nox_epetra_V =
01004 dynamic_cast<const NOX::Epetra::MultiVector&>(VV);
01005 const Epetra_MultiVector& epetra_V = nox_epetra_V.getEpetraMultiVector();
01006 epetra_U.ExtractView(&U_values, &U_LDA);
01007 epetra_V.ExtractView(&V_values, &V_LDA);
01008
01009 const Epetra_BlockMap& v_map = epetra_V.Map();
01010 const Epetra_BlockMap& u_map = epetra_V.Map();
01011 const Epetra_BlockMap& row_map = jac.RowMap();
01012 const Epetra_BlockMap& col_map = jac.ColMap();
01013
01014 for (int row=0; row<numMyRows; row++) {
01015
01016
01017 jac.ExtractMyRowView(row, numEntries, row_values, row_indices);
01018
01019 int row_gid = row_map.GID(row);
01020 int u_row_lid = u_map.LID(row_gid);
01021
01022 for (int col=0; col<numEntries; col++) {
01023
01024
01025 int col_gid = col_map.GID(row_indices[col]);
01026 if (v_map.MyGID(col_gid)) {
01027 int v_row_lid = v_map.LID(col_gid);
01028
01029
01030 val = 0;
01031 for (int k=0; k<numConstraints; k++)
01032 val += U_values[k*U_LDA+u_row_lid]*V_values[k*V_LDA+v_row_lid];
01033
01034
01035 row_values[col] += val;
01036
01037 }
01038
01039 }
01040 }
01041 }
01042
01043 Teuchos::RCP<NOX::Abstract::MultiVector>
01044 LOCA::BorderedSolver::EpetraHouseholder::createBlockMV(
01045 const NOX::Abstract::MultiVector& v) const
01046 {
01047 #ifdef HAVE_NOX_EPETRAEXT
01048 const LOCA::Hopf::ComplexMultiVector& cv =
01049 dynamic_cast<const LOCA::Hopf::ComplexMultiVector&>(v);
01050 Teuchos::RCP<const NOX::Abstract::MultiVector> v_real =
01051 cv.getRealMultiVec();
01052 Teuchos::RCP<const NOX::Abstract::MultiVector> v_imag =
01053 cv.getImagMultiVec();
01054 Teuchos::RCP<const NOX::Epetra::MultiVector> nox_epetra_v_real =
01055 Teuchos::rcp_dynamic_cast<const NOX::Epetra::MultiVector>(v_real);
01056 Teuchos::RCP<const NOX::Epetra::MultiVector> nox_epetra_v_imag =
01057 Teuchos::rcp_dynamic_cast<const NOX::Epetra::MultiVector>(v_imag);
01058
01059 Teuchos::RCP<EpetraExt::BlockMultiVector> epetra_v =
01060 Teuchos::rcp(new EpetraExt::BlockMultiVector(*baseMap, *globalMap,
01061 v.numVectors()));
01062 epetra_v->LoadBlockValues(nox_epetra_v_real->getEpetraMultiVector(), 0);
01063 epetra_v->LoadBlockValues(nox_epetra_v_imag->getEpetraMultiVector(), 1);
01064
01065 return Teuchos::rcp(new NOX::Epetra::MultiVector(
01066 epetra_v,
01067 NOX::DeepCopy,
01068 NOX::Epetra::MultiVector::CreateView));
01069 #else
01070 globalData->locaErrorCheck->throwError("LOCA::BorderedSolver::EpetraHouseholder::createBlockMV()",
01071 "Must have EpetraExt support for Hopf tracking. Configure trilinos with --enable-epetraext");
01072 return Teuchos::null;
01073 #endif
01074 }
01075
01076 void
01077 LOCA::BorderedSolver::EpetraHouseholder::setBlockMV(
01078 const NOX::Abstract::MultiVector& bv,
01079 NOX::Abstract::MultiVector& v) const
01080 {
01081 #ifdef HAVE_NOX_EPETRAEXT
01082 LOCA::Hopf::ComplexMultiVector& cv =
01083 dynamic_cast<LOCA::Hopf::ComplexMultiVector&>(v);
01084 Teuchos::RCP<NOX::Abstract::MultiVector> v_real =
01085 cv.getRealMultiVec();
01086 Teuchos::RCP<NOX::Abstract::MultiVector> v_imag =
01087 cv.getImagMultiVec();
01088 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_v_real =
01089 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(v_real);
01090 Teuchos::RCP<NOX::Epetra::MultiVector> nox_epetra_v_imag =
01091 Teuchos::rcp_dynamic_cast<NOX::Epetra::MultiVector>(v_imag);
01092
01093 const NOX::Epetra::MultiVector& nox_epetra_bv =
01094 dynamic_cast<const NOX::Epetra::MultiVector&>(bv);
01095 const Epetra_MultiVector& epetra_bv = nox_epetra_bv.getEpetraMultiVector();
01096 const EpetraExt::BlockMultiVector& block_bv =
01097 dynamic_cast<const EpetraExt::BlockMultiVector&>(epetra_bv);
01098
01099 block_bv.ExtractBlockValues(nox_epetra_v_real->getEpetraMultiVector(), 0);
01100 block_bv.ExtractBlockValues(nox_epetra_v_imag->getEpetraMultiVector(), 1);
01101 #else
01102 globalData->locaErrorCheck->throwError("LOCA::BorderedSolver::EpetraHouseholder::setBlockMV()",
01103 "Must have EpetraExt support for Hopf tracking. Configure trilinos with --enable-epetraext");
01104 #endif
01105 }