Teko Version of the Day
Teko_StratimikosFactory.cpp
00001 #include "Teko_StratimikosFactory.hpp"
00002 
00003 #include "Teuchos_Time.hpp"
00004 #include "Teuchos_AbstractFactoryStd.hpp"
00005 
00006 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
00007 #include "Thyra_EpetraLinearOp.hpp"
00008 #include "Thyra_DefaultPreconditioner.hpp"
00009 
00010 #include "Teko_InverseLibrary.hpp"
00011 #include "Teko_Preconditioner.hpp"
00012 #include "Teko_InverseFactoryOperator.hpp"
00013 #include "Teko_Utilities.hpp"
00014 #include "Teko_InverseLibrary.hpp"
00015 #include "Teko_StridedEpetraOperator.hpp"
00016 
00017 #include "EpetraExt_RowMatrixOut.h"
00018 
00019 namespace Teko {
00020 
00021 using Teuchos::RCP;
00022 using Teuchos::ParameterList;
00023 
00024 // hide stuff
00025 namespace {
00026    // Simple preconditioner class that adds a counter 
00027    class StratimikosFactoryPreconditioner :  public Thyra::DefaultPreconditioner<double>{
00028    public:
00029       StratimikosFactoryPreconditioner() : iter_(0) {}
00030 
00031       inline void incrIter() { iter_++; }
00032       inline std::size_t getIter() { return iter_; }
00033 
00034    private:
00035       StratimikosFactoryPreconditioner(const StratimikosFactoryPreconditioner &);
00036 
00037       std::size_t iter_;
00038    };
00039 
00040    // factory used to initialize the Teko::StratimikosFactory
00041    // user data
00042    class TekoFactoryBuilder 
00043          : public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
00044    public:
00045       TekoFactoryBuilder(const Teuchos::RCP<Teko::RequestHandler> & rh) : requestHandler_(rh) {}
00046       Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create() const
00047       { return Teuchos::rcp(new StratimikosFactory(requestHandler_)); }
00048  
00049    private:
00050       Teuchos::RCP<Teko::RequestHandler> requestHandler_;
00051    };
00052 }
00053 
00054 // Constructors/initializers/accessors
00055 StratimikosFactory::StratimikosFactory()
00056   :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
00057 {}
00058 
00059 // Constructors/initializers/accessors
00060 StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Teko::RequestHandler> & rh)
00061   :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
00062 {
00063    setRequestHandler(rh);
00064 }
00065 
00066 
00067 // Overridden from PreconditionerFactoryBase
00068 bool StratimikosFactory::isCompatible(const Thyra::LinearOpSourceBase<double> &fwdOpSrc) const
00069 {
00070    using Teuchos::outArg;
00071 
00072    Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00073    Thyra::EOpTransp epetraFwdOpTransp;
00074    Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
00075    Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
00076    double epetraFwdOpScalar;
00077   
00078    // check to make sure this is an epetra CrsMatrix
00079    Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc.getOp();
00080    epetraFwdOpViewExtractor_->getEpetraOpView(
00081            fwdOp,
00082            outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
00083            outArg(epetraFwdOpApplyAs),
00084            outArg(epetraFwdOpAdjointSupport),
00085            outArg(epetraFwdOpScalar));
00086 
00087    if( !dynamic_cast<const Epetra_CrsMatrix*>(&*epetraFwdOp) )
00088       return false;
00089 
00090    return true;
00091 }
00092 
00093 
00094 bool StratimikosFactory::applySupportsConj(Thyra::EConj conj) const
00095 {
00096   return false;
00097 }
00098 
00099 
00100 bool StratimikosFactory::applyTransposeSupportsConj(Thyra::EConj conj) const
00101 {
00102   return false; // See comment below
00103 }
00104 
00105 
00106 Teuchos::RCP<Thyra::PreconditionerBase<double> >
00107 StratimikosFactory::createPrec() const
00108 {
00109   return Teuchos::rcp(new StratimikosFactoryPreconditioner());
00110 }
00111 
00112 
00113 void StratimikosFactory::initializePrec(
00114   const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
00115   Thyra::PreconditionerBase<double> *prec,
00116   const Thyra::ESupportSolveUse supportSolveUse
00117   ) const
00118 {
00119   using Teuchos::outArg;
00120   using Teuchos::RCP;
00121   using Teuchos::rcp;
00122   using Teuchos::rcp_dynamic_cast;
00123   using Teuchos::implicit_cast;
00124 
00125   Teuchos::Time totalTimer(""), timer("");
00126   totalTimer.start(true);
00127 
00128   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00129   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00130 
00131   bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
00132 
00133   Teuchos::OSTab tab(out);
00134   if(mediumVerbosity)
00135     *out << "\nEntering Teko::StratimikosFactory::initializePrec(...) ...\n";
00136 
00137   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00138 #ifdef _DEBUG
00139   TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00140   TEST_FOR_EXCEPT(prec==NULL);
00141 #endif
00142 
00143   //
00144   // Unwrap and get the forward Epetra_Operator object
00145   //
00146   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00147   Thyra::EOpTransp epetraFwdOpTransp;
00148   Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
00149   Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
00150   double epetraFwdOpScalar;
00151   epetraFwdOpViewExtractor_->getEpetraOpView(
00152     fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
00153     outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
00154                                              );
00155   // Get the concrete precondtioner object
00156   StratimikosFactoryPreconditioner & defaultPrec 
00157         = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
00158 
00159   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00160   RCP<Thyra::EpetraLinearOp> epetra_precOp 
00161         = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),true);
00162 
00163   // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
00164   Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
00165   if(epetra_precOp!=Teuchos::null)
00166     teko_precOp = rcp_dynamic_cast<Teko::Epetra::InverseFactoryOperator>(epetra_precOp->epetra_op(),true);
00167 
00168   // Permform initialization if needed
00169   const bool startingOver = (teko_precOp == Teuchos::null);
00170   if(startingOver) {
00171 
00172     if(mediumVerbosity)
00173       *out << "\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
00174 
00175     timer.start(true);
00176 
00177     std::stringstream ss;
00178     ss << paramList_->get<std::string>("Strided Blocking");
00179 
00180     // figure out the decomposition requested by the string
00181     while(not ss.eof()) {
00182        int num=0;
00183        ss >> num;
00184 
00185        TEUCHOS_ASSERT(num>0);
00186 
00187        decomp_.push_back(num);
00188     }
00189 
00190     // build library, and set request handler (user defined!)
00191     invLib_  = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"));
00192     invLib_->setRequestHandler(reqHandler_);
00193 
00194     // build preconditioner factory
00195     invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
00196 
00197     // Create the initial preconditioner: DO NOT compute it yet
00198     teko_precOp = rcp( new Teko::Epetra::InverseFactoryOperator(invFactory_));
00199     
00200     timer.stop();
00201     if(mediumVerbosity)
00202       Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00203   }
00204 
00205   if(mediumVerbosity)
00206     *out << "\nComputing the preconditioner ...\n";
00207 
00208   // construct preconditioner
00209   {
00210      bool writeBlockOps = paramList_->get<bool>("Write Block Operator");
00211      timer.start(true);
00212      teko_precOp->initInverse();
00213 
00214      if(decomp_.size()==1) {
00215         teko_precOp->rebuildInverseOperator(epetraFwdOp);
00216 
00217         // write out to disk
00218         if(writeBlockOps) {
00219            std::stringstream ss;
00220            ss << "block-" << defaultPrec.getIter() << "_00.mm";
00221            EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*rcp_dynamic_cast<const Epetra_CrsMatrix>(epetraFwdOp));
00222         }
00223      }
00224      else {
00225         Teuchos::RCP<Epetra_Operator> wrappedFwdOp 
00226            = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
00227 
00228         // write out to disk
00229         if(writeBlockOps) {
00230            std::stringstream ss;
00231            ss << "block-" << defaultPrec.getIter();
00232            Teuchos::RCP<Teko::Epetra::StridedEpetraOperator> stridedJac
00233                  = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedFwdOp);
00234            if(stridedJac!=Teuchos::null) {
00235               // write out blocks of strided operator
00236               stridedJac->WriteBlocks(ss.str());
00237            }
00238            else TEUCHOS_ASSERT(false);
00239         }
00240    
00241         teko_precOp->rebuildInverseOperator(wrappedFwdOp);
00242      }
00243 
00244      timer.stop();
00245   }
00246 
00247   if(mediumVerbosity)
00248     Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
00249 
00250   // Initialize the output EpetraLinearOp
00251   if(startingOver) {
00252     epetra_precOp = rcp(new Thyra::EpetraLinearOp);
00253   }
00254   epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
00255     ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
00256     ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
00257 
00258   // Initialize the preconditioner
00259   defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
00260   defaultPrec.incrIter();
00261   totalTimer.stop();
00262 
00263   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00264     *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00265   if(mediumVerbosity)
00266     *out << "\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
00267 }
00268 
00269 
00270 void StratimikosFactory::uninitializePrec(
00271   Thyra::PreconditionerBase<double> *prec,
00272   Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > *fwdOp,
00273   Thyra::ESupportSolveUse *supportSolveUse
00274   ) const
00275 {
00276   TEST_FOR_EXCEPT(true);
00277 }
00278 
00279 
00280 // Overridden from ParameterListAcceptor
00281 
00282 
00283 void StratimikosFactory::setParameterList(
00284   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00285   )
00286 {
00287    TEST_FOR_EXCEPT(paramList.get()==NULL);
00288 
00289    paramList->validateParametersAndSetDefaults(*this->getValidParameters(),0);
00290    paramList_ = paramList;
00291 
00292 }
00293 
00294 
00295 Teuchos::RCP<Teuchos::ParameterList>
00296 StratimikosFactory::getNonconstParameterList()
00297 {
00298   return paramList_;
00299 }
00300 
00301 
00302 Teuchos::RCP<Teuchos::ParameterList>
00303 StratimikosFactory::unsetParameterList()
00304 {
00305   Teuchos::RCP<ParameterList> _paramList = paramList_;
00306   paramList_ = Teuchos::null;
00307   return _paramList;
00308 }
00309 
00310 
00311 Teuchos::RCP<const Teuchos::ParameterList>
00312 StratimikosFactory::getParameterList() const
00313 {
00314   return paramList_;
00315 }
00316 
00317 
00318 Teuchos::RCP<const Teuchos::ParameterList>
00319 StratimikosFactory::getValidParameters() const
00320 {
00321 
00322   using Teuchos::rcp;
00323   using Teuchos::tuple;
00324   using Teuchos::implicit_cast;
00325   using Teuchos::rcp_implicit_cast;
00326   typedef Teuchos::ParameterEntryValidator PEV;
00327 
00328   static RCP<const ParameterList> validPL;
00329 
00330   if(is_null(validPL)) {
00331 
00332     Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00333 
00334     pl->set("Test Block Operator",false,
00335             "If Stratiikos/Teko is used to break an operator into its parts,\n"
00336             "then setting this parameter to true will compare applications of the\n"
00337             "segregated operator to the original operator.");
00338     pl->set("Write Block Operator",false,
00339             "Write out the segregated operator to disk with the name \"block-?_xx\"");
00340     pl->set("Strided Blocking","1",
00341             "Assuming that the user wants Strided blocking, break the operating into\n"
00342             "blocks. The syntax can be thought to be associated with the solution\n"
00343             "vector. For example if your variables are [u v w p T], and we want [u v w]\n"
00344             "blocked together, and p and T seperate then the relevant string is \"3 1 1\".\n"
00345             "Meaning put the first 3 unknowns per node together and sperate the v and w\n"
00346             "components.");
00347     pl->set("Reorder Type","",
00348             "This specifies how the blocks are reordered for use in the precondtioner.\n"
00349             "For example, assume the linear system is generated from 3D Navier-Stokes\n"
00350             "with an energy equation, yielding the unknowns [u v w p T]. If the\n"
00351             "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n"
00352             "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n"
00353             "velocity and pressure forming an inner two-by-two block, and then the\n"
00354             "temperature unknowns forming a two-by-two system with the velocity-pressure\n"
00355             "block.");
00356     pl->set("Inverse Type","Amesos", 
00357             "The type of inverse operator the user wants. This can be one of the defaults\n"
00358             "from Stratimikos, or a Teko preconditioner defined in the\n"
00359             "\"Inverse Factory Library\".");
00360     pl->sublist("Inverse Factory Library",false,"Definition of Teko preconditioners."); 
00361 
00362     validPL = pl;
00363   }
00364 
00365   return validPL;
00366 }
00367 
00371 Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
00372                                                                  const Teuchos::RCP<const Epetra_Operator> & Jac,
00373                                                                  const Teuchos::RCP<Epetra_Operator> & wrapInput,
00374                                                                  std::ostream & out
00375                                                                  ) const
00376 {
00377    Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
00378 
00379    // initialize jacobian
00380    if(wrappedOp==Teuchos::null)
00381    {
00382       wrappedOp = Teuchos::rcp(new Teko::Epetra::StridedEpetraOperator(decomp_,Jac));
00383 
00384       // reorder the blocks if requested
00385       std::string reorderType = paramList_->get<std::string>("Reorder Type");
00386       if(reorderType!="") {
00387          RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
00388 
00389          // out << "Teko: Reordering = " << brm->toString() << std::endl;
00390          Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->Reorder(*brm);
00391       }
00392    }
00393    else {
00394       Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->RebuildOps();
00395    }
00396 
00397    // test blocked operator for correctness
00398    if(paramList_->get<bool>("Test Block Operator")) {
00399       bool result
00400          = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
00401 
00402       out << "Teko: Tested operator correctness:  " << (result ? "passed" : "FAILED!") << std::endl;
00403    }
00404 
00405    return wrappedOp;
00406 }
00407 
00408 std::string StratimikosFactory::description() const
00409 {
00410   std::ostringstream oss;
00411   oss << "Teko::StratimikosFactory";
00412   return oss.str();
00413 }
00414 
00415 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
00416                                const std::string & stratName)
00417 {
00418    TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
00419                       "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
00420 
00421    // use default constructor to add Teko::StratimikosFactory
00422    builder.setPreconditioningStrategyFactory(
00423          Teuchos::abstractFactoryStd<Thyra::PreconditionerFactoryBase<double>,Teko::StratimikosFactory>(),
00424          stratName);
00425 }
00426 
00427 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
00428                                const Teuchos::RCP<Teko::RequestHandler> & rh,
00429                                const std::string & stratName)
00430 {
00431    TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
00432                       "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
00433 
00434    // build an instance of a Teuchos::AbsractFactory<Thyra::PFB> so request handler is passed onto
00435    // the resulting StratimikosFactory
00436    Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(rh));
00437    builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
00438 }
00439 
00440 } // namespace Thyra
 All Classes Files Functions Variables