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