Go to the documentation of this file.00001 #include "SundancePCDPreconditioner.hpp"
00002 #include "Sundance.hpp"
00003 #include "TSFLinearSolverBuilder.hpp"
00004 #include "TSFLinearSolverImpl.hpp"
00005 #include "TSFSimpleIdentityOpDecl.hpp"
00006 #include "TSFSimpleComposedOpDecl.hpp"
00007 #include "TSFSimpleBlockOpDecl.hpp"
00008 #include "TSFSimpleScaledOpDecl.hpp"
00009 #include "TSFGenericRightPreconditioner.hpp"
00010 #include "TSFLinearCombinationImpl.hpp"
00011
00012 using namespace TSFExtended;
00013 using namespace Sundance;
00014
00015 PCDPreconditionerFactory::PCDPreconditionerFactory(
00016 const ParameterList& params,
00017 const LinearProblem& MpProb,
00018 const LinearProblem& ApProb,
00019 const LinearProblem& FpProb
00020 )
00021 : MpProb_(MpProb),
00022 ApProb_(ApProb),
00023 FpProb_(FpProb),
00024 MpSolver_(),
00025 ApSolver_(),
00026 FSolver_()
00027 {
00028 ParameterList msParams = params.sublist("MpSolver");
00029 MpSolver_ = LinearSolverBuilder::createSolver(msParams);
00030 ParameterList asParams = params.sublist("ApSolver");
00031 ApSolver_ = LinearSolverBuilder::createSolver(asParams);
00032 ParameterList fsParams = params.sublist("FSolver");
00033 FSolver_ = LinearSolverBuilder::createSolver(fsParams);
00034 }
00035
00036 Preconditioner<double>
00037 PCDPreconditionerFactory::
00038 createPreconditioner(const LinearOperator<double>& K) const
00039 {
00040 Tabs tab;
00041
00042 LinearOperator<double> F = K.getBlock(0,0);
00043 F.setName("F");
00044 LinearOperator<double> FInv = inverse(F, FSolver_);
00045 FInv.setName("FInv");
00046 LinearOperator<double> Bt = K.getBlock(0,1);
00047 Bt.setName("Bt");
00048
00049
00050 LinearOperator<double> Fp = FpProb_.getOperator();
00051
00052 LinearOperator<double> Mp = MpProb_.getOperator();
00053 Mp.setName("Mp");
00054
00055 LinearOperator<double> MpInv = inverse(Mp, MpSolver_);
00056 MpInv.setName("MpInv");
00057
00058 LinearOperator<double> Ap = ApProb_.getOperator();
00059 Ap.setName("Ap");
00060
00061 LinearOperator<double> ApInv = inverse(Ap, ApSolver_);
00062 ApInv.setName("ApInv");
00063
00064
00065 VectorSpace<double> pDomain = Bt.domain();
00066 VectorSpace<double> uDomain = F.domain();
00067
00068 LinearOperator<double> Iu = identityOperator(uDomain);
00069 Iu.setName("Iu");
00070 LinearOperator<double> Ip = identityOperator(pDomain);
00071 Ip.setName("Ip");
00072
00073 LinearOperator<double> XInv = MpInv * Fp * ApInv;
00074
00075 VectorSpace<double> rowSpace = K.range();
00076 VectorSpace<double> colSpace = K.domain();
00077
00078 LinearOperator<double> Q1 = makeBlockOperator(colSpace, rowSpace);
00079 Q1.setName("Q1");
00080 LinearOperator<double> Q2 = makeBlockOperator(colSpace, rowSpace);
00081 Q2.setName("Q2");
00082 LinearOperator<double> Q3 = makeBlockOperator(colSpace, rowSpace);
00083 Q3.setName("Q3");
00084
00085 Q1.setBlock(0, 0, FInv);
00086 Q1.setBlock(1, 1, Ip);
00087 Q1.endBlockFill();
00088
00089 Q2.setBlock(0, 0, Iu);
00090 Q2.setBlock(0, 1, -1.0*Bt);
00091 Q2.setBlock(1, 1, Ip);
00092 Q2.endBlockFill();
00093
00094 Q3.setBlock(0, 0, Iu);
00095 Q3.setBlock(1, 1, -1.0*XInv);
00096 Q3.endBlockFill();
00097
00098 LinearOperator<double> P1 = Q2 * Q3;
00099 LinearOperator<double> PInv = Q1 * P1;
00100
00101 return new GenericRightPreconditioner<double>(PInv);
00102 }
00103
00104
