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 #include "Teko_BlockedEpetraOperator.hpp"
00017 #include "Teko_ReorderedLinearOp.hpp"
00018 
00019 #include "EpetraExt_RowMatrixOut.h"
00020 
00021 namespace Teko {
00022 
00023 using Teuchos::RCP;
00024 using Teuchos::ParameterList;
00025 
00026 // hide stuff
00027 namespace {
00028    // Simple preconditioner class that adds a counter 
00029    class StratimikosFactoryPreconditioner :  public Thyra::DefaultPreconditioner<double>{
00030    public:
00031       StratimikosFactoryPreconditioner() : iter_(0) {}
00032 
00033       inline void incrIter() { iter_++; }
00034       inline std::size_t getIter() { return iter_; }
00035 
00036    private:
00037       StratimikosFactoryPreconditioner(const StratimikosFactoryPreconditioner &);
00038 
00039       std::size_t iter_;
00040    };
00041 
00042    // factory used to initialize the Teko::StratimikosFactory
00043    // user data
00044    class TekoFactoryBuilder 
00045          : public Teuchos::AbstractFactory<Thyra::PreconditionerFactoryBase<double> > {
00046    public:
00047       TekoFactoryBuilder(const Teuchos::RCP<Teko::RequestHandler> & rh) : requestHandler_(rh) {}
00048       Teuchos::RCP<Thyra::PreconditionerFactoryBase<double> > create() const
00049       { return Teuchos::rcp(new StratimikosFactory(requestHandler_)); }
00050  
00051    private:
00052       Teuchos::RCP<Teko::RequestHandler> requestHandler_;
00053    };
00054 }
00055 
00056 // Constructors/initializers/accessors
00057 StratimikosFactory::StratimikosFactory()
00058   :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
00059 {}
00060 
00061 // Constructors/initializers/accessors
00062 StratimikosFactory::StratimikosFactory(const Teuchos::RCP<Teko::RequestHandler> & rh)
00063   :epetraFwdOpViewExtractor_(Teuchos::rcp(new Thyra::EpetraOperatorViewExtractorStd()))
00064 {
00065    setRequestHandler(rh);
00066 }
00067 
00068 
00069 // Overridden from PreconditionerFactoryBase
00070 bool StratimikosFactory::isCompatible(const Thyra::LinearOpSourceBase<double> &fwdOpSrc) const
00071 {
00072    using Teuchos::outArg;
00073 
00074    TEUCHOS_ASSERT(false); // what you doing?
00075 
00076    Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00077    Thyra::EOpTransp epetraFwdOpTransp;
00078    Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
00079    Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
00080    double epetraFwdOpScalar;
00081   
00082    // check to make sure this is an epetra CrsMatrix
00083    Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc.getOp();
00084    epetraFwdOpViewExtractor_->getEpetraOpView(
00085            fwdOp,
00086            outArg(epetraFwdOp),outArg(epetraFwdOpTransp),
00087            outArg(epetraFwdOpApplyAs),
00088            outArg(epetraFwdOpAdjointSupport),
00089            outArg(epetraFwdOpScalar));
00090 
00091    if( !dynamic_cast<const Epetra_CrsMatrix*>(&*epetraFwdOp) )
00092       return false;
00093 
00094    return true;
00095 }
00096 
00097 
00098 bool StratimikosFactory::applySupportsConj(Thyra::EConj conj) const
00099 {
00100   return false;
00101 }
00102 
00103 
00104 bool StratimikosFactory::applyTransposeSupportsConj(Thyra::EConj conj) const
00105 {
00106   return false; // See comment below
00107 }
00108 
00109 
00110 Teuchos::RCP<Thyra::PreconditionerBase<double> >
00111 StratimikosFactory::createPrec() const
00112 {
00113   return Teuchos::rcp(new StratimikosFactoryPreconditioner());
00114 }
00115 
00116 
00117 void StratimikosFactory::initializePrec(
00118   const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
00119   Thyra::PreconditionerBase<double> *prec,
00120   const Thyra::ESupportSolveUse supportSolveUse
00121   ) const
00122 {
00123   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00124 
00125   if(epetraFwdOpViewExtractor_->isCompatible(*fwdOp))
00126      initializePrec_Epetra(fwdOpSrc,prec,supportSolveUse);
00127   else
00128      initializePrec_Thyra(fwdOpSrc,prec,supportSolveUse);
00129 }
00130 
00131 void StratimikosFactory::initializePrec_Thyra(
00132   const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
00133   Thyra::PreconditionerBase<double> *prec,
00134   const Thyra::ESupportSolveUse supportSolveUse
00135   ) const
00136 {
00137   using Teuchos::RCP;
00138   using Teuchos::rcp;
00139   using Teuchos::rcp_dynamic_cast;
00140   using Teuchos::implicit_cast;
00141   using Thyra::LinearOpBase;
00142 
00143   Teuchos::Time totalTimer(""), timer("");
00144   totalTimer.start(true);
00145 
00146   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00147   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00148 
00149   bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
00150 
00151   Teuchos::OSTab tab(out);
00152   if(mediumVerbosity)
00153     *out << "\nEntering Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
00154 
00155   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00156 
00157   // Get the concrete preconditioner object
00158   StratimikosFactoryPreconditioner & defaultPrec 
00159         = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
00160   Teuchos::RCP<LinearOpBase<double> > prec_Op = defaultPrec.getNonconstUnspecifiedPrecOp();
00161 
00162   // Permform initialization if needed
00163   const bool startingOver = (prec_Op == Teuchos::null);
00164   if(startingOver) {
00165     // start over
00166     invLib_ = Teuchos::null;
00167     invFactory_ = Teuchos::null;
00168 
00169     if(mediumVerbosity)
00170       *out << "\nCreating the initial Teko Operator object...\n";
00171 
00172     timer.start(true);
00173 
00174     // build library, and set request handler (user defined!)
00175     invLib_  = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"));
00176     invLib_->setRequestHandler(reqHandler_);
00177 
00178     // build preconditioner factory
00179     invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
00180 
00181     timer.stop();
00182     if(mediumVerbosity)
00183       Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00184   }
00185 
00186   if(mediumVerbosity)
00187     *out << "\nComputing the preconditioner ...\n";
00188 
00189   // setup reordering if required
00190   std::string reorderType = paramList_->get<std::string>("Reorder Type");
00191   if(reorderType!="") {
00192      
00193      Teuchos::RCP<const Thyra::BlockedLinearOpBase<double> > blkFwdOp =
00194          Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(fwdOp,true);
00195      RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
00196      Teko::LinearOp blockedFwdOp = Teko::buildReorderedLinearOp(*brm,blkFwdOp);
00197 
00198      if(prec_Op==Teuchos::null) {
00199         Teko::ModifiableLinearOp reorderedPrec = Teko::buildInverse(*invFactory_,blockedFwdOp);
00200         prec_Op = Teuchos::rcp(new ReorderedLinearOp(brm,reorderedPrec));
00201      }
00202      else {
00203         Teko::ModifiableLinearOp reorderedPrec = Teuchos::rcp_dynamic_cast<ReorderedLinearOp>(prec_Op,true)->getBlockedOp();
00204         Teko::rebuildInverse(*invFactory_,blockedFwdOp,reorderedPrec);
00205      }
00206   }
00207   else {
00208      // no reordering required
00209      if(prec_Op==Teuchos::null)
00210         prec_Op = Teko::buildInverse(*invFactory_,fwdOp);
00211      else
00212         Teko::rebuildInverse(*invFactory_,fwdOp,prec_Op);
00213   }
00214 
00215   // construct preconditioner
00216   timer.stop();
00217 
00218   if(mediumVerbosity)
00219     Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
00220 
00221 
00222   // Initialize the preconditioner
00223   defaultPrec.initializeUnspecified( prec_Op);
00224   defaultPrec.incrIter();
00225   totalTimer.stop();
00226 
00227   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00228     *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00229   if(mediumVerbosity)
00230     *out << "\nLeaving Teko::StratimikosFactory::initializePrec_Thyra(...) ...\n";
00231 }
00232 
00233 void StratimikosFactory::initializePrec_Epetra(
00234   const Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > &fwdOpSrc,
00235   Thyra::PreconditionerBase<double> *prec,
00236   const Thyra::ESupportSolveUse supportSolveUse
00237   ) const
00238 {
00239   using Teuchos::outArg;
00240   using Teuchos::RCP;
00241   using Teuchos::rcp;
00242   using Teuchos::rcp_dynamic_cast;
00243   using Teuchos::implicit_cast;
00244 
00245   Teuchos::Time totalTimer(""), timer("");
00246   totalTimer.start(true);
00247 
00248   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00249   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00250 
00251   bool mediumVerbosity = (out.get() && implicit_cast<int>(verbLevel) > implicit_cast<int>(Teuchos::VERB_LOW));
00252 
00253   Teuchos::OSTab tab(out);
00254   if(mediumVerbosity)
00255     *out << "\nEntering Teko::StratimikosFactory::initializePrec_Epetra(...) ...\n";
00256 
00257   Teuchos::RCP<const LinearOpBase<double> > fwdOp = fwdOpSrc->getOp();
00258 #ifdef _DEBUG
00259   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00260   TEUCHOS_TEST_FOR_EXCEPT(prec==NULL);
00261 #endif
00262  
00263   //
00264   // Unwrap and get the forward Epetra_Operator object
00265   //
00266   Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
00267   Thyra::EOpTransp epetraFwdOpTransp;
00268   Thyra::EApplyEpetraOpAs epetraFwdOpApplyAs;
00269   Thyra::EAdjointEpetraOp epetraFwdOpAdjointSupport;
00270   double epetraFwdOpScalar;
00271   epetraFwdOpViewExtractor_->getEpetraOpView(
00272     fwdOp,outArg(epetraFwdOp),outArg(epetraFwdOpTransp),outArg(epetraFwdOpApplyAs),
00273     outArg(epetraFwdOpAdjointSupport),outArg(epetraFwdOpScalar)
00274                                              );
00275   // Get the concrete preconditioner object
00276   StratimikosFactoryPreconditioner & defaultPrec 
00277         = Teuchos::dyn_cast<StratimikosFactoryPreconditioner>(*prec);
00278 
00279   // Get the EpetraLinearOp object that is used to implement the preconditoner linear op
00280   RCP<Thyra::EpetraLinearOp> epetra_precOp 
00281         = rcp_dynamic_cast<Thyra::EpetraLinearOp>(defaultPrec.getNonconstUnspecifiedPrecOp(),true);
00282 
00283   // Get the embedded ML_Epetra::MultiLevelPreconditioner object if it exists
00284   Teuchos::RCP<Teko::Epetra::InverseFactoryOperator> teko_precOp;
00285   if(epetra_precOp!=Teuchos::null)
00286     teko_precOp = rcp_dynamic_cast<Teko::Epetra::InverseFactoryOperator>(epetra_precOp->epetra_op(),true);
00287 
00288   // Permform initialization if needed
00289   const bool startingOver = (teko_precOp == Teuchos::null);
00290   if(startingOver) {
00291     // start over
00292     invLib_ = Teuchos::null;
00293     invFactory_ = Teuchos::null;
00294     decomp_.clear();
00295 
00296     if(mediumVerbosity)
00297       *out << "\nCreating the initial Teko::Epetra::InverseFactoryOperator object...\n";
00298 
00299     timer.start(true);
00300 
00301     std::stringstream ss;
00302     ss << paramList_->get<std::string>("Strided Blocking");
00303 
00304     // figure out the decomposition requested by the string
00305     while(not ss.eof()) {
00306        int num=0;
00307        ss >> num;
00308 
00309        TEUCHOS_ASSERT(num>0);
00310 
00311        decomp_.push_back(num);
00312     }
00313 
00314     // build library, and set request handler (user defined!)
00315     invLib_  = Teko::InverseLibrary::buildFromParameterList(paramList_->sublist("Inverse Factory Library"));
00316     invLib_->setRequestHandler(reqHandler_);
00317 
00318     // build preconditioner factory
00319     invFactory_ = invLib_->getInverseFactory(paramList_->get<std::string>("Inverse Type"));
00320 
00321     // Create the initial preconditioner: DO NOT compute it yet
00322     teko_precOp = rcp( new Teko::Epetra::InverseFactoryOperator(invFactory_));
00323     
00324     timer.stop();
00325     if(mediumVerbosity)
00326       Teuchos::OSTab(out).o() <<"> Creation time = "<<timer.totalElapsedTime()<<" sec\n";
00327   }
00328 
00329   if(mediumVerbosity)
00330     *out << "\nComputing the preconditioner ...\n";
00331 
00332   // construct preconditioner
00333   {
00334      bool writeBlockOps = paramList_->get<bool>("Write Block Operator");
00335      timer.start(true);
00336      teko_precOp->initInverse();
00337 
00338      if(decomp_.size()==1) {
00339         teko_precOp->rebuildInverseOperator(epetraFwdOp);
00340 
00341         // write out to disk
00342         if(writeBlockOps) {
00343            std::stringstream ss;
00344            ss << "block-" << defaultPrec.getIter() << "_00.mm";
00345            EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*rcp_dynamic_cast<const Epetra_CrsMatrix>(epetraFwdOp));
00346         }
00347      }
00348      else {
00349         Teuchos::RCP<Epetra_Operator> wrappedFwdOp 
00350            = buildWrappedEpetraOperator(epetraFwdOp,teko_precOp->getNonconstForwardOp(),*out);
00351 
00352         // write out to disk
00353         if(writeBlockOps) {
00354            std::stringstream ss;
00355            ss << "block-" << defaultPrec.getIter();
00356            // Teuchos::RCP<Teko::Epetra::StridedEpetraOperator> stridedJac
00357            //       = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedFwdOp);
00358            Teuchos::RCP<Teko::Epetra::BlockedEpetraOperator> stridedJac
00359                  = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedFwdOp);
00360            if(stridedJac!=Teuchos::null) {
00361               // write out blocks of strided operator
00362               stridedJac->WriteBlocks(ss.str());
00363            }
00364            else TEUCHOS_ASSERT(false);
00365         }
00366    
00367         teko_precOp->rebuildInverseOperator(wrappedFwdOp);
00368      }
00369 
00370      timer.stop();
00371   }
00372 
00373   if(mediumVerbosity)
00374     Teuchos::OSTab(out).o() <<"=> Preconditioner construction time = "<<timer.totalElapsedTime()<<" sec\n";
00375 
00376   // Initialize the output EpetraLinearOp
00377   if(startingOver) {
00378     epetra_precOp = rcp(new Thyra::EpetraLinearOp);
00379   }
00380   epetra_precOp->initialize(teko_precOp,epetraFwdOpTransp
00381     ,Thyra::EPETRA_OP_APPLY_APPLY_INVERSE
00382     ,Thyra::EPETRA_OP_ADJOINT_UNSUPPORTED);
00383 
00384   // Initialize the preconditioner
00385   defaultPrec.initializeUnspecified( Teuchos::rcp_implicit_cast<LinearOpBase<double> >(epetra_precOp));
00386   defaultPrec.incrIter();
00387   totalTimer.stop();
00388 
00389   if(out.get() && implicit_cast<int>(verbLevel) >= implicit_cast<int>(Teuchos::VERB_LOW))
00390     *out << "\nTotal time in Teko::StratimikosFactory = "<<totalTimer.totalElapsedTime()<<" sec\n";
00391   if(mediumVerbosity)
00392     *out << "\nLeaving Teko::StratimikosFactory::initializePrec(...) ...\n";
00393 }
00394 
00395 
00396 void StratimikosFactory::uninitializePrec(
00397   Thyra::PreconditionerBase<double> *prec,
00398   Teuchos::RCP<const Thyra::LinearOpSourceBase<double> > *fwdOp,
00399   Thyra::ESupportSolveUse *supportSolveUse
00400   ) const
00401 {
00402   TEUCHOS_TEST_FOR_EXCEPT(true);
00403 }
00404 
00405 
00406 // Overridden from ParameterListAcceptor
00407 
00408 
00409 void StratimikosFactory::setParameterList(
00410   Teuchos::RCP<Teuchos::ParameterList> const& paramList
00411   )
00412 {
00413    TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
00414 
00415    paramList->validateParametersAndSetDefaults(*this->getValidParameters(),0);
00416    paramList_ = paramList;
00417 
00418 }
00419 
00420 
00421 Teuchos::RCP<Teuchos::ParameterList>
00422 StratimikosFactory::getNonconstParameterList()
00423 {
00424   return paramList_;
00425 }
00426 
00427 
00428 Teuchos::RCP<Teuchos::ParameterList>
00429 StratimikosFactory::unsetParameterList()
00430 {
00431   Teuchos::RCP<ParameterList> _paramList = paramList_;
00432   paramList_ = Teuchos::null;
00433   return _paramList;
00434 }
00435 
00436 
00437 Teuchos::RCP<const Teuchos::ParameterList>
00438 StratimikosFactory::getParameterList() const
00439 {
00440   return paramList_;
00441 }
00442 
00443 
00444 Teuchos::RCP<const Teuchos::ParameterList>
00445 StratimikosFactory::getValidParameters() const
00446 {
00447 
00448   using Teuchos::rcp;
00449   using Teuchos::tuple;
00450   using Teuchos::implicit_cast;
00451   using Teuchos::rcp_implicit_cast;
00452   typedef Teuchos::ParameterEntryValidator PEV;
00453 
00454   static RCP<const ParameterList> validPL;
00455 
00456   if(is_null(validPL)) {
00457 
00458     Teuchos::RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
00459 
00460     pl->set("Test Block Operator",false,
00461             "If Stratiikos/Teko is used to break an operator into its parts,\n"
00462             "then setting this parameter to true will compare applications of the\n"
00463             "segregated operator to the original operator.");
00464     pl->set("Write Block Operator",false,
00465             "Write out the segregated operator to disk with the name \"block-?_xx\"");
00466     pl->set("Strided Blocking","1",
00467             "Assuming that the user wants Strided blocking, break the operating into\n"
00468             "blocks. The syntax can be thought to be associated with the solution\n"
00469             "vector. For example if your variables are [u v w p T], and we want [u v w]\n"
00470             "blocked together, and p and T seperate then the relevant string is \"3 1 1\".\n"
00471             "Meaning put the first 3 unknowns per node together and sperate the v and w\n"
00472             "components.");
00473     pl->set("Reorder Type","",
00474             "This specifies how the blocks are reordered for use in the preconditioner.\n"
00475             "For example, assume the linear system is generated from 3D Navier-Stokes\n"
00476             "with an energy equation, yielding the unknowns [u v w p T]. If the\n"
00477             "\"Strided Blocking\" string is \"3 1 1\", then setting this parameter to\n"
00478             "\"[2 [0 1]]\" will reorder the blocked operator so its nested with the\n"
00479             "velocity and pressure forming an inner two-by-two block, and then the\n"
00480             "temperature unknowns forming a two-by-two system with the velocity-pressure\n"
00481             "block.");
00482     pl->set("Inverse Type","Amesos", 
00483             "The type of inverse operator the user wants. This can be one of the defaults\n"
00484             "from Stratimikos, or a Teko preconditioner defined in the\n"
00485             "\"Inverse Factory Library\".");
00486     pl->sublist("Inverse Factory Library",false,"Definition of Teko preconditioners."); 
00487 
00488     validPL = pl;
00489   }
00490 
00491   return validPL;
00492 }
00493 
00497 Teuchos::RCP<Epetra_Operator> StratimikosFactory::buildWrappedEpetraOperator(
00498                                                                  const Teuchos::RCP<const Epetra_Operator> & Jac,
00499                                                                  const Teuchos::RCP<Epetra_Operator> & wrapInput,
00500                                                                  std::ostream & out
00501                                                                  ) const
00502 {
00503    Teuchos::RCP<Epetra_Operator> wrappedOp = wrapInput;
00504 
00505 //    // initialize jacobian
00506 //    if(wrappedOp==Teuchos::null)
00507 //    {
00508 //       wrappedOp = Teuchos::rcp(new Teko::Epetra::StridedEpetraOperator(decomp_,Jac));
00509 // 
00510 //       // reorder the blocks if requested
00511 //       std::string reorderType = paramList_->get<std::string>("Reorder Type");
00512 //       if(reorderType!="") {
00513 //          RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
00514 // 
00515 //          // out << "Teko: Reordering = " << brm->toString() << std::endl;
00516 //          Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->Reorder(*brm);
00517 //       }
00518 //    }
00519 //    else {
00520 //       Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->RebuildOps();
00521 //    }
00522 // 
00523 //    // test blocked operator for correctness
00524 //    if(paramList_->get<bool>("Test Block Operator")) {
00525 //       bool result
00526 //          = Teuchos::rcp_dynamic_cast<Teko::Epetra::StridedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
00527 // 
00528 //       out << "Teko: Tested operator correctness:  " << (result ? "passed" : "FAILED!") << std::endl;
00529 //    }
00530 
00531    // initialize jacobian
00532    if(wrappedOp==Teuchos::null)
00533    {
00534       // build strided vector
00535       std::vector<std::vector<int> > vars;
00536       buildStridedVectors(*Jac,decomp_,vars);
00537       wrappedOp = Teuchos::rcp(new Teko::Epetra::BlockedEpetraOperator(vars,Jac));
00538 
00539       // reorder the blocks if requested
00540       std::string reorderType = paramList_->get<std::string>("Reorder Type");
00541       if(reorderType!="") {
00542          RCP<const Teko::BlockReorderManager> brm = Teko::blockedReorderFromString(reorderType);
00543 
00544          // out << "Teko: Reordering = " << brm->toString() << std::endl;
00545          Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->Reorder(*brm);
00546       }
00547    }
00548    else {
00549       Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->RebuildOps();
00550    }
00551 
00552    // test blocked operator for correctness
00553    if(paramList_->get<bool>("Test Block Operator")) {
00554       bool result
00555          = Teuchos::rcp_dynamic_cast<Teko::Epetra::BlockedEpetraOperator>(wrappedOp)->testAgainstFullOperator(600,1e-14);
00556 
00557       out << "Teko: Tested operator correctness:  " << (result ? "passed" : "FAILED!") << std::endl;
00558    }
00559 
00560    return wrappedOp;
00561 }
00562 
00563 std::string StratimikosFactory::description() const
00564 {
00565   std::ostringstream oss;
00566   oss << "Teko::StratimikosFactory";
00567   return oss.str();
00568 }
00569 
00570 void StratimikosFactory::buildStridedVectors(const Epetra_Operator & Jac,
00571                                              const std::vector<int> & decomp,
00572                                              std::vector<std::vector<int> > & vars) const
00573 {
00574    const Epetra_Map & rangeMap = Jac.OperatorRangeMap();
00575 
00576    // compute total number of variables
00577    int numVars = 0;
00578    for(std::size_t i=0;i<decomp.size();i++)
00579       numVars += decomp[i];
00580 
00581    // verify that the decomposition is appropriate for this matrix
00582    TEUCHOS_ASSERT((rangeMap.NumMyElements() % numVars)==0);
00583    TEUCHOS_ASSERT((rangeMap.NumGlobalElements() % numVars)==0);
00584   
00585    int * globalIds = rangeMap.MyGlobalElements();
00586 
00587    vars.resize(decomp.size());
00588    for(int i=0;i<rangeMap.NumMyElements();) {
00589 
00590       // for each "node" copy global ids to vectors
00591       for(std::size_t d=0;d<decomp.size();d++) {
00592          // for this variable copy global ids to variable arrays
00593          int current = decomp[d];
00594          for(int v=0;v<current;v++,i++)
00595             vars[d].push_back(globalIds[i]);
00596       }
00597    }
00598 }
00599 
00600 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
00601                                const std::string & stratName)
00602 {
00603    TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
00604                       "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
00605 
00606    // use default constructor to add Teko::StratimikosFactory
00607    builder.setPreconditioningStrategyFactory(
00608          Teuchos::abstractFactoryStd<Thyra::PreconditionerFactoryBase<double>,Teko::StratimikosFactory>(),
00609          stratName);
00610 }
00611 
00612 void addTekoToStratimikosBuilder(Stratimikos::DefaultLinearSolverBuilder & builder,
00613                                const Teuchos::RCP<Teko::RequestHandler> & rh,
00614                                const std::string & stratName)
00615 {
00616    TEUCHOS_TEST_FOR_EXCEPTION(builder.getValidParameters()->sublist("Preconditioner Types").isParameter(stratName),std::logic_error,
00617                       "Teko::addTekoToStratimikosBuilder cannot add \"" + stratName +"\" because it is already included in builder!");
00618 
00619    // build an instance of a Teuchos::AbsractFactory<Thyra::PFB> so request handler is passed onto
00620    // the resulting StratimikosFactory
00621    Teuchos::RCP<TekoFactoryBuilder> tekoFactoryBuilder = Teuchos::rcp(new TekoFactoryBuilder(rh));
00622    builder.setPreconditioningStrategyFactory(tekoFactoryBuilder,stratName);
00623 }
00624 
00625 } // namespace Thyra
 All Classes Files Functions Variables