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_LAPACKDirectSolve.H"
00043 #include "LOCA_GlobalData.H"
00044 #include "LOCA_ErrorCheck.H"
00045 #include "LOCA_MultiContinuation_ConstraintInterfaceMVDX.H"
00046 #include "LOCA_LAPACK.H"
00047 #include "Teuchos_LAPACK.hpp"
00048 #include "LOCA_BorderedSolver_JacobianOperator.H"
00049 #include "LOCA_BorderedSolver_ComplexOperator.H"
00050 #include "LOCA_Hopf_ComplexMultiVector.H"
00051
00052 LOCA::BorderedSolver::LAPACKDirectSolve::LAPACKDirectSolve(
00053 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00054 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00055 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams):
00056 globalData(global_data),
00057 solverParams(slvrParams),
00058 grp(),
00059 op(),
00060 A(),
00061 B(),
00062 C(),
00063 augJacSolver(),
00064 augComplexSolver(),
00065 n(0),
00066 m(0),
00067 N(0),
00068 isZeroA(true),
00069 isZeroB(true),
00070 isZeroC(true),
00071 isZeroF(true),
00072 isZeroG(true),
00073 isComplex(false)
00074 {
00075 }
00076
00077 LOCA::BorderedSolver::LAPACKDirectSolve::~LAPACKDirectSolve()
00078 {
00079 }
00080
00081 void
00082 LOCA::BorderedSolver::LAPACKDirectSolve::setMatrixBlocks(
00083 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& oper,
00084 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00085 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00086 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC)
00087 {
00088 string callingFunction =
00089 "LOCA::BorderedSolver::LAPACKDirectSolve::setMatrixBlocks()";
00090
00091
00092 op = oper;
00093 A = blockA;
00094
00095 B = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(blockB);
00096 if (B.get() == NULL)
00097 globalData->locaErrorCheck->throwError(
00098 callingFunction,
00099 string("Constraint argument is not of type\n") +
00100 string("LOCA::MultiContinuation::ConstraintInterfaceMVDX!\n") +
00101 string("The LAPACK Direct Solve bordered solver method can\n") +
00102 string("only be used with constraints that support obtaining\n") +
00103 string("the constraint derivative as a multivector."));
00104
00105 C = blockC;
00106
00107 isZeroA = (A.get() == NULL);
00108 isZeroB = B->isDXZero();
00109 isZeroC = (C.get() == NULL);
00110
00111
00112 if (isZeroB && isZeroC)
00113 globalData->locaErrorCheck->throwError(
00114 callingFunction,
00115 "Blocks B and C cannot both be zero");
00116
00117
00118 if (isZeroA && isZeroC)
00119 globalData->locaErrorCheck->throwError(
00120 callingFunction,
00121 "Blocks A and C cannot both be zero");
00122
00123
00124 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> jacOp =
00125 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::JacobianOperator>(op);
00126 Teuchos::RCP<const LOCA::BorderedSolver::ComplexOperator> complexOp =
00127 Teuchos::rcp_dynamic_cast<const LOCA::BorderedSolver::ComplexOperator>(op);
00128
00129 if (jacOp != Teuchos::null) {
00130 const NOX::LAPACK::Vector *v;
00131 const NOX::Abstract::MultiVector *BV;
00132
00133 isComplex = false;
00134
00135 grp = Teuchos::rcp_dynamic_cast<const LOCA::LAPACK::Group>(jacOp->getGroup());
00136 if (grp.get() == NULL)
00137 globalData->locaErrorCheck->throwError(
00138 callingFunction,
00139 string("Group argument is not of type LOCA::LAPACK::Group!\n") +
00140 string("The LAPACK Direct Solve bordered solver method can\n") +
00141 string("only be used with LAPACK groups."));
00142
00143
00144 const NOX::LAPACK::Matrix<double>& J = grp->getJacobianMatrix();
00145 n = J.numRows();
00146
00147
00148 if (!isZeroA)
00149 m = A->numVectors();
00150 else
00151 m = C->numCols();
00152
00153
00154 if (n+m != N) {
00155 N = n+m;
00156 augJacSolver = Teuchos::rcp(new NOX::LAPACK::LinearSolver<double>(N));
00157 }
00158 else {
00159 augJacSolver->reset();
00160 }
00161 NOX::LAPACK::Matrix<double>& augmentedJ = augJacSolver->getMatrix();
00162
00163
00164 for (int j=0; j<n; j++)
00165 for (int i=0; i<n; i++)
00166 augmentedJ(i,j) = J(i,j);
00167
00168
00169 if (isZeroA) {
00170 for (int j=0; j<m; j++)
00171 for (int i=0; i<n; i++)
00172 augmentedJ(i,j+n) = 0.0;
00173 }
00174 else {
00175 for (int j=0; j<m; j++) {
00176 v = dynamic_cast<const NOX::LAPACK::Vector*>(&(*A)[j]);
00177 for (int i=0; i<n; i++)
00178 augmentedJ(i,j+n) = (*v)(i);
00179 }
00180 }
00181
00182
00183 if (isZeroB) {
00184 for (int i=0; i<m; i++)
00185 for (int j=0; j<n; j++)
00186 augmentedJ(i+n,j) = 0.0;
00187 }
00188 else {
00189 BV = B->getDX();
00190 for (int i=0; i<m; i++) {
00191 v = dynamic_cast<const NOX::LAPACK::Vector*>(&(*BV)[i]);
00192 for (int j=0; j<n; j++)
00193 augmentedJ(i+n,j) = (*v)(j);
00194 }
00195 }
00196
00197
00198 if (isZeroC) {
00199 for (int j=0; j<m; j++)
00200 for (int i=0; i<m; i++)
00201 augmentedJ(i+n,j+n) = 0.0;
00202 }
00203 else {
00204 for (int j=0; j<m; j++)
00205 for (int i=0; i<m; i++)
00206 augmentedJ(i+n,j+n) = (*C)(i,j);
00207 }
00208
00209 }
00210 else if (complexOp != Teuchos::null) {
00211 Teuchos::RCP<const LOCA::Hopf::ComplexMultiVector> cA;
00212 Teuchos::RCP<const LOCA::Hopf::ComplexMultiVector> cB;
00213 Teuchos::RCP<const NOX::Abstract::MultiVector> A_real;
00214 Teuchos::RCP<const NOX::Abstract::MultiVector> A_imag;
00215 Teuchos::RCP<const NOX::Abstract::MultiVector> B_real;
00216 Teuchos::RCP<const NOX::Abstract::MultiVector> B_imag;
00217 const NOX::LAPACK::Vector *v1, *v2;
00218 Teuchos::RCP<const NOX::Abstract::MultiVector> BV;
00219
00220 isComplex = true;
00221
00222 if (!isZeroA) {
00223 cA =
00224 Teuchos::rcp_dynamic_cast<const LOCA::Hopf::ComplexMultiVector>(A,
00225 true);
00226 A_real = cA->getRealMultiVec();
00227 A_imag = cA->getImagMultiVec();
00228 }
00229 if (!isZeroB) {
00230 BV = Teuchos::rcp(B->getDX(),false);
00231 cB =
00232 Teuchos::rcp_dynamic_cast<const LOCA::Hopf::ComplexMultiVector>(BV,
00233 true);
00234 B_real = cB->getRealMultiVec();
00235 B_imag = cB->getImagMultiVec();
00236 }
00237
00238 grp = Teuchos::rcp_dynamic_cast<const LOCA::LAPACK::Group>(complexOp->getGroup());
00239 if (grp.get() == NULL)
00240 globalData->locaErrorCheck->throwError(
00241 callingFunction,
00242 string("Group argument is not of type LOCA::LAPACK::Group!\n") +
00243 string("The LAPACK Direct Solve bordered solver method can\n") +
00244 string("only be used with LAPACK groups."));
00245
00246
00247 if (!isZeroA)
00248 m = A->numVectors()/2;
00249 else
00250 m = C->numCols()/2;
00251
00252
00253 const NOX::LAPACK::Matrix< std::complex<double> >& mat =
00254 grp->getComplexMatrix();
00255 n = mat.numRows();
00256
00257
00258 if (n+m != N) {
00259 N = n+m;
00260 augComplexSolver =
00261 Teuchos::rcp(new NOX::LAPACK::LinearSolver< std::complex<double> >(N));
00262 }
00263 else {
00264 augComplexSolver->reset();
00265 }
00266 NOX::LAPACK::Matrix< std::complex<double> >& augmentedMat =
00267 augComplexSolver->getMatrix();
00268
00269
00270 for (int j=0; j<n; j++)
00271 for (int i=0; i<n; i++)
00272 augmentedMat(i,j) = mat(i,j);
00273
00274
00275 if (isZeroA) {
00276 for (int j=0; j<m; j++)
00277 for (int i=0; i<n; i++)
00278 augmentedMat(i,j+n) = 0.0;
00279 }
00280 else {
00281 for (int j=0; j<m; j++) {
00282 v1 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*A_real)[2*j]);
00283 v2 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*A_imag)[2*j]);
00284 for (int i=0; i<n; i++)
00285 augmentedMat(i,j+n) = std::complex<double>((*v1)(i), (*v2)(i));
00286 }
00287 }
00288
00289
00290 if (isZeroB) {
00291 for (int i=0; i<m; i++)
00292 for (int j=0; j<n; j++)
00293 augmentedMat(i+n,j) = 0.0;
00294 }
00295 else {
00296 for (int i=0; i<m; i++) {
00297 v1 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*B_real)[2*i]);
00298 v2 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*B_imag)[2*i]);
00299 for (int j=0; j<n; j++)
00300 augmentedMat(i+n,j) = std::complex<double>((*v1)(j), -(*v2)(j));
00301 }
00302 }
00303
00304
00305 if (isZeroC) {
00306 for (int j=0; j<m; j++)
00307 for (int i=0; i<m; i++)
00308 augmentedMat(i+n,j+n) = 0.0;
00309 }
00310 else {
00311 for (int j=0; j<m; j++)
00312 for (int i=0; i<m; i++)
00313 augmentedMat(i+n,j+n) = std::complex<double>((*C)(i,2*j),
00314 (*C)(i+m,2*j));
00315 }
00316
00317 }
00318
00319 else {
00320 globalData->locaErrorCheck->throwError(
00321 callingFunction,
00322 string("Op argument must be of type !\n") +
00323 string("LOCA::BorderedSolver::JacobianOperator or \n") +
00324 string("LOCA::BorderedSolver::ComplexOperator."));
00325 }
00326 }
00327
00328 NOX::Abstract::Group::ReturnType
00329 LOCA::BorderedSolver::LAPACKDirectSolve::initForSolve()
00330 {
00331 return NOX::Abstract::Group::Ok;
00332 }
00333
00334 NOX::Abstract::Group::ReturnType
00335 LOCA::BorderedSolver::LAPACKDirectSolve::initForTransposeSolve()
00336 {
00337 return NOX::Abstract::Group::Ok;
00338 }
00339
00340 NOX::Abstract::Group::ReturnType
00341 LOCA::BorderedSolver::LAPACKDirectSolve::apply(
00342 const NOX::Abstract::MultiVector& X,
00343 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00344 NOX::Abstract::MultiVector& U,
00345 NOX::Abstract::MultiVector::DenseMatrix& V) const
00346 {
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427 NOX::Abstract::Group::ReturnType status = op->apply(X, U);
00428
00429
00430 if (!isZeroA)
00431 U.update(Teuchos::NO_TRANS, 1.0, *A, Y, 1.0);
00432
00433
00434 if (!isZeroB)
00435 B->multiplyDX(1.0, X, V);
00436
00437
00438 if (!isZeroC) {
00439 int e;
00440 if (isZeroB)
00441 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00442 else
00443 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00444 if (e < 0)
00445 status = NOX::Abstract::Group::Failed;
00446 }
00447
00448 return status;
00449 }
00450
00451 NOX::Abstract::Group::ReturnType
00452 LOCA::BorderedSolver::LAPACKDirectSolve::applyTranspose(
00453 const NOX::Abstract::MultiVector& X,
00454 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00455 NOX::Abstract::MultiVector& U,
00456 NOX::Abstract::MultiVector::DenseMatrix& V) const
00457 {
00458
00459 NOX::Abstract::Group::ReturnType status = op->applyTranspose(X, U);
00460
00461
00462 if (!isZeroA)
00463 B->addDX(Teuchos::NO_TRANS, 1.0, Y, 1.0, U);
00464
00465
00466 if (!isZeroB)
00467 X.multiply(1.0, *A, V);
00468
00469
00470 if (!isZeroC) {
00471 int e;
00472 if (isZeroB)
00473 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00474 else
00475 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00476 if (e < 0)
00477 status = NOX::Abstract::Group::Failed;
00478 }
00479
00480 return status;
00481 }
00482
00483 NOX::Abstract::Group::ReturnType
00484 LOCA::BorderedSolver::LAPACKDirectSolve::applyInverse(
00485 Teuchos::ParameterList& params,
00486 const NOX::Abstract::MultiVector* F,
00487 const NOX::Abstract::MultiVector::DenseMatrix* G,
00488 NOX::Abstract::MultiVector& X,
00489 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00490 {
00491 return solve(false, params, F, G, X, Y);
00492 }
00493
00494 NOX::Abstract::Group::ReturnType
00495 LOCA::BorderedSolver::LAPACKDirectSolve::applyInverseTranspose(
00496 Teuchos::ParameterList& params,
00497 const NOX::Abstract::MultiVector* F,
00498 const NOX::Abstract::MultiVector::DenseMatrix* G,
00499 NOX::Abstract::MultiVector& X,
00500 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00501 {
00502 return solve(true, params, F, G, X, Y);
00503 }
00504
00505 NOX::Abstract::Group::ReturnType
00506 LOCA::BorderedSolver::LAPACKDirectSolve::solve(
00507 bool trans,
00508 Teuchos::ParameterList& params,
00509 const NOX::Abstract::MultiVector* F,
00510 const NOX::Abstract::MultiVector::DenseMatrix* G,
00511 NOX::Abstract::MultiVector& X,
00512 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00513 {
00514 bool isZeroF = (F == NULL);
00515 bool isZeroG = (G == NULL);
00516
00517
00518 if (isZeroF && isZeroG) {
00519 X.init(0.0);
00520 Y.putScalar(0.0);
00521 return NOX::Abstract::Group::Ok;
00522 }
00523
00524 int numColsRHS;
00525 if (!isZeroF)
00526 numColsRHS = F->numVectors();
00527 else
00528 numColsRHS = G->numCols();
00529
00530 bool res;
00531 if (!isComplex) {
00532 const NOX::LAPACK::Vector *v;
00533 NOX::LAPACK::Vector *w;
00534
00535
00536 NOX::LAPACK::Matrix<double> RHS(N,numColsRHS);
00537 if (isZeroF) {
00538 for (int j=0; j<numColsRHS; j++)
00539 for (int i=0; i<n; i++)
00540 RHS(i,j) = 0.0;
00541 }
00542 else {
00543 for (int j=0; j<numColsRHS; j++) {
00544 v = dynamic_cast<const NOX::LAPACK::Vector*>(&(*F)[j]);
00545 for (int i=0; i<n; i++)
00546 RHS(i,j) = (*v)(i);
00547 }
00548 }
00549 if (isZeroG) {
00550 for (int j=0; j<numColsRHS; j++)
00551 for (int i=0; i<m; i++)
00552 RHS(i+n,j) = 0.0;
00553 }
00554 else {
00555 for (int j=0; j<numColsRHS; j++)
00556 for (int i=0; i<m; i++)
00557 RHS(i+n,j) = (*G)(i,j);
00558 }
00559
00560
00561 res = augJacSolver->solve(trans, numColsRHS, &RHS(0,0));
00562
00563
00564 for (int j=0; j<numColsRHS; j++) {
00565 w = dynamic_cast<NOX::LAPACK::Vector*>(&X[j]);
00566 for (int i=0; i<n; i++)
00567 (*w)(i) = RHS(i,j);
00568 for (int i=0; i<m; i++)
00569 Y(i,j) = RHS(n+i,j);
00570 }
00571 }
00572
00573 else {
00574 const LOCA::Hopf::ComplexMultiVector* cF;
00575 LOCA::Hopf::ComplexMultiVector* cX;
00576 Teuchos::RCP<const NOX::Abstract::MultiVector> F_real;
00577 Teuchos::RCP<const NOX::Abstract::MultiVector> F_imag;
00578 Teuchos::RCP<NOX::Abstract::MultiVector> X_real;
00579 Teuchos::RCP<NOX::Abstract::MultiVector> X_imag;
00580 const NOX::LAPACK::Vector *v1, *v2;
00581 NOX::LAPACK::Vector *w1, *w2;
00582
00583 if (!isZeroF) {
00584 cF = dynamic_cast<const LOCA::Hopf::ComplexMultiVector*>(F);
00585 F_real = cF->getRealMultiVec();
00586 F_imag = cF->getImagMultiVec();
00587 }
00588 cX = dynamic_cast<LOCA::Hopf::ComplexMultiVector*>(&X);
00589 X_real = cX->getRealMultiVec();
00590 X_imag = cX->getImagMultiVec();
00591
00592
00593 NOX::LAPACK::Matrix< std::complex<double> > RHS(N,numColsRHS);
00594 if (isZeroF) {
00595 for (int j=0; j<numColsRHS; j++)
00596 for (int i=0; i<n; i++)
00597 RHS(i,j) = 0.0;
00598 }
00599 else {
00600 for (int j=0; j<numColsRHS; j++) {
00601 v1 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*F_real)[j]);
00602 v2 = dynamic_cast<const NOX::LAPACK::Vector*>(&(*F_imag)[j]);
00603 for (int i=0; i<n; i++)
00604 RHS(i,j) = std::complex<double>((*v1)(i), (*v2)(i));
00605 }
00606 }
00607 if (isZeroG) {
00608 for (int j=0; j<numColsRHS; j++)
00609 for (int i=0; i<m; i++)
00610 RHS(i+n,j) = 0.0;
00611 }
00612 else {
00613 for (int j=0; j<numColsRHS; j++)
00614 for (int i=0; i<m; i++)
00615 RHS(i+n,j) = std::complex<double>((*G)(i,j), (*G)(i+m,j));
00616 }
00617
00618
00619 res = augComplexSolver->solve(trans, numColsRHS, &RHS(0,0));
00620
00621
00622 for (int j=0; j<numColsRHS; j++) {
00623 w1 = dynamic_cast<NOX::LAPACK::Vector*>(&(*X_real)[j]);
00624 w2 = dynamic_cast<NOX::LAPACK::Vector*>(&(*X_imag)[j]);
00625 for (int i=0; i<n; i++) {
00626 (*w1)(i) = RHS(i,j).real();
00627 (*w2)(i) = RHS(i,j).imag();
00628 }
00629 for (int i=0; i<m; i++) {
00630 Y(i,j) = RHS(n+i,j).real();
00631 Y(i+m,j) = RHS(n+i,j).imag();
00632 }
00633 }
00634 }
00635
00636 if (!res)
00637 return NOX::Abstract::Group::Failed;
00638 else
00639 return NOX::Abstract::Group::Ok;
00640 }