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"
25 using 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)
51 : builder_(builder), requestHandler_(rh) {}
52 Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create()
const {
53 return Teuchos::rcp(
new StratimikosFactory(builder_, requestHandler_));
57 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builder_;
58 Teuchos::RCP<Teko::RequestHandler> requestHandler_;
62 #ifdef TEKO_HAVE_EPETRA
66 : epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())) {}
70 : epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())) {
71 setRequestHandler(rh);
74 #endif // TEKO_HAVE_EPETRA
77 const Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> &builder,
78 const Teuchos::RCP<Teko::RequestHandler> &rh)
80 #ifdef TEKO_HAVE_EPETRA
81 epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd())),
84 setRequestHandler(rh);
89 const Thyra::LinearOpSourceBase<double> & )
const {
90 using Teuchos::outArg;
92 TEUCHOS_ASSERT(
false);
124 return Teuchos::rcp(
new StratimikosFactoryPreconditioner());
128 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
129 Thyra::PreconditionerBase<double> *prec,
const Thyra::ESupportSolveUse supportSolveUse)
const {
130 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
132 #ifdef TEKO_HAVE_EPETRA
133 if (epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
134 initializePrec_Epetra(fwdOpSrc, prec, supportSolveUse);
141 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
142 Thyra::PreconditionerBase<double> *prec,
const Thyra::ESupportSolveUse
144 using Teuchos::implicit_cast;
147 using Teuchos::rcp_dynamic_cast;
148 using Thyra::LinearOpBase;
150 Teuchos::Time totalTimer(
""), timer(
"");
151 totalTimer.start(
true);
153 const RCP<Teuchos::FancyOStream> out = this->getOStream();
154 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
156 bool mediumVerbosity =
157 (out.get() && implicit_cast<
int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
159 Teuchos::OSTab tab(out);
161 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
163 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
166 StratimikosFactoryPreconditioner &defaultPrec =
167 Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
168 Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
171 const bool startingOver = (prec_Op == Teuchos::null);
174 invLib_ = Teuchos::null;
175 invFactory_ = Teuchos::null;
177 if (mediumVerbosity) *out <<
"\nCreating the initial Teko Operator object...\n";
182 invLib_ = Teko::InverseLibrary::buildFromParameterList(
183 paramList_->sublist(
"Inverse Factory Library"), builder_);
184 invLib_->setRequestHandler(reqHandler_);
187 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
191 Teuchos::OSTab(out).o() <<
"> Creation time = " << timer.totalElapsedTime() <<
" sec\n";
194 if (mediumVerbosity) *out <<
"\nComputing the preconditioner ...\n";
197 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
198 if (reorderType !=
"") {
199 Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
200 Teuchos::rcp_dynamic_cast<
const Thyra::BlockedLinearOpBase<double> >(fwdOp,
true);
202 Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm, blkFwdOp);
204 if (prec_Op == Teuchos::null) {
205 Teko::ModifiableLinearOp reorderedPrec =
Teko::buildInverse(*invFactory_, blockedFwdOp);
208 Teko::ModifiableLinearOp reorderedPrec =
214 if (prec_Op == Teuchos::null)
224 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = " << timer.totalElapsedTime()
228 defaultPrec.initializeUnspecified(prec_Op);
229 defaultPrec.incrIter();
232 if (out.get() && implicit_cast<
int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
233 *out <<
"\nTotal time in Teko::StratimikosFactory = " << totalTimer.totalElapsedTime()
236 *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
239 #ifdef TEKO_HAVE_EPETRA
240 void StratimikosFactory::initializePrec_Epetra(
241 const Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
242 Thyra::PreconditionerBase<double> *prec,
const Thyra::ESupportSolveUse
244 using Teuchos::implicit_cast;
245 using Teuchos::outArg;
248 using Teuchos::rcp_dynamic_cast;
250 Teuchos::Time totalTimer(
""), timer(
"");
251 totalTimer.start(
true);
253 const RCP<Teuchos::FancyOStream> out = this->getOStream();
254 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
256 bool mediumVerbosity =
257 (out.get() && implicit_cast<
int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
259 Teuchos::OSTab tab(out);
261 *out <<
"\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
263 Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
265 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get() == NULL);
266 TEUCHOS_TEST_FOR_EXCEPT(prec == NULL);
272 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
273 Thyra::EOpTransp epetraFwdOpTransp;
274 Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
275 Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
276 double epetraFwdOpScalar;
277 epetraFwdOpViewExtractor_->getEpetraOpView(
278 fwdOp, outArg(epetraFwdOp), outArg(epetraFwdOpTransp), outArg(epetraFwdOpApplyAs),
279 outArg(epetraFwdOpAdjointSupport), outArg(epetraFwdOpScalar));
281 StratimikosFactoryPreconditioner &defaultPrec =
282 Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
285 RCP<Thyra::EpetraLinearOp> epetra_precOp =
286 rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),
true);
289 Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
290 if (epetra_precOp != Teuchos::null)
295 const bool startingOver = (teko_precOp == Teuchos::null);
298 invLib_ = Teuchos::null;
299 invFactory_ = Teuchos::null;
303 *out <<
"\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
307 std::stringstream ss;
308 ss << paramList_->get<std::string>(
"Strided Blocking");
311 while (not ss.eof()) {
315 TEUCHOS_ASSERT(num > 0);
317 decomp_.push_back(num);
321 invLib_ = Teko::InverseLibrary::buildFromParameterList(
322 paramList_->sublist(
"Inverse Factory Library"), builder_);
323 invLib_->setRequestHandler(reqHandler_);
326 invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>(
"Inverse Type"));
333 Teuchos::OSTab(out).o() <<
"> Creation time = " << timer.totalElapsedTime() <<
" sec\n";
336 if (mediumVerbosity) *out <<
"\nComputing the preconditioner ...\n";
340 bool writeBlockOps = paramList_->get<
bool>(
"Write Block Operator");
342 teko_precOp->initInverse();
344 if (decomp_.size() == 1) {
345 teko_precOp->rebuildInverseOperator(epetraFwdOp);
349 std::stringstream ss;
350 ss <<
"block-" << defaultPrec.getIter() <<
"_00.mm";
351 EpetraExt::RowMatrixToMatrixMarketFile(
352 ss.str().c_str(), *rcp_dynamic_cast<
const Epetra_CrsMatrix>(epetraFwdOp));
355 Teuchos::RCP<Epetra_Operator> wrappedFwdOp =
356 buildWrappedEpetraOperator(epetraFwdOp, teko_precOp->getNonconstForwardOp(), *out);
360 std::stringstream ss;
361 ss <<
"block-" << defaultPrec.getIter();
364 Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac =
366 if (stridedJac != Teuchos::null) {
370 TEUCHOS_ASSERT(
false);
373 teko_precOp->rebuildInverseOperator(wrappedFwdOp);
380 Teuchos::OSTab(out).o() <<
"=> Preconditioner construction time = " << timer.totalElapsedTime()
385 epetra_precOp = rcp(
new Thyra::EpetraLinearOp);
387 epetra_precOp->initialize(teko_precOp, epetraFwdOpTransp, Thyra::EPETRA_OP_APPLY_APPLY_INVERSE,
388 Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
391 defaultPrec.initializeUnspecified(
392 Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
393 defaultPrec.incrIter();
396 if (out.get() && implicit_cast<
int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
397 *out <<
"\nTotal time in Teko::StratimikosFactory = " << totalTimer.totalElapsedTime()
399 if (mediumVerbosity) *out <<
"\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
401 #endif // TEKO_HAVE_EPETRA
404 Thyra::PreconditionerBase<double> * ,
405 Teuchos::RCP<
const Thyra::LinearOpSourceBase<double> > * ,
406 Thyra::ESupportSolveUse *
408 TEUCHOS_TEST_FOR_EXCEPT(
true);
414 TEUCHOS_TEST_FOR_EXCEPT(paramList.get() == NULL);
417 paramList_ = paramList;
425 Teuchos::RCP<ParameterList> _paramList = paramList_;
426 paramList_ = Teuchos::null;
435 using Teuchos::implicit_cast;
437 using Teuchos::rcp_implicit_cast;
438 using Teuchos::tuple;
440 static RCP<const ParameterList> validPL;
442 if (is_null(validPL)) {
443 Teuchos::RCP<Teuchos::ParameterList> pl = rcp(
new Teuchos::ParameterList());
445 pl->set(
"Test Block Operator",
false,
446 "If Stratiikos/Teko is used to break an operator into its parts,\n"
447 "then setting this parameter to true will compare applications of the\n"
448 "segregated operator to the original operator.");
449 pl->set(
"Write Block Operator",
false,
450 "Write out the segregated operator to disk with the name \"block-?_xx\"");
451 pl->set(
"Strided Blocking",
"1",
452 "Assuming that the user wants Strided blocking, break the operator into\n"
453 "blocks. The syntax can be thought to be associated with the solution\n"
454 "vector. For example if your variables are [u v w p T], and we want [u v w]\n"
455 "blocked together, and p and T separate then the relevant string is \"3 1 1\".\n"
456 "Meaning put the first 3 unknowns per node together and separate the v and w\n"
458 pl->set(
"Reorder Type",
"",
459 "This specifies how the blocks are reordered for use in the preconditioner.\n"
460 "For example, assume the linear system is generated from 3D Navier-Stokes\n"
461 "with an energy equation, yielding the unknowns [u v w p T]. If the\n"
462 "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n"
463 "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n"
464 "velocity and pressure forming an inner two-by-two block, and then the\n"
465 "temperature unknowns forming a two-by-two system with the velocity-pressure\n"
467 std::string defaultInverseType =
"";
468 #if defined(Teko_ENABLE_Amesos)
469 defaultInverseType =
"Amesos";
470 #elif defined(Teko_ENABLE_Amesos2)
471 defaultInverseType =
"Amesos2";
473 pl->set(
"Inverse Type", defaultInverseType,
474 "The type of inverse operator the user wants. This can be one of the defaults\n"
475 "from Stratimikos, or a Teko preconditioner defined in the\n"
476 "\"Inverse Factory Library\".");
477 pl->sublist(
"Inverse Factory Library",
false,
"Definition of Teko preconditioners.");
485 #ifdef TEKO_HAVE_EPETRA
489 Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
490 const Teuchos::RCP<const Epetra_Operator> &Jac,
const Teuchos::RCP<Epetra_Operator> &wrapInput,
491 std::ostream &out)
const {
492 Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
524 if (wrappedOp == Teuchos::null) {
526 std::vector<std::vector<int> > vars;
527 buildStridedVectors(*Jac, decomp_, vars);
531 std::string reorderType = paramList_->get<std::string>(
"Reorder Type");
532 if (reorderType !=
"") {
543 if (paramList_->get<
bool>(
"Test Block Operator")) {
545 ->testAgainstFullOperator(600, 1e-14);
547 out <<
"Teko: Tested operator correctness: " << (result ?
"passed" :
"FAILED!") << std::endl;
552 #endif // TEKO_HAVE_EPETRA
555 std::ostringstream oss;
556 oss <<
"Teko::StratimikosFactory";
560 #ifdef TEKO_HAVE_EPETRA
561 void StratimikosFactory::buildStridedVectors(
const Epetra_Operator &Jac,
562 const std::vector<int> &decomp,
563 std::vector<std::vector<int> > &vars)
const {
564 const Epetra_Map &rangeMap = Jac.OperatorRangeMap();
568 for (std::size_t i = 0; i < decomp.size(); i++) numVars += decomp[i];
571 TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars) == 0);
572 TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars) == 0);
574 int *globalIds = rangeMap.MyGlobalElements();
576 vars.resize(decomp.size());
577 for (
int i = 0; i < rangeMap.NumMyElements();) {
579 for (std::size_t d = 0; d < decomp.size(); d++) {
581 int current = decomp[d];
582 for (
int v = 0; v < current; v++, i++) vars[d].push_back(globalIds[i]);
586 #endif // TEKO_HAVE_EPETRA
588 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder &builder,
589 const std::string &stratName) {
590 TEUCHOS_TEST_FOR_EXCEPTION(
591 builder.getValidParameters()->sublist(
"Preconditioner Types").isParameter(stratName),
593 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +
594 "\" because it is already included in builder!");
596 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy =
597 Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
600 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder =
601 Teuchos::rcp(
new TekoFactoryBuilder(builderCopy, Teuchos::null));
602 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder, stratName);
605 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder &builder,
606 const Teuchos::RCP<Teko::RequestHandler> &rh,
607 const std::string &stratName) {
608 TEUCHOS_TEST_FOR_EXCEPTION(
609 builder.getValidParameters()->sublist(
"Preconditioner Types").isParameter(stratName),
611 "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +
612 "\" because it is already included in builder!");
614 Teuchos::RCP<Stratimikos::DefaultLinearSolverBuilder> builderCopy =
615 Teuchos::rcp(
new Stratimikos::DefaultLinearSolverBuilder(builder));
619 Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder =
620 Teuchos::rcp(
new TekoFactoryBuilder(builderCopy, rh));
621 builder.setPreconditioningStrategyFactory(tekoFactoryBuilder, stratName);
bool applySupportsConj(Thyra::EConj conj) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void rebuildInverse(const InverseFactory &factory, const LinearOp &A, InverseLinearOp &invA)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void initializePrec(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
void uninitializePrec(Thyra::PreconditionerBase< double > *prec, Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > *fwdOp, Thyra::ESupportSolveUse *supportSolveUse) const
bool applyTransposeSupportsConj(Thyra::EConj conj) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
void initializePrec_Thyra(const Teuchos::RCP< const Thyra::LinearOpSourceBase< double > > &fwdOp, Thyra::PreconditionerBase< double > *prec, const Thyra::ESupportSolveUse supportSolveUse) const
A single Epetra wrapper for all operators constructed from an inverse operator.
bool isCompatible(const Thyra::LinearOpSourceBase< double > &fwdOp) const
Teuchos::RCP< Thyra::PreconditionerBase< double > > createPrec() const
InverseLinearOp buildInverse(const InverseFactory &factory, const LinearOp &A)
Build an inverse operator using a factory and a linear operator.
Teuchos::RCP< const BlockReorderManager > blockedReorderFromString(std::string &reorder)
Convert a string to a block reorder manager object.
This class takes a blocked linear op and represents it in a flattened form.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const ¶mList)
Tear about a user specified Epetra_Operator (CrsMatrix) using a vector of vectors of GIDs for each bl...
std::string description() const
virtual void WriteBlocks(const std::string &prefix) const