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_Bordering.H"
00043 #include "LOCA_BorderedSolver_AbstractOperator.H"
00044 #include "LOCA_GlobalData.H"
00045 #include "LOCA_ErrorCheck.H"
00046 #include "LOCA_MultiContinuation_ConstraintInterface.H"
00047 #include "LOCA_MultiContinuation_ConstraintInterfaceMVDX.H"
00048 #include "Teuchos_LAPACK.hpp"
00049 #include "LOCA_BorderedSolver_LowerTriangularBlockElimination.H"
00050 #include "LOCA_BorderedSolver_UpperTriangularBlockElimination.H"
00051 #include "LOCA_Abstract_TransposeSolveGroup.H"
00052
00053 LOCA::BorderedSolver::Bordering::Bordering(
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 op(),
00060 A(),
00061 B(),
00062 C(),
00063 isZeroA(true),
00064 isZeroB(true),
00065 isZeroC(true),
00066 isZeroF(true),
00067 isZeroG(true)
00068 {
00069 }
00070
00071 LOCA::BorderedSolver::Bordering::~Bordering()
00072 {
00073 }
00074
00075 void
00076 LOCA::BorderedSolver::Bordering::setMatrixBlocks(
00077 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& oper,
00078 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00079 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00080 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC)
00081 {
00082 op = oper;
00083 A = blockA;
00084 B = blockB;
00085 C = blockC;
00086
00087 isZeroA = (A.get() == NULL);
00088 isZeroB = B->isDXZero();
00089 isZeroC = (C.get() == NULL);
00090
00091
00092 if (isZeroB && isZeroC)
00093 globalData->locaErrorCheck->throwError(
00094 "LOCA::BorderedSolver::Bordering::setMatrixBlocks",
00095 "Blocks B and C cannot both be zero");
00096
00097
00098 if (isZeroA && isZeroC)
00099 globalData->locaErrorCheck->throwError(
00100 "LOCA::BorderedSolver::Bordering::setMatrixBlocks",
00101 "Blocks A and C cannot both be zero");
00102 }
00103
00104 NOX::Abstract::Group::ReturnType
00105 LOCA::BorderedSolver::Bordering::initForSolve()
00106 {
00107 return NOX::Abstract::Group::Ok;
00108 }
00109
00110 NOX::Abstract::Group::ReturnType
00111 LOCA::BorderedSolver::Bordering::initForTransposeSolve()
00112 {
00113 return NOX::Abstract::Group::Ok;
00114 }
00115
00116 NOX::Abstract::Group::ReturnType
00117 LOCA::BorderedSolver::Bordering::apply(
00118 const NOX::Abstract::MultiVector& X,
00119 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00120 NOX::Abstract::MultiVector& U,
00121 NOX::Abstract::MultiVector::DenseMatrix& V) const
00122 {
00123
00124 NOX::Abstract::Group::ReturnType status =
00125 op->apply(X, U);
00126
00127
00128 if (!isZeroA)
00129 U.update(Teuchos::NO_TRANS, 1.0, *A, Y, 1.0);
00130
00131
00132 if (!isZeroB)
00133 B->multiplyDX(1.0, X, V);
00134
00135
00136 if (!isZeroC) {
00137 int e;
00138 if (isZeroB)
00139 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00140 else
00141 e = V.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00142 if (e < 0)
00143 status = NOX::Abstract::Group::Failed;
00144 }
00145
00146 return status;
00147 }
00148
00149 NOX::Abstract::Group::ReturnType
00150 LOCA::BorderedSolver::Bordering::applyTranspose(
00151 const NOX::Abstract::MultiVector& X,
00152 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00153 NOX::Abstract::MultiVector& U,
00154 NOX::Abstract::MultiVector::DenseMatrix& V) const
00155 {
00156
00157 NOX::Abstract::Group::ReturnType status =
00158 op->applyTranspose(X, U);
00159
00160
00161 if (!isZeroA)
00162 B->addDX(Teuchos::NO_TRANS, 1.0, Y, 1.0, U);
00163
00164
00165 if (!isZeroB)
00166 X.multiply(1.0, *A, V);
00167
00168
00169 if (!isZeroC) {
00170 int e;
00171 if (isZeroB)
00172 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 0.0);
00173 else
00174 e = V.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *C, Y, 1.0);
00175 if (e < 0)
00176 status = NOX::Abstract::Group::Failed;
00177 }
00178
00179 return status;
00180 }
00181
00182 NOX::Abstract::Group::ReturnType
00183 LOCA::BorderedSolver::Bordering::applyInverse(
00184 Teuchos::ParameterList& params,
00185 const NOX::Abstract::MultiVector* F,
00186 const NOX::Abstract::MultiVector::DenseMatrix* G,
00187 NOX::Abstract::MultiVector& X,
00188 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00189 {
00190 string callingFunction =
00191 "LOCA::BorderedSolver::Bordering::applyInverse()";
00192 NOX::Abstract::Group::ReturnType status;
00193
00194 isZeroF = (F == NULL);
00195 isZeroG = (G == NULL);
00196
00197 if (isZeroA) {
00198 LOCA::BorderedSolver::LowerTriangularBlockElimination ltbe(globalData);
00199 status = ltbe.solve(params, *op, *B, *C, F, G, X, Y);
00200 }
00201
00202 else if (isZeroB) {
00203 LOCA::BorderedSolver::UpperTriangularBlockElimination utbe(globalData);
00204 status = utbe.solve(params, *op, A.get(), *C, F, G, X, Y);
00205
00206 }
00207
00208 else if (isZeroF)
00209 status = solveFZero(params, A.get(), B.get(), C.get(), G, X, Y);
00210
00211 else {
00212
00213 int numColsA = A->numVectors();
00214 int numColsF = F->numVectors();
00215
00216
00217 vector<int> indexF(numColsF);
00218 vector<int> indexA(numColsA);
00219 for (int i=0; i<numColsF; i++)
00220 indexF[i] = i;
00221 for (int i=0; i<numColsA; i++)
00222 indexA[i] = numColsF + i;
00223 int numColsRHS = numColsF + numColsA;
00224
00225
00226 Teuchos::RCP<NOX::Abstract::MultiVector> RHS =
00227 F->clone(numColsRHS);
00228 Teuchos::RCP<NOX::Abstract::MultiVector> LHS =
00229 X.clone(numColsRHS);
00230 Teuchos::RCP<NOX::Abstract::MultiVector> X1 =
00231 LHS->subView(indexF);
00232 RHS->setBlock(*F, indexF);
00233 RHS->setBlock(*A, indexA);
00234
00235
00236 status = solveContiguous(params, A.get(), B.get(), C.get(),
00237 indexF, indexA, RHS.get(), G, *LHS, Y);
00238 X = *X1;
00239
00240 }
00241
00242 return status;
00243 }
00244
00245 NOX::Abstract::Group::ReturnType
00246 LOCA::BorderedSolver::Bordering::applyInverseTranspose(
00247 Teuchos::ParameterList& params,
00248 const NOX::Abstract::MultiVector* F,
00249 const NOX::Abstract::MultiVector::DenseMatrix* G,
00250 NOX::Abstract::MultiVector& X,
00251 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00252 {
00253 string callingFunction =
00254 "LOCA::BorderedSolver::Bordering::applyInverseTranspose()";
00255 NOX::Abstract::Group::ReturnType status;
00256
00257 isZeroF = (F == NULL);
00258 isZeroG = (G == NULL);
00259
00260
00261 Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterfaceMVDX> B_mvdx;
00262 const NOX::Abstract::MultiVector* BB = NULL;
00263
00264 if (!isZeroB) {
00265 B_mvdx = Teuchos::rcp_dynamic_cast<const LOCA::MultiContinuation::ConstraintInterfaceMVDX>(B);
00266 if (B == Teuchos::null)
00267 globalData->locaErrorCheck->throwError(
00268 callingFunction,
00269 "Constraints object must be of type ConstraintInterfaceMVDX");
00270 BB = B_mvdx->getDX();
00271 }
00272
00273 if (isZeroA) {
00274 LOCA::BorderedSolver::UpperTriangularBlockElimination utbe(globalData);
00275 status = utbe.solveTranspose(params, *op, BB, *C, F, G, X, Y);
00276 }
00277
00278 else if (isZeroB) {
00279 LOCA::BorderedSolver::LowerTriangularBlockElimination ltbe(globalData);
00280 status = ltbe.solveTranspose(params, *op, *A, *C, F, G, X, Y);
00281
00282 }
00283
00284 else if (isZeroF)
00285 status = solveFZeroTrans(params, A.get(), BB, C.get(), G, X, Y);
00286
00287 else {
00288
00289 int numColsB = BB->numVectors();
00290 int numColsF = F->numVectors();
00291
00292
00293 vector<int> indexF(numColsF);
00294 vector<int> indexB(numColsB);
00295 for (int i=0; i<numColsF; i++)
00296 indexF[i] = i;
00297 for (int i=0; i<numColsB; i++)
00298 indexB[i] = numColsF + i;
00299 int numColsRHS = numColsF + numColsB;
00300
00301
00302 Teuchos::RCP<NOX::Abstract::MultiVector> RHS =
00303 F->clone(numColsRHS);
00304 Teuchos::RCP<NOX::Abstract::MultiVector> LHS =
00305 X.clone(numColsRHS);
00306 Teuchos::RCP<NOX::Abstract::MultiVector> X1 =
00307 LHS->subView(indexF);
00308 RHS->setBlock(*F, indexF);
00309 RHS->setBlock(*BB, indexB);
00310
00311
00312 status = solveContiguousTrans(params, A.get(), BB, C.get(),
00313 indexF, indexB, RHS.get(), G, *LHS, Y);
00314 X = *X1;
00315
00316 }
00317
00318 return status;
00319 }
00320
00321
00322
00323
00324
00325
00326 NOX::Abstract::Group::ReturnType
00327 LOCA::BorderedSolver::Bordering::solveFZero(
00328 Teuchos::ParameterList& params,
00329 const NOX::Abstract::MultiVector* AA,
00330 const LOCA::MultiContinuation::ConstraintInterface* BB,
00331 const NOX::Abstract::MultiVector::DenseMatrix* CC,
00332 const NOX::Abstract::MultiVector::DenseMatrix* G,
00333 NOX::Abstract::MultiVector& X,
00334 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00335 {
00336 string callingFunction =
00337 "LOCA::BorderedSolver::Bordering::solveFZero()";
00338 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00339 NOX::Abstract::Group::ReturnType status;
00340
00341
00342 if (isZeroG) {
00343 X.init(0.0);
00344 Y.putScalar(0.0);
00345 return finalStatus;
00346 }
00347
00348 Teuchos::RCP<NOX::Abstract::MultiVector> Xt =
00349 AA->clone(NOX::ShapeCopy);
00350
00351
00352 status = op->applyInverse(params, *AA, *Xt);
00353 finalStatus =
00354 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00355 callingFunction);
00356
00357
00358 NOX::Abstract::MultiVector::DenseMatrix t(BB->numConstraints(),
00359 Xt->numVectors());
00360 BB->multiplyDX(-1.0, *Xt, t);
00361
00362
00363 if (!isZeroC)
00364 t += *CC;
00365
00366
00367 Y.assign(*G);
00368 Teuchos::LAPACK<int,double> L;
00369 int *ipiv = new int[t.numRows()];
00370 int info;
00371 L.GESV(t.numRows(), Y.numCols(), t.values(), t.stride(), ipiv,
00372 Y.values(), Y.stride(), &info);
00373 delete [] ipiv;
00374 if (info != 0) {
00375 status = NOX::Abstract::Group::Failed;
00376 finalStatus =
00377 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00378 finalStatus,
00379 callingFunction);
00380 }
00381
00382
00383 X.update(Teuchos::NO_TRANS, -1.0, *Xt, Y, -0.0);
00384
00385 return finalStatus;
00386 }
00387
00388
00389
00390 NOX::Abstract::Group::ReturnType
00391 LOCA::BorderedSolver::Bordering::solveContiguous(
00392 Teuchos::ParameterList& params,
00393 const NOX::Abstract::MultiVector* AA,
00394 const LOCA::MultiContinuation::ConstraintInterface* BB,
00395 const NOX::Abstract::MultiVector::DenseMatrix* CC,
00396 vector<int>& indexF,
00397 vector<int>& indexA,
00398 const NOX::Abstract::MultiVector* F,
00399 const NOX::Abstract::MultiVector::DenseMatrix* G,
00400 NOX::Abstract::MultiVector& X,
00401 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00402 {
00403 string callingFunction =
00404 "LOCA::BorderedSolver::Bordering::solveContiguous()";
00405 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00406 NOX::Abstract::Group::ReturnType status;
00407
00408
00409 status = op->applyInverse(params, *F, X);
00410 finalStatus =
00411 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00412 callingFunction);
00413 Teuchos::RCP<NOX::Abstract::MultiVector> X1 = X.subView(indexF);
00414 Teuchos::RCP<NOX::Abstract::MultiVector> X2 = X.subView(indexA);
00415
00416
00417 BB->multiplyDX(-1.0, *X1, Y);
00418
00419
00420 NOX::Abstract::MultiVector::DenseMatrix t2(BB->numConstraints(),
00421 X2->numVectors());
00422 BB->multiplyDX(-1.0, *X2, t2);
00423
00424
00425 if (!isZeroG)
00426 Y += *G;
00427
00428
00429 if (!isZeroC)
00430 t2 += *CC;
00431
00432
00433 Teuchos::LAPACK<int,double> L;
00434 int *ipiv = new int[t2.numRows()];
00435 int info;
00436 L.GESV(t2.numRows(), Y.numCols(), t2.values(), t2.stride(), ipiv,
00437 Y.values(), Y.stride(), &info);
00438 delete [] ipiv;
00439 if (info != 0) {
00440 status = NOX::Abstract::Group::Failed;
00441 finalStatus =
00442 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00443 finalStatus,
00444 callingFunction);
00445 }
00446
00447
00448 X1->update(Teuchos::NO_TRANS, -1.0, *X2, Y, 1.0);
00449
00450 return finalStatus;
00451 }
00452
00453
00454
00455
00456
00457
00458 NOX::Abstract::Group::ReturnType
00459 LOCA::BorderedSolver::Bordering::solveFZeroTrans(
00460 Teuchos::ParameterList& params,
00461 const NOX::Abstract::MultiVector* AA,
00462 const NOX::Abstract::MultiVector* BB,
00463 const NOX::Abstract::MultiVector::DenseMatrix* CC,
00464 const NOX::Abstract::MultiVector::DenseMatrix* G,
00465 NOX::Abstract::MultiVector& X,
00466 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00467 {
00468 string callingFunction =
00469 "LOCA::BorderedSolver::Bordering::solveFTransZero()";
00470 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00471 NOX::Abstract::Group::ReturnType status;
00472
00473
00474 if (isZeroG) {
00475 X.init(0.0);
00476 Y.putScalar(0.0);
00477 return finalStatus;
00478 }
00479
00480 Teuchos::RCP<NOX::Abstract::MultiVector> Xt =
00481 BB->clone(NOX::ShapeCopy);
00482
00483
00484 status = op->applyInverseTranspose(params, *BB, *Xt);
00485 finalStatus =
00486 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00487 callingFunction);
00488
00489
00490 NOX::Abstract::MultiVector::DenseMatrix t(AA->numVectors(),
00491 Xt->numVectors());
00492 Xt->multiply(-1.0, *AA, t);
00493
00494
00495 if (!isZeroC)
00496 for (int i=0; i<t.numRows(); i++)
00497 for (int j=0; j<t.numCols(); j++)
00498 t(i,j) += (*CC)(j,i);
00499
00500
00501 Y.assign(*G);
00502 Teuchos::LAPACK<int,double> L;
00503 int *ipiv = new int[t.numRows()];
00504 int info;
00505 L.GESV(t.numRows(), Y.numCols(), t.values(), t.stride(), ipiv,
00506 Y.values(), Y.stride(), &info);
00507 delete [] ipiv;
00508 if (info != 0) {
00509 status = NOX::Abstract::Group::Failed;
00510 finalStatus =
00511 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00512 finalStatus,
00513 callingFunction);
00514 }
00515
00516
00517 X.update(Teuchos::NO_TRANS, -1.0, *Xt, Y, -0.0);
00518
00519 return finalStatus;
00520 }
00521
00522
00523
00524 NOX::Abstract::Group::ReturnType
00525 LOCA::BorderedSolver::Bordering::solveContiguousTrans(
00526 Teuchos::ParameterList& params,
00527 const NOX::Abstract::MultiVector* AA,
00528 const NOX::Abstract::MultiVector* BB,
00529 const NOX::Abstract::MultiVector::DenseMatrix* CC,
00530 vector<int>& indexF,
00531 vector<int>& indexB,
00532 const NOX::Abstract::MultiVector* F,
00533 const NOX::Abstract::MultiVector::DenseMatrix* G,
00534 NOX::Abstract::MultiVector& X,
00535 NOX::Abstract::MultiVector::DenseMatrix& Y) const
00536 {
00537 string callingFunction =
00538 "LOCA::BorderedSolver::Bordering::solveContiguousTrans()";
00539 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00540 NOX::Abstract::Group::ReturnType status;
00541
00542
00543 status = op->applyInverseTranspose(params, *F, X);
00544 finalStatus =
00545 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00546 callingFunction);
00547 Teuchos::RCP<NOX::Abstract::MultiVector> X1 = X.subView(indexF);
00548 Teuchos::RCP<NOX::Abstract::MultiVector> X2 = X.subView(indexB);
00549
00550
00551 X1->multiply(-1.0, *AA, Y);
00552
00553
00554 NOX::Abstract::MultiVector::DenseMatrix t2(AA->numVectors(),
00555 X2->numVectors());
00556 X2->multiply(-1.0, *AA, t2);
00557
00558
00559 if (!isZeroG)
00560 Y += *G;
00561
00562
00563 if (!isZeroC)
00564 for (int i=0; i<t2.numRows(); i++)
00565 for (int j=0; j<t2.numCols(); j++)
00566 t2(i,j) += (*CC)(j,i);
00567
00568
00569 Teuchos::LAPACK<int,double> L;
00570 int *ipiv = new int[t2.numRows()];
00571 int info;
00572 L.GESV(t2.numRows(), Y.numCols(), t2.values(), t2.stride(), ipiv,
00573 Y.values(), Y.stride(), &info);
00574 delete [] ipiv;
00575 if (info != 0) {
00576 status = NOX::Abstract::Group::Failed;
00577 finalStatus =
00578 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00579 finalStatus,
00580 callingFunction);
00581 }
00582
00583
00584 X1->update(Teuchos::NO_TRANS, -1.0, *X2, Y, 1.0);
00585
00586 return finalStatus;
00587 }