Teko Version of the Day
Loading...
Searching...
No Matches
Teko_MLPreconditionerFactory.cpp
1#include "Teko_MLPreconditionerFactory.hpp"
2
3#include "Teko_MLLinearOp.hpp"
4
5#include "ml_include.h"
6#include "ml_MultiLevelPreconditioner.h"
7#include "ml_epetra_utils.h"
8#include "ml_RowMatrix.h"
9
10#include "Thyra_get_Epetra_Operator.hpp"
11
12namespace Teko {
13
15// MLPreconditionerState
17
18MLPreconditionerState::MLPreconditionerState() : isFilled_(false) {}
19
20void MLPreconditionerState::setMLComm(ML_Comm * comm)
21{
22 mlComm_ = Teuchos::rcpWithDealloc(comm,Teuchos::deallocFunctorDelete<ML_Comm>(cleanup_ML_Comm));
23}
24
25void MLPreconditionerState::setMLOperator(ML_Operator * op)
26{
27 mlOp_ = Teuchos::rcpWithDealloc(op,Teuchos::deallocFunctorDelete<ML_Operator>(cleanup_ML_Operator));
28}
29
30void MLPreconditionerState::setIsFilled(bool value)
31{
32 isFilled_ = value;
33}
34
35bool MLPreconditionerState::isFilled() const
36{
37 return isFilled_;
38}
39
40void MLPreconditionerState::cleanup_ML_Comm(ML_Comm * mlComm)
41{
42 ML_Comm_Destroy(&mlComm);
43}
44
45void MLPreconditionerState::cleanup_ML_Operator(ML_Operator * mlOp)
46{
47 ML_Operator_Destroy(&mlOp);
48}
49
50void MLPreconditionerState::setAggregationMatrices(const std::vector<Epetra_RowMatrix *> & diags)
51{
52 diagonalOps_ = diags;
53}
54
55Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner>
56MLPreconditionerState::constructMLPreconditioner(const Teuchos::ParameterList & mainList,
57 const std::vector<Teuchos::RCP<const Teuchos::ParameterList> > & coarseningParams)
58{
59 TEUCHOS_ASSERT(isFilled());
60 std::vector<Teuchos::ParameterList> cpls(coarseningParams.size());
61 for(std::size_t i=0;i<coarseningParams.size();i++)
62 cpls[i] = *coarseningParams[i];
63
64 mlPreconditioner_ = rcp(new ML_Epetra::MultiLevelPreconditioner(&*mlOp_, mainList, &diagonalOps_[0], &cpls[0], diagonalOps_.size()));
65
66 return mlPreconditioner_;
67}
68
70// MLPreconditionerFactory
72
73MLPreconditionerFactory::MLPreconditionerFactory()
74{
75}
76
77LinearOp MLPreconditionerFactory::buildPreconditionerOperator(BlockedLinearOp & blo,BlockPreconditionerState & state) const
78{
79 MLPreconditionerState & mlState = Teuchos::dyn_cast<MLPreconditionerState>(state);
80
81 if(not mlState.isFilled())
82 fillMLPreconditionerState(blo,mlState);
83 TEUCHOS_ASSERT(mlState.isFilled());
84
85 Teuchos::RCP<ML_Epetra::MultiLevelPreconditioner> mlPrecOp =
86 mlState.constructMLPreconditioner(mainParams_,coarseningParams_);
87
88 // return Thyra::epetraLinearOp(mlPrecOp);
89 return Teuchos::rcp(new MLLinearOp(mlPrecOp));
90}
91
92Teuchos::RCP<PreconditionerState> MLPreconditionerFactory::buildPreconditionerState() const
93{
94 return Teuchos::rcp(new MLPreconditionerState());
95}
96
97void MLPreconditionerFactory::fillMLPreconditionerState(const BlockedLinearOp & blo,MLPreconditionerState & mlState) const
98{
99 TEUCHOS_ASSERT(not mlState.isFilled());
100
101 EpetraExt::CrsMatrix_SolverMap RowMatrixColMapTrans;
102
103 int rowCnt = blockRowCount(blo);
104 int colCnt = blockRowCount(blo);
105
106 // construct comm object...add to state
107 ML_Comm * mlComm;
108 ML_Comm_Create(&mlComm);
109 mlState.setMLComm(mlComm);
110
111 ML_Operator * mlBlkMat = ML_Operator_Create(mlComm);
112 ML_Operator_BlkMatInit(mlBlkMat,mlComm,rowCnt,colCnt,ML_DESTROY_EVERYTHING);
113
114 std::vector<Epetra_RowMatrix *> aggMats;
115 for(int r=0;r<rowCnt;r++) {
116 for(int c=0;c<colCnt;c++) {
117 ML_Operator * tmp = 0;
118 Epetra_RowMatrix * EMat = 0;
119 std::stringstream ss;
120
121 ss << "fine_"<< r << "_" << c;
122
123 // extract crs matrix
124 Teuchos::RCP<Epetra_CrsMatrix> crsMat
125 = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(Thyra::get_Epetra_Operator(*blo->getNonconstBlock(r,c)));
126 EMat = ML_Epetra::ModifyEpetraMatrixColMap(*crsMat,RowMatrixColMapTrans, ss.str().c_str());
127 if(r==c) // setup diagonal scaling matrices
128 aggMats.push_back(EMat);
129
130 // extract ml sub operator
131 tmp = ML_Operator_Create(mlComm);
132 ML_Operator_WrapEpetraMatrix(EMat, tmp);
133
134 // add it to the block
135 ML_Operator_BlkMatInsert(mlBlkMat, tmp, r,c);
136 }
137 }
138 ML_Operator_BlkMatFinalize(mlBlkMat);
139
140 // finish setting up state object
141 mlState.setMLOperator(mlBlkMat);
142
143 // now set aggregation matrices
144 mlState.setAggregationMatrices(aggMats);
145
146 mlState.setIsFilled(true); // register state object as filled
147}
148
152void MLPreconditionerFactory::initializeFromParameterList(const Teuchos::ParameterList & settings)
153{
154 using Teuchos::RCP;
155 using Teuchos::rcp;
156
157 Teko_DEBUG_SCOPE("MLPreconditionerFactory::initializeFromParameterList",10);
158
159 // clear initial state
160 coarseningParams_.clear();
161 blockRowCount_ = 0;
162
163 blockRowCount_ = settings.get<int>("Block Row Count");
164
165 // read in main parameter list: with smoothing information
167 mainParams_ = settings.sublist("Smoothing Parameters");
168 mainParams_.set<Teuchos::RCP<const Teko::InverseLibrary> >("smoother: teko inverse library",getInverseLibrary());
169
170 // read in aggregation sub lists
172 const Teuchos::ParameterList & aggSublist = settings.sublist("Block Aggregation");
173
174 for(int block=0;block<blockRowCount_;block++) {
175 // write out sub list name: "Block #"
176 std::stringstream ss;
177 ss << "Block " << block;
178 std::string sublistName = ss.str();
179
180 // grab sublist
181 RCP<Teuchos::ParameterList> userSpec = rcp(new Teuchos::ParameterList(aggSublist.sublist(sublistName)));
182 coarseningParams_.push_back(userSpec);
183 }
184}
185
186} // end namespace Teko
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
An implementation of a state object for block preconditioners.
Contains operator internals need for ML.
bool isFilled() const
Has this object been filled yet.
void setAggregationMatrices(const std::vector< Epetra_RowMatrix * > &diags)
Set matrices to build multigrid hierarcy from.
void setIsFilled(bool value)
Set if ML internals been constructed yet.
void setMLOperator(ML_Operator *op)
set ML Operator pointer...it will be automatically cleaned up
void setMLComm(ML_Comm *comm)
set ML Comm pointer...it will be automatically cleaned up
Teuchos::RCP< ML_Epetra::MultiLevelPreconditioner > constructMLPreconditioner(const Teuchos::ParameterList &mainList, const std::vector< Teuchos::RCP< const Teuchos::ParameterList > > &coarseningParams)
Build an ML preconditioner object using the set of coarsening parameters.