Sundance 2.1

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 //  F.setName("F");
00050   LinearOperator<double> FInv = inverse(F, FSolver_);
00051 //  FInv.setName("FInv");
00052   LinearOperator<double> Bt = K.getBlock(0,1);
00053 //  Bt.setName("Bt");
00054 
00055 
00056   LinearOperator<double> Fp = FpProb_.getOperator();
00057 
00058   LinearOperator<double> Mp = MpProb_.getOperator();
00059 //  Mp.setName("Mp");
00060 
00061   LinearOperator<double> MpInv = inverse(Mp, MpSolver_);
00062 //  MpInv.setName("MpInv");
00063 
00064   LinearOperator<double> Ap = ApProb_.getOperator();
00065 //  Ap.setName("Ap");
00066 
00067   LinearOperator<double> ApInv = inverse(Ap, ApSolver_);
00068 //  ApInv.setName("ApInv");
00069 
00070 
00071   VectorSpace<double> pDomain = Bt.domain();
00072   VectorSpace<double> uDomain = F.domain();
00073 
00074   LinearOperator<double> Iu = identityOperator(uDomain);
00075 //  Iu.setName("Iu");
00076   LinearOperator<double> Ip = identityOperator(pDomain);
00077 //  Ip.setName("Ip");
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 //  Q1.setName("Q1");
00086   LinearOperator<double> Q2 = makeBlockOperator(colSpace, rowSpace);
00087   // Q2.setName("Q2");
00088   LinearOperator<double> Q3 = makeBlockOperator(colSpace, rowSpace);
00089   //Q3.setName("Q3");
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
SundancePCDPreconditioner.cpp
