1 #include "Teko_StratimikosFactory.hpp" 3 #include "Teuchos_Time.hpp" 4 #include "Teuchos_AbstractFactoryStd.hpp" 6 #include "Thyra_EpetraOperatorViewExtractorStd.hpp" 7 #include "Thyra_EpetraLinearOp.hpp" 8 #include "Thyra_DefaultPreconditioner.hpp" 10 #include "Teko_InverseLibrary.hpp" 11 #include "Teko_Preconditioner.hpp" 12 #include "Teko_InverseFactoryOperator.hpp" 14 #include "Teko_InverseLibrary.hpp" 15 #include "Teko_StridedEpetraOperator.hpp" 16 #include "Teko_BlockedEpetraOperator.hpp" 17 #include "Teko_ReorderedLinearOp.hpp" 19 #include "EpetraExt_RowMatrixOut.h" 24 using Teuchos::ParameterList;
29 class StratimikosFactoryPreconditioner :
public Thyra::DefaultPreconditioner<double>{
31 StratimikosFactoryPreconditioner() : iter_(0) {}
33 inline void incrIter() { iter_++; }
34 inline std::size_t getIter() {
return iter_; }
37 StratimikosFactoryPreconditioner(
const StratimikosFactoryPreconditioner &);
44 class TekoFactoryBuilder
45 :
public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
47 TekoFactoryBuilder(
const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> & builder,
48 const Teuchos::RCP<Teko::RequestHandler> & rh) : builder_(builder), requestHandler_(rh) {}
49 Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create()
const 50 {
return Teuchos::rcp(
new StratimikosFactory(builder_,requestHandler_)); }
53 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builder_;
54 Teuchos::RCP<Teko::RequestHandler> requestHandler_;
60 :epetraFwdOpViewExtractor_(
Teuchos::rcp(new
Thyra::EpetraOperatorViewExtractorStd()))
65 :epetraFwdOpViewExtractor_(
Teuchos::rcp(new
Thyra::EpetraOperatorViewExtractorStd()))
71 const Teuchos::RCP<Teko::RequestHandler> & rh)
72 :epetraFwdOpViewExtractor_(
Teuchos::rcp(new
Thyra::EpetraOperatorViewExtractorStd())), builder_(builder)
81 using Teuchos::outArg;
83 TEUCHOS_ASSERT(
false);
85 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
86 Thyra::EOpTransp epetraFwdOpTransp;
87 Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
88 Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
89 double epetraFwdOpScalar;
92 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc.getOp();
93 epetraFwdOpViewExtractor_->getEpetraOpView(
95 outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
96 outArg(epetraFwdOpApplyAs),
97 outArg(epetraFwdOpAdjointSupport),
98 outArg(epetraFwdOpScalar));
100 if( !dynamic_cast<const Epetra_CrsMatrix*>(&*epetraFwdOp) )
119 Teuchos::RCP<Thyra::PreconditionerBase<double> >
122 return Teuchos::rcp(
new StratimikosFactoryPreconditioner());
127 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
128 Thyra::PreconditionerBase<double> *prec,
129 const Thyra::ESupportSolveUse supportSolveUse
132 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
134 if(epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
141 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
142 Thyra::PreconditionerBase<double> *prec,
143 const Thyra::ESupportSolveUse supportSolveUse
148 using Teuchos::rcp_dynamic_cast;
149 using Teuchos::implicit_cast;
150 using Thyra::LinearOpBase;
152 Teuchos::Time totalTimer(
""), timer(
"");
153 totalTimer.start(
true);
155 const RCP<Teuchos::FancyOStream> out = this->getOStream();
156 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
158 bool mediumVerbosity = (out.get() && implicit_cast<
int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
160 Teuchos::OSTab tab(out);
162 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
164 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
167 StratimikosFactoryPreconditioner & defaultPrec
168 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
169 Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
172 const bool startingOver = (prec_Op == Teuchos::null);
175 invLib_ = Teuchos::null;
176 invFactory_ = Teuchos::null;
179 *out <<
"\nCreating the initial Teko Operator object...\n";
184 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist(
"Inverse Factory Library"),builder_);
185 invLib_->setRequestHandler(reqHandler_);
188 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
192 Teuchos::OSTab(out).o() <<
"> Creation time = "<<timer.totalElapsedTime()<<
" sec\n";
196 *out <<
"\nComputing the preconditioner ...\n";
199 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
200 if(reorderType!=
"") {
202 Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
203 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(fwdOp,
true);
205 Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm,blkFwdOp);
207 if(prec_Op==Teuchos::null) {
208 Teko::ModifiableLinearOp reorderedPrec =
Teko::buildInverse(*invFactory_,blockedFwdOp);
212 Teko::ModifiableLinearOp reorderedPrec = Teuchos::rcp_dynamic_cast<
ReorderedLinearOp>(prec_Op,
true)->getBlockedOp();
218 if(prec_Op==Teuchos::null)
228 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<
" sec\n";
232 defaultPrec.initializeUnspecified( prec_Op);
233 defaultPrec.incrIter();
236 if(out.get() && implicit_cast<
int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
237 *out <<
"\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<
" sec\n";
239 *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
243 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
244 Thyra::PreconditionerBase<double> *prec,
245 const Thyra::ESupportSolveUse supportSolveUse
248 using Teuchos::outArg;
251 using Teuchos::rcp_dynamic_cast;
252 using Teuchos::implicit_cast;
254 Teuchos::Time totalTimer(
""), timer(
"");
255 totalTimer.start(
true);
257 const RCP<Teuchos::FancyOStream> out = this->getOStream();
258 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
260 bool mediumVerbosity = (out.get() && implicit_cast<
int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
262 Teuchos::OSTab tab(out);
264 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
266 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
268 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
269 TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
275 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
276 Thyra::EOpTransp epetraFwdOpTransp;
277 Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
278 Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
279 double epetraFwdOpScalar;
280 epetraFwdOpViewExtractor_->getEpetraOpView(
281 fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
282 outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
285 StratimikosFactoryPreconditioner & defaultPrec
286 = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
289 RCP<Thyra::EpetraLinearOp> epetra_precOp
290 = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),
true);
293 Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
294 if(epetra_precOp!=Teuchos::null)
298 const bool startingOver = (teko_precOp == Teuchos::null);
301 invLib_ = Teuchos::null;
302 invFactory_ = Teuchos::null;
306 *out <<
"\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
310 std::stringstream ss;
311 ss << paramList_->get<std::string>(
"Strided Blocking");
314 while(not ss.eof()) {
318 TEUCHOS_ASSERT(num>0);
320 decomp_.push_back(num);
324 invLib_ = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist(
"Inverse Factory Library"),builder_);
325 invLib_->setRequestHandler(reqHandler_);
328 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
335 Teuchos::OSTab(out).o() <<
"> Creation time = "<<timer.totalElapsedTime()<<
" sec\n";
339 *out <<
"\nComputing the preconditioner ...\n";
343 bool writeBlockOps = paramList_->get<
bool>(
"Write Block Operator");
345 teko_precOp->initInverse();
347 if(decomp_.size()==1) {
348 teko_precOp->rebuildInverseOperator(epetraFwdOp);
352 std::stringstream ss;
353 ss <<
"block-" << defaultPrec.getIter() <<
"_00.mm";
354 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*rcp_dynamic_cast<
const Epetra_CrsMatrix>(epetraFwdOp));
358 Teuchos::RCP<Epetra_Operator> wrappedFwdOp
359 = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
363 std::stringstream ss;
364 ss <<
"block-" << defaultPrec.getIter();
367 Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac
369 if(stridedJac!=Teuchos::null) {
373 else TEUCHOS_ASSERT(
false);
376 teko_precOp->rebuildInverseOperator(wrappedFwdOp);
383 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<
" sec\n";
387 epetra_precOp = rcp(
new Thyra::EpetraLinearOp);
389 epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
390 ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
391 ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
394 defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
395 defaultPrec.incrIter();
398 if(out.get() && implicit_cast<
int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
399 *out <<
"\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<
" sec\n";
401 *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
406 Thyra::PreconditionerBase<double> *prec,
407 Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > *fwdOp,
408 Thyra::ESupportSolveUse *supportSolveUse
411 TEUCHOS_TEST_FOR_EXCEPT(
true);
419 Teuchos::RCP<Teuchos::ParameterList>
const& paramList
422 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
425 paramList_ = paramList;
430 Teuchos::RCP<Teuchos::ParameterList>
437 Teuchos::RCP<Teuchos::ParameterList>
440 Teuchos::RCP<ParameterList> _paramList = paramList_;
441 paramList_ = Teuchos::null;
446 Teuchos::RCP<const Teuchos::ParameterList>
453 Teuchos::RCP<const Teuchos::ParameterList>
458 using Teuchos::tuple;
459 using Teuchos::implicit_cast;
460 using Teuchos::rcp_implicit_cast;
462 static RCP<const ParameterList> validPL;
464 if(is_null(validPL)) {
466 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
468 pl->set(
"Test Block Operator",
false,
469 "If Stratiikos/Teko is used to break an operator into its parts,\n" 470 "then setting this parameter to true will compare applications of the\n" 471 "segregated operator to the original operator.");
472 pl->set(
"Write Block Operator",
false,
473 "Write out the segregated operator to disk with the name \"block-?_xx\"");
474 pl->set(
"Strided Blocking",
"1",
475 "Assuming that the user wants Strided blocking, break the operator into\n" 476 "blocks. The syntax can be thought to be associated with the solution\n" 477 "vector. For example if your variables are [u v w p T], and we want [u v w]\n" 478 "blocked together, and p and T separate then the relevant string is \"3 1 1\".\n" 479 "Meaning put the first 3 unknowns per node together and separate the v and w\n" 481 pl->set(
"Reorder Type",
"",
482 "This specifies how the blocks are reordered for use in the preconditioner.\n" 483 "For example, assume the linear system is generated from 3D Navier-Stokes\n" 484 "with an energy equation, yielding the unknowns [u v w p T]. If the\n" 485 "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n" 486 "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n" 487 "velocity and pressure forming an inner two-by-two block, and then the\n" 488 "temperature unknowns forming a two-by-two system with the velocity-pressure\n" 490 pl->set(
"Inverse Type",
"Amesos",
491 "The type of inverse operator the user wants. This can be one of the defaults\n" 492 "from Stratimikos, or a Teko preconditioner defined in the\n" 493 "\"Inverse Factory Library\".");
494 pl->sublist(
"Inverse Factory Library",
false,
"Definition of Teko preconditioners.");
505 Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
506 const Teuchos::RCP<const Epetra_Operator> & Jac,
507 const Teuchos::RCP<Epetra_Operator> & wrapInput,
511 Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
540 if(wrappedOp==Teuchos::null)
543 std::vector<std::vector<int> > vars;
544 buildStridedVectors(*Jac,decomp_,vars);
548 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
549 if(reorderType!=
"") {
561 if(paramList_->get<
bool>(
"Test Block Operator")) {
565 out <<
"Teko: Tested operator correctness: " << (result ?
"passed" :
"FAILED!") << std::endl;
573 std::ostringstream oss;
574 oss <<
"Teko::StratimikosFactory";
578 void StratimikosFactory::buildStridedVectors(
const Epetra_Operator & Jac,
579 const std::vector<int> & decomp,
580 std::vector<std::vector<int> > & vars)
const 582 const Epetra_Map & rangeMap = Jac.OperatorRangeMap();
586 for(std::size_t i=0;i<decomp.size();i++)
587 numVars += decomp[i];
590 TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars)==0);
591 TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars)==0);
593 int * globalIds = rangeMap.MyGlobalElements();
595 vars.resize(decomp.size());
596 for(
int i=0;i<rangeMap.NumMyElements();) {
599 for(std::size_t d=0;d<decomp.size();d++) {
601 int current = decomp[d];
602 for(
int v=0;v<current;v++,i++)
603 vars[d].push_back(globalIds[i]);
608 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
609 const std::string & stratName)
611 TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist(
"Preconditioner Types").isParameter(stratName),std::logic_error,
612 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +
"\" because it is already included in builder!");
614 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
615 = Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
618 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(
new TekoFactoryBuilder(builderCopy,Teuchos::null));
619 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
622 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
623 const Teuchos::RCP<Teko::RequestHandler> & rh,
624 const std::string & stratName)
626 TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist(
"Preconditioner Types").isParameter(stratName),std::logic_error,
627 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +
"\" because it is already included in builder!");
629 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy
630 = Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
634 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(
new TekoFactoryBuilder(builderCopy,rh));
635 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
void uninitializePrec(Thyra::PreconditionerBase< double > *prec, Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > *fwdOp, Thyra::ESupportSolveUse *supportSolveUse) const
void setRequestHandler(const Teuchos::RCP< Teko::RequestHandler > &rh)
bool applyTransposeSupportsConj(Thyra::EConj conj) const
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
bool applySupportsConj(Thyra::EConj conj) const
virtual void WriteBlocks(const std::string &prefix) const
bool isCompatible(const Thyra::LinearOpSourceBase< double > &fwdOp) const
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
std::string description() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Thyra::PreconditionerBase< double > > createPrec() const
A single Epetra wrapper for all operators constructed from an inverse operator.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void initializePrec_Thyra(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
This class takes a blocked linear op and represents it in a flattened form.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
void initializePrec(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
void initializePrec_Epetra(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
Tear about a user specified Epetra_Operator (CrsMatrix) using a vector of vectors of GIDs for each bl...
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.