Teko Version of the Day
Loading...
Searching...
No Matches
Teko_mlutils.cpp
1#include "Teko_mlutils.hpp"
2
3#include <vector>
4
5#ifdef HAVE_MPI
6 #include "mpi.h"
7 #include "Epetra_MpiComm.h"
8#else
9 #include "Epetra_SerialComm.h"
10#endif
11#include "Epetra_Vector.h"
12
13#include "ml_epetra_utils.h"
14#include "ml_op_utils.h"
15
16#include "Thyra_EpetraLinearOp.hpp"
17
18#include "EpetraExt_RowMatrixOut.h"
19
20#include "Teuchos_XMLParameterListHelpers.hpp"
21
22#include "Teko_InverseFactory.hpp"
23#include "Teko_InverseLibrary.hpp"
24#include "Teko_Utilities.hpp"
25#include "Teko_PreconditionerFactory.hpp"
26#include "Teko_EpetraOperatorWrapper.hpp"
27#include "Teko_SmootherPreconditionerFactory.hpp"
28
29namespace Teko {
30namespace mlutils {
31
32// build a very simple row map from the ML_Operator
33Teuchos::RCP<Epetra_Map> buildRowMap(ML_Operator * mlOp)
34{
35 ML_Comm *comm = mlOp->comm;
36#ifdef ML_MPI
37 MPI_Comm mpi_comm;
38 mpi_comm = comm->USR_comm;
39 Epetra_MpiComm eComm( mpi_comm ) ;
40#else
41 Epetra_SerialComm eComm ;
42#endif
43
44 // get operators row count, build map...who cares what GIDs are
45 int rowCnt = mlOp->outvec_leng;
46 return Teuchos::rcp(new Epetra_Map(-1,rowCnt,0,eComm));
47}
48
52Teuchos::RCP<Epetra_CrsMatrix> convertToCrsMatrix(ML_Operator * mlOp,
53 const Teuchos::RCP<Epetra_Map> & rowMapArg)
54{
55 ML_Comm *comm = mlOp->comm;
56#ifdef ML_MPI
57 MPI_Comm mpi_comm;
58 mpi_comm = comm->USR_comm;
59 Epetra_MpiComm eComm( mpi_comm ) ;
60#else
61 Epetra_SerialComm eComm ;
62#endif
63
64 // build a default row map
65 Teuchos::RCP<Epetra_Map> rowMap = rowMapArg;
66 if(rowMapArg==Teuchos::null)
67 rowMap = buildRowMap(mlOp);
68
69 // build lightweight CrsMatrix wrapper
70 Epetra_CrsMatrix * crsMatrixPtr = 0;
71 int maxNumNonzeros = 0;
72 double cpuTime = 0.0;
73 bool verbose = false;
74 ML_Operator2EpetraCrsMatrix(mlOp,crsMatrixPtr,maxNumNonzeros,false,cpuTime,verbose);
75
76 return Teuchos::rcp(crsMatrixPtr);
77}
78
79Teko::LinearOp buildTekoBlockOp(ML_Operator * mlOp,int level)
80{
81 Teko_DEBUG_SCOPE("Teko::mlutils::buildTekoBlockOp",0);
82
83 using Teuchos::RCP;
84
85 int numRows = ML_Operator_BlkMatNumBlockRows(mlOp);
86 int numCols = ML_Operator_BlkMatNumBlockCols(mlOp);
87
88 Teko::BlockedLinearOp tekoOp = Teko::createBlockedOp();
89 Teko::beginBlockFill(tekoOp,numRows,numCols);
90 for(int i=0;i<numRows;i++) {
91 for(int j=0;j<numCols;j++) {
92 ML_Operator * subBlock = ML_Operator_BlkMatExtract(mlOp,i,j);
93
94 // if required construct and add a block to tekoOp
95 if(subBlock!=0) {
96 // create a CRS matrix from ML operator
97 RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(subBlock);
98
99 #if 0
100 std::stringstream ss;
101 ss << "A(" << level << ")_" << i << "_" << j << ".mm";
102 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*eCrsMat);
103 #endif
104
105 // build Teko operator, add to Teko blocked operator
106 Teko::LinearOp tekoSubBlock = Thyra::epetraLinearOp(eCrsMat);
107 Teko::setBlock(i,j,tekoOp,tekoSubBlock);
108 }
109 }
110 }
111 Teko::endBlockFill(tekoOp);
112
113 return tekoOp.getConst();
114}
115
116int smoother(ML_Smoother *mydata, int leng1, double x[], int leng2,
117 double rhs[])
118{
119 Teko_DEBUG_SCOPE("Teko::mlutils::smoother",10);
120
121 // std::cout << "init guess = " << mydata->init_guess << std::endl;
122
123 // grab data object
124 SmootherData * smootherData = (struct SmootherData *) ML_Get_MySmootherData(mydata);
125
126 Epetra_Vector X(View, smootherData->Amat->OperatorDomainMap(), x);
127 Epetra_Vector Y(View, smootherData->Amat->OperatorRangeMap(),rhs);
128
129 smootherData->smootherOperator->Apply(Y,X);
130
131 return 0;
132}
133
134Teuchos::RCP<Teko::InverseFactory> ML_Gen_Init_Teko_Prec(const std::string smoothName,const Teko::InverseLibrary & invLib)
135{
136 Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Init_Teko_Prec",10);
137
138 // Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
139 // Teuchos::RCP<Teko::InverseLibrary> invLib
140 // = Teko::InverseLibrary::buildFromParameterList(pl);
141 // invLib->setRequestHandler(rh);
142
143 Teuchos::RCP<Teko::PreconditionerFactory> precFact;
144 Teuchos::RCP<Teko::InverseFactory> invFact = invLib.getInverseFactory(smoothName);
145
146 return invFact;
147}
148
149extern "C"
150void ML_Destroy_Smoother_Teko(void * data)
151{
153 delete sData;
154}
155
156extern "C"
157int ML_Gen_Smoother_Teko(ML *ml, int level, int pre_or_post, int ntimes,const Teuchos::RCP<const Teuchos::ParameterList> & tekoPL,
158 const Teuchos::RCP<const Teko::InverseLibrary> & invLib_in,const std::string & inverse, bool isBlocked)
159{
160 Teko_DEBUG_SCOPE("Teko::mlutils::ML_Gen_Smoother_Teko",10);
161 ML_Operator * BlkMat = &(ml->Amat[level]);
162
163 Teuchos::RCP<const Teko::InverseLibrary> invLib;
164 if(invLib_in ==Teuchos::null) {
165 Teuchos::RCP<Teko::RequestHandler> rh = Teuchos::rcp(new Teko::RequestHandler());
166 Teuchos::RCP<Teko::InverseLibrary> ncInvLib = Teko::InverseLibrary::buildFromParameterList(*tekoPL);
167 ncInvLib->setRequestHandler(rh);
168
169 invLib = ncInvLib.getConst();
170 }
171 else {
172 // this assumes a request handler has already been set
173 invLib = invLib_in;
174 }
175
176 Teko::LinearOp tekoAmat;
177 if(isBlocked) {
178 // build teko blocked operator
179 tekoAmat = Teko::mlutils::buildTekoBlockOp(BlkMat,level);
180 }
181 else {
182 // build unblocked operator
183 Teuchos::RCP<const Epetra_CrsMatrix> eCrsMat = convertToCrsMatrix(BlkMat);
184 tekoAmat = Thyra::epetraLinearOp(eCrsMat);
185 }
186
187 // Teuchos::RCP<Teko::mlutils::SmootherData> data = Teuchos::rcp(new Teko::mlutils::SmootherData);
189 data->Amat = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(tekoAmat));
190
191 // Teuchos::ParameterList pl = *Teuchos::getParametersFromXmlFile(filename);
192 Teuchos::RCP<Teko::InverseFactory> invFact = ML_Gen_Init_Teko_Prec(inverse,*invLib);
193 Teuchos::RCP<Teko::RequestHandler> rh = invFact->getRequestHandler();
194
195 // build smoother operator
196 Teko::LinearOp precOp = Teko::buildInverse(*invFact,tekoAmat);
197 Teko::LinearOp smootherOp = Teko::buildSmootherLinearOp(tekoAmat,precOp,ntimes,true);
198
199 // get an epetra operator wrapper
200 data->smootherOperator = Teuchos::rcp(new Teko::Epetra::EpetraOperatorWrapper(smootherOp));
201
202 int ret_val = ML_Set_Smoother(ml,level, pre_or_post, (void*) data, Teko::mlutils::smoother, NULL);
203 ml->post_smoother[level].data_destroy = ML_Destroy_Smoother_Teko;
204
205 return ret_val;
206}
207
208}
209}
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void beginBlockFill(BlockedLinearOp &blo, int rowCnt, int colCnt)
Let the blocked operator know that you are going to set the sub blocks.
void setBlock(int i, int j, BlockedLinearOp &blo, const LinearOp &lo)
Set the i,j block in a BlockedLinearOp object.
BlockedLinearOp createBlockedOp()
Build a new blocked linear operator.
void endBlockFill(BlockedLinearOp &blo)
Notify the blocked operator that the fill stage is completed.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...