Teko Version of the Day
Loading...
Searching...
No Matches
Teko_StratimikosFactory.cpp
1#include "Teko_StratimikosFactory.hpp"
2
3#include "Teuchos_Time.hpp"
4#include "Teuchos_AbstractFactoryStd.hpp"
5
6#include "Thyra_DefaultPreconditioner.hpp"
7
8#include "Teko_InverseLibrary.hpp"
9#include "Teko_Preconditioner.hpp"
10#include "Teko_Utilities.hpp"
11#include "Teko_InverseLibrary.hpp"
12#include "Teko_ReorderedLinearOp.hpp"
13
14#ifdef TEKO_HAVE_EPETRA
15#include "Teko_InverseFactoryOperator.hpp" // an epetra specific object
16#include "Thyra_EpetraOperatorViewExtractorStd.hpp"
17#include "Teko_StridedEpetraOperator.hpp"
18#include "Teko_BlockedEpetraOperator.hpp"
19#include "Thyra_EpetraLinearOp.hpp"
20#include "EpetraExt_RowMatrixOut.h"
21#endif
22
23namespace Teko {
24
25using Teuchos::RCP;
26using Teuchos::ParameterList;
27
28// hide stuff
29namespace {
30 // Simple preconditioner class that adds a counter
31 class StratimikosFactoryPreconditioner : public Thyra::DefaultPreconditioner<double>{
32 public:
33 StratimikosFactoryPreconditioner() : iter_(0) {}
34
35 inline void incrIter() { iter_++; }
36 inline std::size_t getIter() { return iter_; }
37
38 private:
39 StratimikosFactoryPreconditioner(const StratimikosFactoryPreconditioner &);
40
41 std::size_t iter_;
42 };
43
44 // factory used to initialize the Teko::StratimikosFactory
45 // user data
46 class TekoFactoryBuilder
47 : public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
48 public:
49 TekoFactoryBuilder(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
50 const Teuchos::RCP<Teko::RequestHandler> & rh) : builder_(builder), requestHandler_(rh) {}
51 Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create() const
52 { return Teuchos::rcp(new StratimikosFactory(builder_,requestHandler_)); }
53
54 private:
55 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builder_;
56 Teuchos::RCP<Teko::RequestHandler> requestHandler_;
57 };
58}
59
60#ifdef TEKO_HAVE_EPETRA
61
62// Constructors/initializers/accessors
64 :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
65{}
66
67// Constructors/initializers/accessors
68StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Teko::RequestHandler> & rh)
69 :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
70{
71 setRequestHandler(rh);
72}
73
74#endif // TEKO_HAVE_EPETRA
75
76StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
77 const Teuchos::RCP<Teko::RequestHandler> & rh)
78 :
79#ifdef TEKO_HAVE_EPETRA
80 epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())),
81#endif
82 builder_(builder)
83{
84 setRequestHandler(rh);
85}
86
87// Overridden from PreconditionerFactoryBase
88bool StratimikosFactory::isCompatible(const Thyra::LinearOpSourceBase<double> &/* fwdOpSrc */) const
89{
90 using Teuchos::outArg;
91
92 TEUCHOS_ASSERT(false); // what you doing?
93
94/*
95 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
96 Thyra::EOpTransp epetraFwdOpTransp;
97 Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
98 Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
99 double epetraFwdOpScalar;
100
101 // check to make sure this is an epetra CrsMatrix
102 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc.getOp();
103 epetraFwdOpViewExtractor_->getEpetraOpView(
104 fwdOp,
105 outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
106 outArg(epetraFwdOpApplyAs),
107 outArg(epetraFwdOpAdjointSupport),
108 outArg(epetraFwdOpScalar));
109
110 if( !dynamic_cast<const Epetra_CrsMatrix*>(&*epetraFwdOp) )
111 return false;
112*/
113
114 return true;
115}
116
117
118bool StratimikosFactory::applySupportsConj(Thyra::EConj /* conj */) const
119{
120 return false;
121}
122
123
124bool StratimikosFactory::applyTransposeSupportsConj(Thyra::EConj /* conj */) const
125{
126 return false; // See comment below
127}
128
129
130Teuchos::RCP<Thyra::PreconditionerBase<double> >
131StratimikosFactory::createPrec() const
132{
133 return Teuchos::rcp(new StratimikosFactoryPreconditioner());
134}
135
136
137void StratimikosFactory::initializePrec(
138 const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
139 Thyra::PreconditionerBase<double> *prec,
140 const Thyra::ESupportSolveUse supportSolveUse
141 ) const
142{
143 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
144
145#ifdef TEKO_HAVE_EPETRA
146 if(epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
147 initializePrec_Epetra(fwdOpSrc,prec,supportSolveUse);
148 else
149#endif
150 initializePrec_Thyra(fwdOpSrc,prec,supportSolveUse);
151}
152
153void StratimikosFactory::initializePrec_Thyra(
154 const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
155 Thyra::PreconditionerBase<double> *prec,
156 const Thyra::ESupportSolveUse /* supportSolveUse */
157 ) const
158{
159 using Teuchos::RCP;
160 using Teuchos::rcp;
161 using Teuchos::rcp_dynamic_cast;
162 using Teuchos::implicit_cast;
163 using Thyra::LinearOpBase;
164
165 Teuchos::Time totalTimer(""), timer("");
166 totalTimer.start(true);
167
168 const RCP<Teuchos::FancyOStream> out = this->getOStream();
169 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
170
171 bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
172
173 Teuchos::OSTab tab(out);
174 if(mediumVerbosity)
175 *out << "\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
176
177 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
178
179 // Get the concrete preconditioner object
180 StratimikosFactoryPreconditioner & defaultPrec
181 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
182 Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
183
184 // Permform initialization if needed
185 const bool startingOver = (prec_Op == Teuchos::null);
186 if(startingOver) {
187 // start over
188 invLib_ = Teuchos::null;
189 invFactory_ = Teuchos::null;
190
191 if(mediumVerbosity)
192 *out << "\nCreating the initial Teko Operator object...\n";
193
194 timer.start(true);
195
196 // build library, and set request handler (user defined!)
197 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"),builder_);
198 invLib_->setRequestHandler(reqHandler_);
199
200 // build preconditioner factory
201 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
202
203 timer.stop();
204 if(mediumVerbosity)
205 Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
206 }
207
208 if(mediumVerbosity)
209 *out << "\nComputing the preconditioner ...\n";
210
211 // setup reordering if required
212 std::string reorderType = paramList_->get<std::string>("Reorder Type");
213 if(reorderType!="") {
214
215 Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
216 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(fwdOp,true);
217 RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
218 Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm,blkFwdOp);
219
220 if(prec_Op==Teuchos::null) {
221 Teko::ModifiableLinearOp reorderedPrec = Teko::buildInverse(*invFactory_,blockedFwdOp);
222 prec_Op = Teuchos::rcp(new ReorderedLinearOp(brm,reorderedPrec));
223 }
224 else {
225 Teko::ModifiableLinearOp reorderedPrec = Teuchos::rcp_dynamic_cast<ReorderedLinearOp>(prec_Op,true)->getBlockedOp();
226 Teko::rebuildInverse(*invFactory_,blockedFwdOp,reorderedPrec);
227 }
228 }
229 else {
230 // no reordering required
231 if(prec_Op==Teuchos::null)
232 prec_Op = Teko::buildInverse(*invFactory_,fwdOp);
233 else
234 Teko::rebuildInverse(*invFactory_,fwdOp,prec_Op);
235 }
236
237 // construct preconditioner
238 timer.stop();
239
240 if(mediumVerbosity)
241 Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
242
243
244 // Initialize the preconditioner
245 defaultPrec.initializeUnspecified( prec_Op);
246 defaultPrec.incrIter();
247 totalTimer.stop();
248
249 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
250 *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
251 if(mediumVerbosity)
252 *out << "\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
253}
254
255#ifdef TEKO_HAVE_EPETRA
256void StratimikosFactory::initializePrec_Epetra(
257 const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
258 Thyra::PreconditionerBase<double> *prec,
259 const Thyra::ESupportSolveUse /* supportSolveUse */
260 ) const
261{
262 using Teuchos::outArg;
263 using Teuchos::RCP;
264 using Teuchos::rcp;
265 using Teuchos::rcp_dynamic_cast;
266 using Teuchos::implicit_cast;
267
268 Teuchos::Time totalTimer(""), timer("");
269 totalTimer.start(true);
270
271 const RCP<Teuchos::FancyOStream> out = this->getOStream();
272 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
273
274 bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
275
276 Teuchos::OSTab tab(out);
277 if(mediumVerbosity)
278 *out << "\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
279
280 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
281#ifdef _DEBUG
282 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
283 TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
284#endif
285
286 //
287 // Unwrap and get the forward Epetra_Operator object
288 //
289 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
290 Thyra::EOpTransp epetraFwdOpTransp;
291 Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
292 Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
293 double epetraFwdOpScalar;
294 epetraFwdOpViewExtractor_->getEpetraOpView(
295 fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
296 outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
297 );
298 // Get the concrete preconditioner object
299 StratimikosFactoryPreconditioner & defaultPrec
300 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
301
302 // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
303 RCP<Thyra::EpetraLinearOp> epetra_precOp
304 = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),true);
305
306 // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
307 Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
308 if(epetra_precOp!=Teuchos::null)
309 teko_precOp = rcp_dynamic_cast<Teko::Epetra::InverseFactoryOperator>(epetra_precOp->epetra_op(),true);
310
311 // Permform initialization if needed
312 const bool startingOver = (teko_precOp == Teuchos::null);
313 if(startingOver) {
314 // start over
315 invLib_ = Teuchos::null;
316 invFactory_ = Teuchos::null;
317 decomp_.clear();
318
319 if(mediumVerbosity)
320 *out << "\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
321
322 timer.start(true);
323
324 std::stringstream ss;
325 ss << paramList_->get<std::string>("Strided Blocking");
326
327 // figure out the decomposition requested by the string
328 while(not ss.eof()) {
329 int num=0;
330 ss >> num;
331
332 TEUCHOS_ASSERT(num>0);
333
334 decomp_.push_back(num);
335 }
336
337 // build library, and set request handler (user defined!)
338 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"),builder_);
339 invLib_->setRequestHandler(reqHandler_);
340
341 // build preconditioner factory
342 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
343
344 // Create the initial preconditioner: DO NOT compute it yet
345 teko_precOp = rcp( new Teko::Epetra::InverseFactoryOperator(invFactory_));
346
347 timer.stop();
348 if(mediumVerbosity)
349 Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
350 }
351
352 if(mediumVerbosity)
353 *out << "\nComputing the preconditioner ...\n";
354
355 // construct preconditioner
356 {
357 bool writeBlockOps = paramList_->get<bool>("Write Block Operator");
358 timer.start(true);
359 teko_precOp->initInverse();
360
361 if(decomp_.size()==1) {
362 teko_precOp->rebuildInverseOperator(epetraFwdOp);
363
364 // write out to disk
365 if(writeBlockOps) {
366 std::stringstream ss;
367 ss << "block-" << defaultPrec.getIter() << "_00.mm";
368 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*rcp_dynamic_cast<const Epetra_CrsMatrix>(epetraFwdOp));
369 }
370 }
371 else {
372 Teuchos::RCP<Epetra_Operator> wrappedFwdOp
373 = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
374
375 // write out to disk
376 if(writeBlockOps) {
377 std::stringstream ss;
378 ss << "block-" << defaultPrec.getIter();
379 // Teuchos::RCP<Teko::Epetra::StridedEpetraOperator> stridedJac
380 // = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedFwdOp);
381 Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac
382 = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedFwdOp);
383 if(stridedJac!=Teuchos::null) {
384 // write out blocks of strided operator
385 stridedJac->WriteBlocks(ss.str());
386 }
387 else TEUCHOS_ASSERT(false);
388 }
389
390 teko_precOp->rebuildInverseOperator(wrappedFwdOp);
391 }
392
393 timer.stop();
394 }
395
396 if(mediumVerbosity)
397 Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
398
399 // Initialize the output EpetraLinearOp
400 if(startingOver) {
401 epetra_precOp = rcp(new Thyra::EpetraLinearOp);
402 }
403 epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
404 ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
405 ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
406
407 // Initialize the preconditioner
408 defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
409 defaultPrec.incrIter();
410 totalTimer.stop();
411
412 if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
413 *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
414 if(mediumVerbosity)
415 *out << "\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
416}
417#endif // TEKO_HAVE_EPETRA
418
419void StratimikosFactory::uninitializePrec(
420 Thyra::PreconditionerBase<double> * /* prec */,
421 Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > * /* fwdOp */,
422 Thyra::ESupportSolveUse * /* supportSolveUse */
423 ) const
424{
425 TEUCHOS_TEST_FOR_EXCEPT(true);
426}
427
428
429// Overridden from ParameterListAcceptor
430
431
432void StratimikosFactory::setParameterList(
433 Teuchos::RCP<Teuchos::ParameterList> const& paramList
434 )
435{
436 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
437
438 paramList->validateParametersAndSetDefaults(*this->getValidParameters(),0);
439 paramList_ = paramList;
440
441}
442
443
444Teuchos::RCP<Teuchos::ParameterList>
445StratimikosFactory::getNonconstParameterList()
446{
447 return paramList_;
448}
449
450
451Teuchos::RCP<Teuchos::ParameterList>
452StratimikosFactory::unsetParameterList()
453{
454 Teuchos::RCP<ParameterList> _paramList = paramList_;
455 paramList_ = Teuchos::null;
456 return _paramList;
457}
458
459
460Teuchos::RCP<const Teuchos::ParameterList>
461StratimikosFactory::getParameterList() const
462{
463 return paramList_;
464}
465
466
467Teuchos::RCP<const Teuchos::ParameterList>
468StratimikosFactory::getValidParameters() const
469{
470
471 using Teuchos::rcp;
472 using Teuchos::tuple;
473 using Teuchos::implicit_cast;
474 using Teuchos::rcp_implicit_cast;
475
476 static RCP<const ParameterList> validPL;
477
478 if(is_null(validPL)) {
479
480 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
481
482 pl->set("Test Block Operator",false,
483 "If Stratiikos/Teko is used to break an operator into its parts,\n"
484 "then setting this parameter to true will compare applications of the\n"
485 "segregated operator to the original operator.");
486 pl->set("Write Block Operator",false,
487 "Write out the segregated operator to disk with the name \"block-?_xx\"");
488 pl->set("Strided Blocking","1",
489 "Assuming that the user wants Strided blocking, break the operator into\n"
490 "blocks. The syntax can be thought to be associated with the solution\n"
491 "vector. For example if your variables are [u v w p T], and we want [u v w]\n"
492 "blocked together, and p and T separate then the relevant string is \"3 1 1\".\n"
493 "Meaning put the first 3 unknowns per node together and separate the v and w\n"
494 "components.");
495 pl->set("Reorder Type","",
496 "This specifies how the blocks are reordered for use in the preconditioner.\n"
497 "For example, assume the linear system is generated from 3D Navier-Stokes\n"
498 "with an energy equation, yielding the unknowns [u v w p T]. If the\n"
499 "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n"
500 "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n"
501 "velocity and pressure forming an inner two-by-two block, and then the\n"
502 "temperature unknowns forming a two-by-two system with the velocity-pressure\n"
503 "block.");
504 pl->set("Inverse Type","Amesos",
505 "The type of inverse operator the user wants. This can be one of the defaults\n"
506 "from Stratimikos, or a Teko preconditioner defined in the\n"
507 "\"Inverse Factory Library\".");
508 pl->sublist("Inverse Factory Library",false,"Definition of Teko preconditioners.");
509
510 validPL = pl;
511 }
512
513 return validPL;
514}
515
516#ifdef TEKO_HAVE_EPETRA
520Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
521 const Teuchos::RCP<const Epetra_Operator> & Jac,
522 const Teuchos::RCP<Epetra_Operator> & wrapInput,
523 std::ostream & out
524 ) const
525{
526 Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
527
528// // initialize jacobian
529// if(wrappedOp==Teuchos::null)
530// {
531// wrappedOp = Teuchos::rcp(new Teko::Epetra::StridedEpetraOperator(decomp_,Jac));
532//
533// // reorder the blocks if requested
534// std::string reorderType = paramList_->get<std::string>("Reorder Type");
535// if(reorderType!="") {
536// RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
537//
538// // out << "Teko: Reordering = " << brm->toString() << std::endl;
539// Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->Reorder(*brm);
540// }
541// }
542// else {
543// Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->RebuildOps();
544// }
545//
546// // test blocked operator for correctness
547// if(paramList_->get<bool>("Test Block Operator")) {
548// bool result
549// = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
550//
551// out << "Teko: Tested operator correctness: " << (result ? "passed" : "FAILED!") << std::endl;
552// }
553
554 // initialize jacobian
555 if(wrappedOp==Teuchos::null)
556 {
557 // build strided vector
558 std::vector<std::vector<int> > vars;
559 buildStridedVectors(*Jac,decomp_,vars);
560 wrappedOp = Teuchos::rcp(new Teko::Epetra::BlockedEpetraOperator(vars,Jac));
561
562 // reorder the blocks if requested
563 std::string reorderType = paramList_->get<std::string>("Reorder Type");
564 if(reorderType!="") {
565 RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
566
567 // out << "Teko: Reordering = " << brm->toString() << std::endl;
568 Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->Reorder(*brm);
569 }
570 }
571 else {
572 Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->RebuildOps();
573 }
574
575 // test blocked operator for correctness
576 if(paramList_->get<bool>("Test Block Operator")) {
577 bool result
578 = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
579
580 out << "Teko: Tested operator correctness: " << (result ? "passed" : "FAILED!") << std::endl;
581 }
582
583 return wrappedOp;
584}
585#endif // TEKO_HAVE_EPETRA
586
587std::string StratimikosFactory::description() const
588{
589 std::ostringstream oss;
590 oss << "Teko::StratimikosFactory";
591 return oss.str();
592}
593
594#ifdef TEKO_HAVE_EPETRA
595void StratimikosFactory::buildStridedVectors(const Epetra_Operator & Jac,
596 const std::vector<int> & decomp,
597 std::vector<std::vector<int> > & vars) const
598{
599 const Epetra_Map & rangeMap = Jac.OperatorRangeMap();
600
601 // compute total number of variables
602 int numVars = 0;
603 for(std::size_t i=0;i<decomp.size();i++)
604 numVars += decomp[i];
605
606 // verify that the decomposition is appropriate for this matrix
607 TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars)==0);
608 TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars)==0);
609
610 int * globalIds = rangeMap.MyGlobalElements();
611
612 vars.resize(decomp.size());
613 for(int i=0;i<rangeMap.NumMyElements();) {
614
615 // for each "node" copy global ids to vectors
616 for(std::size_t d=0;d<decomp.size();d++) {
617 // for this variable copy global ids to variable arrays
618 int current = decomp[d];
619 for(int v=0;v<current;v++,i++)
620 vars[d].push_back(globalIds[i]);
621 }
622 }
623}
624#endif // TEKO_HAVE_EPETRA
625
626void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
627 const std::string & stratName)
628{
629 TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
630 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
631
632 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
633 = Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder(builder));
634
635 // use default constructor to add Teko::StratimikosFactory
636 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(builderCopy,Teuchos::null));
637 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
638}
639
640void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
641 const Teuchos::RCP<Teko::RequestHandler> & rh,
642 const std::string & stratName)
643{
644 TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
645 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
646
647 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
648 = Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder(builder));
649
650 // build an instance of a Teuchos::AbsractFactory<Thyra::PFB> so request handler is passed onto
651 // the resulting StratimikosFactory
652 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(builderCopy,rh));
653 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
654}
655
656} // namespace Thyra
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
Tear about a user specified Epetra_Operator (CrsMatrix) using a vector of vectors of GIDs for each bl...
A single Epetra wrapper for all operators constructed from an inverse operator.
This class takes a blocked linear op and represents it in a flattened form.