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 #ifndef LOCA_BORDEREDSOLVER_EPETRAHOUSEHOLDER_H
00043 #define LOCA_BORDEREDSOLVER_EPETRAHOUSEHOLDER_H
00044
00045 #include "LOCA_BorderedSolver_AbstractStrategy.H"
00046
00047 #include "LOCA_BorderedSolver_HouseholderQR.H"
00048 #include "Teuchos_BLAS.hpp"
00049
00050
00051 class Epetra_Operator;
00052 class Epetra_CrsMatrix;
00053 class Epetra_BlockMap;
00054 namespace LOCA {
00055 class GlobalData;
00056 namespace Parameter {
00057 class SublistParser;
00058 }
00059 namespace MultiContinuation {
00060 class ConstraintInterfaceMVDX;
00061 }
00062 namespace Epetra {
00063 class Group;
00064 }
00065 }
00066 namespace NOX {
00067 namespace Epetra {
00068 class LinearSystem;
00069 }
00070 }
00071 namespace EpetraExt {
00072 class BlockMultiVector;
00073 }
00074
00075 namespace LOCA {
00076
00077 namespace BorderedSolver {
00078
00080
00257 class EpetraHouseholder : public LOCA::BorderedSolver::AbstractStrategy {
00258
00259 public:
00260
00262
00267 EpetraHouseholder(
00268 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00269 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00270 const Teuchos::RCP<Teuchos::ParameterList>& solverParams);
00271
00273 virtual ~EpetraHouseholder();
00274
00276
00281 virtual void setMatrixBlocks(
00282 const Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator>& op,
00283 const Teuchos::RCP<const NOX::Abstract::MultiVector>& blockA,
00284 const Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterface>& blockB,
00285 const Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix>& blockC);
00286
00288
00292 virtual NOX::Abstract::Group::ReturnType
00293 initForSolve();
00294
00296
00300 virtual NOX::Abstract::Group::ReturnType
00301 initForTransposeSolve();
00302
00327 virtual NOX::Abstract::Group::ReturnType
00328 apply(const NOX::Abstract::MultiVector& X,
00329 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00330 NOX::Abstract::MultiVector& U,
00331 NOX::Abstract::MultiVector::DenseMatrix& V) const;
00332
00357 virtual NOX::Abstract::Group::ReturnType
00358 applyTranspose(const NOX::Abstract::MultiVector& X,
00359 const NOX::Abstract::MultiVector::DenseMatrix& Y,
00360 NOX::Abstract::MultiVector& U,
00361 NOX::Abstract::MultiVector::DenseMatrix& V) const;
00362
00376 virtual NOX::Abstract::Group::ReturnType
00377 applyInverse(Teuchos::ParameterList& params,
00378 const NOX::Abstract::MultiVector* F,
00379 const NOX::Abstract::MultiVector::DenseMatrix* G,
00380 NOX::Abstract::MultiVector& X,
00381 NOX::Abstract::MultiVector::DenseMatrix& Y) const;
00382
00389 virtual NOX::Abstract::Group::ReturnType
00390 applyInverseTranspose(Teuchos::ParameterList& params,
00391 const NOX::Abstract::MultiVector* F,
00392 const NOX::Abstract::MultiVector::DenseMatrix* G,
00393 NOX::Abstract::MultiVector& X,
00394 NOX::Abstract::MultiVector::DenseMatrix& Y) const;
00395
00396 protected:
00397
00402 virtual NOX::Abstract::Group::ReturnType
00403 solve(Teuchos::ParameterList& params,
00404 const NOX::Abstract::MultiVector* F,
00405 const NOX::Abstract::MultiVector::DenseMatrix* G,
00406 NOX::Abstract::MultiVector& X,
00407 NOX::Abstract::MultiVector::DenseMatrix& Y) const;
00408
00412 virtual NOX::Abstract::Group::ReturnType
00413 solveTranspose(Teuchos::ParameterList& params,
00414 const NOX::Abstract::MultiVector* F,
00415 const NOX::Abstract::MultiVector::DenseMatrix* G,
00416 NOX::Abstract::MultiVector& X,
00417 NOX::Abstract::MultiVector::DenseMatrix& Y) const;
00418
00420 NOX::Abstract::Group::ReturnType
00421 computeUV(const NOX::Abstract::MultiVector::DenseMatrix& Y1,
00422 const NOX::Abstract::MultiVector& Y2,
00423 const NOX::Abstract::MultiVector::DenseMatrix& T,
00424 const NOX::Abstract::MultiVector& A,
00425 NOX::Abstract::MultiVector& U,
00426 NOX::Abstract::MultiVector& V,
00427 bool use_jac_transpose);
00428
00433 void updateJacobianForPreconditioner(const NOX::Abstract::MultiVector& U,
00434 const NOX::Abstract::MultiVector& V,
00435 Epetra_CrsMatrix& jac) const;
00436
00437 Teuchos::RCP<NOX::Abstract::MultiVector>
00438 createBlockMV(const NOX::Abstract::MultiVector& v) const;
00439
00440 void
00441 setBlockMV(const NOX::Abstract::MultiVector& bv,
00442 NOX::Abstract::MultiVector& v) const;
00443
00444 private:
00445
00447 EpetraHouseholder(const EpetraHouseholder&);
00448
00450 EpetraHouseholder& operator = (const EpetraHouseholder&);
00451
00452 protected:
00453
00455 enum PRECONDITIONER_METHOD {
00456 JACOBIAN,
00457 SMW
00458 };
00459
00461 Teuchos::RCP<LOCA::GlobalData> globalData;
00462
00464 Teuchos::RCP<Teuchos::ParameterList> solverParams;
00465
00467 Teuchos::RCP<LOCA::Epetra::Group> grp;
00468
00469
00470 Teuchos::RCP<const LOCA::BorderedSolver::AbstractOperator> op;
00471
00473 Teuchos::RCP<const NOX::Abstract::MultiVector> A;
00474
00476 Teuchos::RCP<const NOX::Abstract::MultiVector> B;
00477
00479 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> C;
00480
00482 Teuchos::RCP<const LOCA::MultiContinuation::ConstraintInterfaceMVDX> constraints;
00483
00485 LOCA::BorderedSolver::HouseholderQR qrFact;
00486
00488 Teuchos::RCP<NOX::Abstract::MultiVector> house_x;
00489
00491 NOX::Abstract::MultiVector::DenseMatrix house_p;
00492
00494 NOX::Abstract::MultiVector::DenseMatrix T;
00495
00497 NOX::Abstract::MultiVector::DenseMatrix R;
00498
00500 Teuchos::RCP<NOX::Abstract::MultiVector> U;
00501
00503 Teuchos::RCP<NOX::Abstract::MultiVector> V;
00504
00506 Teuchos::RCP<NOX::Abstract::MultiVector> house_x_trans;
00507
00509 NOX::Abstract::MultiVector::DenseMatrix house_p_trans;
00510
00512 NOX::Abstract::MultiVector::DenseMatrix T_trans;
00513
00515 NOX::Abstract::MultiVector::DenseMatrix R_trans;
00516
00518 Teuchos::RCP<NOX::Abstract::MultiVector> U_trans;
00519
00521 Teuchos::RCP<NOX::Abstract::MultiVector> V_trans;
00522
00524 Teuchos::RCP<const NOX::Abstract::MultiVector> Ablock;
00525
00527 Teuchos::RCP<const NOX::Abstract::MultiVector> Bblock;
00528
00530 Teuchos::RCP<NOX::Abstract::MultiVector> Ascaled;
00531
00533 Teuchos::RCP<NOX::Abstract::MultiVector> Bscaled;
00534
00536 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> Cscaled;
00537
00539 Teuchos::RCP<NOX::Epetra::LinearSystem> linSys;
00540
00542 Teuchos::RCP<Epetra_Operator> epetraOp;
00543
00545 Teuchos::RCP<const Epetra_BlockMap> baseMap;
00546
00548 Teuchos::RCP<const Epetra_BlockMap> globalMap;
00549
00551 int numConstraints;
00552
00554 bool isZeroA;
00555
00557 bool isZeroB;
00558
00560 bool isZeroC;
00561
00566 bool isValidForSolve;
00567
00572 bool isValidForTransposeSolve;
00573
00575 Teuchos::BLAS<int,double> dblas;
00576
00578 bool scale_rows;
00579
00581 std::vector<double> scale_vals;
00582
00584 PRECONDITIONER_METHOD precMethod;
00585
00587 bool includeUV;
00588
00590 bool use_P_For_Prec;
00591
00593 bool isComplex;
00594
00596 double omega;
00597
00598 };
00599 }
00600 }
00601
00602 #endif