1#include "Teko_StratimikosFactory.hpp"
3#include "Teuchos_Time.hpp"
4#include "Teuchos_AbstractFactoryStd.hpp"
6#include "Thyra_DefaultPreconditioner.hpp"
8#include "Teko_InverseLibrary.hpp"
9#include "Teko_Preconditioner.hpp"
11#include "Teko_InverseLibrary.hpp"
12#include "Teko_ReorderedLinearOp.hpp"
14#ifdef TEKO_HAVE_EPETRA
15#include "Teko_InverseFactoryOperator.hpp"
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"
26using Teuchos::ParameterList;
31 class StratimikosFactoryPreconditioner :
public Thyra::DefaultPreconditioner<double>{
33 StratimikosFactoryPreconditioner() : iter_(0) {}
35 inline void incrIter() { iter_++; }
36 inline std::size_t getIter() {
return iter_; }
39 StratimikosFactoryPreconditioner(
const StratimikosFactoryPreconditioner &);
46 class TekoFactoryBuilder
47 :
public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
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_)); }
55 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builder_;
56 Teuchos::RCP<Teko::RequestHandler> requestHandler_;
60#ifdef TEKO_HAVE_EPETRA
64 :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
68StratimikosFactory::StratimikosFactory(
const Teuchos::RCP<Teko::RequestHandler> & rh)
69 :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
71 setRequestHandler(rh);
76StratimikosFactory::StratimikosFactory(
const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
77 const Teuchos::RCP<Teko::RequestHandler> & rh)
79#ifdef TEKO_HAVE_EPETRA
80 epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())),
84 setRequestHandler(rh);
88bool StratimikosFactory::isCompatible(
const Thyra::LinearOpSourceBase<double> &)
const
90 using Teuchos::outArg;
92 TEUCHOS_ASSERT(
false);
118bool StratimikosFactory::applySupportsConj(Thyra::EConj )
const
124bool StratimikosFactory::applyTransposeSupportsConj(Thyra::EConj )
const
130Teuchos::RCP<Thyra::PreconditionerBase<double> >
131StratimikosFactory::createPrec()
const
133 return Teuchos::rcp(
new StratimikosFactoryPreconditioner());
137void StratimikosFactory::initializePrec(
138 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
139 Thyra::PreconditionerBase<double> *prec,
140 const Thyra::ESupportSolveUse supportSolveUse
143 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
145#ifdef TEKO_HAVE_EPETRA
146 if(epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
147 initializePrec_Epetra(fwdOpSrc,prec,supportSolveUse);
150 initializePrec_Thyra(fwdOpSrc,prec,supportSolveUse);
153void StratimikosFactory::initializePrec_Thyra(
154 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
155 Thyra::PreconditionerBase<double> *prec,
156 const Thyra::ESupportSolveUse
161 using Teuchos::rcp_dynamic_cast;
162 using Teuchos::implicit_cast;
163 using Thyra::LinearOpBase;
165 Teuchos::Time totalTimer(
""), timer(
"");
166 totalTimer.start(
true);
168 const RCP<Teuchos::FancyOStream> out = this->getOStream();
169 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
171 bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
173 Teuchos::OSTab tab(out);
175 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
177 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
180 StratimikosFactoryPreconditioner & defaultPrec
181 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
182 Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
185 const bool startingOver = (prec_Op == Teuchos::null);
188 invLib_ = Teuchos::null;
189 invFactory_ = Teuchos::null;
192 *out <<
"\nCreating the initial Teko Operator object...\n";
197 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist(
"Inverse Factory Library"),builder_);
198 invLib_->setRequestHandler(reqHandler_);
201 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
205 Teuchos::OSTab(out).o() <<
"> Creation time = "<<timer.totalElapsedTime()<<
" sec\n";
209 *out <<
"\nComputing the preconditioner ...\n";
212 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
213 if(reorderType!=
"") {
215 Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
216 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(fwdOp,
true);
218 Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm,blkFwdOp);
220 if(prec_Op==Teuchos::null) {
221 Teko::ModifiableLinearOp reorderedPrec =
Teko::buildInverse(*invFactory_,blockedFwdOp);
225 Teko::ModifiableLinearOp reorderedPrec = Teuchos::rcp_dynamic_cast<ReorderedLinearOp>(prec_Op,
true)->getBlockedOp();
231 if(prec_Op==Teuchos::null)
241 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<
" sec\n";
245 defaultPrec.initializeUnspecified( prec_Op);
246 defaultPrec.incrIter();
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";
252 *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
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
262 using Teuchos::outArg;
265 using Teuchos::rcp_dynamic_cast;
266 using Teuchos::implicit_cast;
268 Teuchos::Time totalTimer(
""), timer(
"");
269 totalTimer.start(
true);
271 const RCP<Teuchos::FancyOStream> out = this->getOStream();
272 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
274 bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
276 Teuchos::OSTab tab(out);
278 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
280 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
282 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
283 TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
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)
299 StratimikosFactoryPreconditioner & defaultPrec
300 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
303 RCP<Thyra::EpetraLinearOp> epetra_precOp
304 = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),
true);
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);
312 const bool startingOver = (teko_precOp == Teuchos::null);
315 invLib_ = Teuchos::null;
316 invFactory_ = Teuchos::null;
320 *out <<
"\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
324 std::stringstream ss;
325 ss << paramList_->get<std::string>(
"Strided Blocking");
328 while(not ss.eof()) {
332 TEUCHOS_ASSERT(num>0);
334 decomp_.push_back(num);
338 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist(
"Inverse Factory Library"),builder_);
339 invLib_->setRequestHandler(reqHandler_);
342 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
349 Teuchos::OSTab(out).o() <<
"> Creation time = "<<timer.totalElapsedTime()<<
" sec\n";
353 *out <<
"\nComputing the preconditioner ...\n";
357 bool writeBlockOps = paramList_->get<
bool>(
"Write Block Operator");
359 teko_precOp->initInverse();
361 if(decomp_.size()==1) {
362 teko_precOp->rebuildInverseOperator(epetraFwdOp);
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));
372 Teuchos::RCP<Epetra_Operator> wrappedFwdOp
373 = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
377 std::stringstream ss;
378 ss <<
"block-" << defaultPrec.getIter();
381 Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac
382 = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedFwdOp);
383 if(stridedJac!=Teuchos::null) {
385 stridedJac->WriteBlocks(ss.str());
387 else TEUCHOS_ASSERT(
false);
390 teko_precOp->rebuildInverseOperator(wrappedFwdOp);
397 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<
" sec\n";
401 epetra_precOp = rcp(
new Thyra::EpetraLinearOp);
403 epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
404 ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
405 ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
408 defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
409 defaultPrec.incrIter();
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";
415 *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
419void StratimikosFactory::uninitializePrec(
420 Thyra::PreconditionerBase<double> * ,
421 Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > * ,
422 Thyra::ESupportSolveUse *
425 TEUCHOS_TEST_FOR_EXCEPT(
true);
432void StratimikosFactory::setParameterList(
433 Teuchos::RCP<Teuchos::ParameterList>
const& paramList
436 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
438 paramList->validateParametersAndSetDefaults(*this->getValidParameters(),0);
439 paramList_ = paramList;
444Teuchos::RCP<Teuchos::ParameterList>
445StratimikosFactory::getNonconstParameterList()
451Teuchos::RCP<Teuchos::ParameterList>
452StratimikosFactory::unsetParameterList()
454 Teuchos::RCP<ParameterList> _paramList = paramList_;
455 paramList_ = Teuchos::null;
460Teuchos::RCP<const Teuchos::ParameterList>
461StratimikosFactory::getParameterList()
const
467Teuchos::RCP<const Teuchos::ParameterList>
468StratimikosFactory::getValidParameters()
const
472 using Teuchos::tuple;
473 using Teuchos::implicit_cast;
474 using Teuchos::rcp_implicit_cast;
476 static RCP<const ParameterList> validPL;
478 if(is_null(validPL)) {
480 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
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"
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"
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.");
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,
526 Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
555 if(wrappedOp==Teuchos::null)
558 std::vector<std::vector<int> > vars;
559 buildStridedVectors(*Jac,decomp_,vars);
563 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
564 if(reorderType!=
"") {
568 Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->Reorder(*brm);
572 Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->RebuildOps();
576 if(paramList_->get<
bool>(
"Test Block Operator")) {
578 = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
580 out <<
"Teko: Tested operator correctness: " << (result ?
"passed" :
"FAILED!") << std::endl;
587std::string StratimikosFactory::description()
const
589 std::ostringstream oss;
590 oss <<
"Teko::StratimikosFactory";
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
599 const Epetra_Map & rangeMap = Jac.OperatorRangeMap();
603 for(std::size_t i=0;i<decomp.size();i++)
604 numVars += decomp[i];
607 TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars)==0);
608 TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars)==0);
610 int * globalIds = rangeMap.MyGlobalElements();
612 vars.resize(decomp.size());
613 for(
int i=0;i<rangeMap.NumMyElements();) {
616 for(std::size_t d=0;d<decomp.size();d++) {
618 int current = decomp[d];
619 for(
int v=0;v<current;v++,i++)
620 vars[d].push_back(globalIds[i]);
626void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
627 const std::string & stratName)
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!");
632 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
633 = Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
636 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(
new TekoFactoryBuilder(builderCopy,Teuchos::null));
637 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
640void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
641 const Teuchos::RCP<Teko::RequestHandler> & rh,
642 const std::string & stratName)
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!");
647 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
648 = Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
652 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(
new TekoFactoryBuilder(builderCopy,rh));
653 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
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.