Thyra_AztecOOLinearOpWithSolveFactory.cpp

Go to the documentation of this file.
00001 /*@Header
00002 // ***********************************************************************
00003 // 
00004 //        AztecOO: An Object-Oriented Aztec Linear Solver Package 
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef __sun
00031 
00032 #include "Thyra_AztecOOLinearOpWithSolveFactory.hpp"
00033 #include "Thyra_AztecOOLinearOpWithSolve.hpp"
00034 #include "Thyra_PreconditionerFactoryHelpers.hpp"
00035 #include "Thyra_EpetraOperatorViewExtractorStd.hpp"
00036 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00037 #include "Thyra_EpetraLinearOpBase.hpp"
00038 #include "Thyra_EpetraOperatorWrapper.hpp"
00039 #include "EpetraExt_ProductOperator.h"
00040 #include "Teuchos_ParameterList.hpp"
00041 #include "Teuchos_dyn_cast.hpp"
00042 #include "AztecOOParameterList.hpp"
00043 
00044 namespace {
00045 
00046 const std::string AOOLOWSF_epetraPrecOp_str = "AOOLOWSF::epetraPrecOp";
00047 const std::string AOOLOWSF_aztec_epetra_epetraFwdOp_str = "AOOLOWSF::aztec_epetra_epetraFwdOp";
00048 const std::string AOOLOWSF_aztec_epetra_epetraAdjOp_str = "AOOLOWSF::aztec_epetra_epetraAdjOp";
00049 const std::string AOOLOWSF_rowmatrix_epetraFwdOp_str = "AOOLOWSF::rowmatrix_epetraFwdOp";
00050 const std::string AOOLOWSF_rowmatrix_epetraPrecOp_str = "AOOLOWSF::rowmatrix_epetraPrecOp";
00051 const std::string AOOLOWSF_aztec_fwd_epetra_epetraPrecOp_str = "AOOLOWSF::aztec_fwd_epetra_epetraPrecOp";
00052 const std::string AOOLOWSF_aztec_adj_epetra_epetraPrecOp_str = "AOOLOWSF::aztec_adj_epetra_epetraPrecOp";
00053 const std::string AOOLOWSF_setPrecondtionerOperator_str = "AOOLOWSF::setPrecondtionerOperator";
00054 const std::string AOOLOWSF_constructedAztecPreconditoner_str = "AOOLOWSF::constructedAztecPreconditoner";
00055 
00056 const std::string  ForwardSolve_name = "Forward Solve";
00057 const std::string  AdjointSolve_name = "Adjoint Solve";
00058 const std::string  MaxIterations_name = "Max Iterations";
00059 const int          MaxIterations_default = 400;
00060 const std::string  Tolerance_name = "Tolerance";
00061 const double       Tolerance_default = 1e-6;
00062 const std::string  OutputEveryRhs_name = "Output Every RHS";
00063 const bool         OutputEveryRhs_default = false;
00064 const std::string  AztecOO_Settings_name = "AztecOO Settings";
00065 
00066 } // namespace
00067 
00068 namespace Thyra {
00069 
00070 // Constructors/initializers/accessors
00071 
00072 AztecOOLinearOpWithSolveFactory::AztecOOLinearOpWithSolveFactory()
00073   :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
00074   ,defaultFwdMaxIterations_(MaxIterations_default)
00075   ,defaultFwdTolerance_(Tolerance_default)
00076   ,defaultAdjMaxIterations_(MaxIterations_default)
00077   ,defaultAdjTolerance_(Tolerance_default)
00078   ,outputEveryRhs_(OutputEveryRhs_default)
00079 {
00080   updateThisValidParamList();
00081 }
00082 
00083 // Overridden from LinearOpWithSolveFactoryBase
00084 
00085 bool AztecOOLinearOpWithSolveFactory::acceptsPreconditionerFactory() const
00086 {
00087   return true;
00088 }
00089 
00090 void AztecOOLinearOpWithSolveFactory::setPreconditionerFactory(
00091   const Teuchos::RefCountPtr<PreconditionerFactoryBase<double> >  &precFactory
00092   ,const std::string                                              &precFactoryName
00093   )
00094 {
00095   TEST_FOR_EXCEPT(!precFactory.get());
00096   Teuchos::RefCountPtr<const Teuchos::ParameterList>
00097     precFactoryValidPL = precFactory->getValidParameters();
00098   const std::string _precFactoryName =
00099     ( precFactoryName != ""
00100       ? precFactoryName
00101       : ( precFactoryValidPL.get() ? precFactoryValidPL->name() : "GENERIC PRECONDITIONER FACTORY" )
00102       );
00103   precFactory_ = precFactory;
00104   precFactoryName_ = _precFactoryName;
00105   updateThisValidParamList();
00106 }
00107 
00108 Teuchos::RefCountPtr<PreconditionerFactoryBase<double> >
00109 AztecOOLinearOpWithSolveFactory::getPreconditionerFactory() const
00110 {
00111   return precFactory_;
00112 }
00113 
00114 void AztecOOLinearOpWithSolveFactory::unsetPreconditionerFactory(
00115   Teuchos::RefCountPtr<PreconditionerFactoryBase<double> >  *precFactory
00116   ,std::string                                              *precFactoryName
00117   )
00118 {
00119   if(precFactory) *precFactory = precFactory_;
00120   if(precFactoryName) *precFactoryName = precFactoryName_;
00121   precFactory_ = Teuchos::null;
00122   precFactoryName_ = "";
00123   updateThisValidParamList();
00124 }
00125 
00126 bool AztecOOLinearOpWithSolveFactory::isCompatible(
00127   const LinearOpSourceBase<double> &fwdOpSrc
00128   ) const
00129 {
00130   return epetraFwdOpViewExtractor_->isCompatible(*fwdOpSrc.getOp());
00131 }
00132 
00133 Teuchos::RefCountPtr<LinearOpWithSolveBase<double> >
00134 AztecOOLinearOpWithSolveFactory::createOp() const
00135 {
00136   return Teuchos::rcp(new AztecOOLinearOpWithSolve());
00137 }
00138 
00139 void AztecOOLinearOpWithSolveFactory::initializeOp(
00140   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >    &fwdOpSrc
00141   ,LinearOpWithSolveBase<double>                                   *Op
00142   ,const ESupportSolveUse                                          supportSolveUse
00143   ) const
00144 {
00145   this->initializeOp_impl(fwdOpSrc,Teuchos::null,Teuchos::null,false,Op);
00146 }
00147 
00148 void AztecOOLinearOpWithSolveFactory::initializeAndReuseOp(
00149   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >    &fwdOpSrc
00150   ,LinearOpWithSolveBase<double>                                   *Op
00151   ) const
00152 {
00153   this->initializeOp_impl(fwdOpSrc,Teuchos::null,Teuchos::null,true,Op);
00154 }
00155 
00156 bool AztecOOLinearOpWithSolveFactory::supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
00157 {
00158   const_cast<bool&>(useAztecPrec_) = (
00159     paramList_.get() && paramList_->sublist(ForwardSolve_name).sublist(AztecOO_Settings_name).get("Aztec Preconditioner","none")!="none"
00160     );
00161   switch(precOpType) {
00162     case PRECONDITIONER_INPUT_TYPE_AS_OPERATOR:
00163       return true;
00164       break;
00165     case PRECONDITIONER_INPUT_TYPE_AS_MATRIX:
00166       return useAztecPrec_;
00167       break;
00168     default:
00169       TEST_FOR_EXCEPT(true);
00170   }
00171   return PRECONDITIONER_INPUT_TYPE_AS_OPERATOR; // Should never be called!
00172 }
00173 
00174 void AztecOOLinearOpWithSolveFactory::initializePreconditionedOp(
00175   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >       &fwdOpSrc
00176   ,const Teuchos::RefCountPtr<const PreconditionerBase<double> >      &prec
00177   ,LinearOpWithSolveBase<double>                                      *Op
00178   ,const ESupportSolveUse                                             supportSolveUse
00179   ) const
00180 {
00181   TEST_FOR_EXCEPT(prec.get()==NULL);
00182   this->initializeOp_impl(fwdOpSrc,prec,Teuchos::null,false,Op);
00183 }
00184 
00185 void AztecOOLinearOpWithSolveFactory::initializeApproxPreconditionedOp(
00186   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >       &fwdOpSrc
00187   ,const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >      &approxFwdOpSrc
00188   ,LinearOpWithSolveBase<double>                                      *Op
00189   ,const ESupportSolveUse                                             supportSolveUse
00190   ) const
00191 {
00192   TEST_FOR_EXCEPT(approxFwdOpSrc.get()==NULL);
00193   TEST_FOR_EXCEPT(approxFwdOpSrc->getOp().get()==NULL);
00194   this->initializeOp_impl(fwdOpSrc,Teuchos::null,approxFwdOpSrc,false,Op);
00195 }
00196 
00197 void AztecOOLinearOpWithSolveFactory::uninitializeOp(
00198   LinearOpWithSolveBase<double>                               *Op
00199   ,Teuchos::RefCountPtr<const LinearOpSourceBase<double> >    *fwdOpSrc
00200   ,Teuchos::RefCountPtr<const PreconditionerBase<double> >    *prec
00201   ,Teuchos::RefCountPtr<const LinearOpSourceBase<double> >    *approxFwdOpSrc
00202   ,ESupportSolveUse                                           *supportSolveUse
00203   ) const
00204 {
00205 #ifdef TEUCHOS_DEBUG
00206   TEST_FOR_EXCEPT(Op==NULL);
00207 #endif
00208   AztecOOLinearOpWithSolve
00209     *aztecOp = &Teuchos::dyn_cast<AztecOOLinearOpWithSolve>(*Op);
00210   // Extract and unset the fwdOP and approxFwdOp objects
00211   Teuchos::RefCountPtr<const LinearOpSourceBase<double> >
00212     _fwdOpSrc = aztecOp->extract_fwdOpSrc(),             // Will be null if not initialized!
00213     _approxFwdOpSrc = aztecOp->extract_approxFwdOpSrc(); // Will be null if no approxFwdOp set
00214   if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc;
00215   if(approxFwdOpSrc) *approxFwdOpSrc = _approxFwdOpSrc;
00216   // Only extract and uset the prec object if it is external.  If it is
00217   // internal, then we need to hold on to this so that we can reinitialize it
00218   // later.
00219   if(aztecOp->isExternalPrec()) {
00220     Teuchos::RefCountPtr<const PreconditionerBase<double> >
00221       _prec = aztecOp->extract_prec(); // Will be null if not external preconditioner was set
00222     if(prec) *prec = _prec;
00223   }
00224   // ToDo: Extract the Epetra_Operator views what where used to initialize the
00225   // forward and adjoint solvers!  This is needed to make this totally
00226   // stateless.
00227 }
00228 
00229 // Overridden from ParameterListAcceptor
00230 
00231 void AztecOOLinearOpWithSolveFactory::setParameterList(Teuchos::RefCountPtr<Teuchos::ParameterList> const& paramList)
00232 {
00233   TEST_FOR_EXCEPT(paramList.get()==NULL);
00234   paramList->validateParameters(*this->getValidParameters());
00235   paramList_ = paramList;
00236   //
00237   outputEveryRhs_ = paramList_->get(OutputEveryRhs_name,OutputEveryRhs_default);
00238   // Foward Solve parameters
00239   Teuchos::ParameterList
00240     &fwdSolvePL = paramList_->sublist(ForwardSolve_name);
00241   defaultFwdMaxIterations_ = fwdSolvePL.get(MaxIterations_name,defaultFwdMaxIterations_);
00242   defaultFwdTolerance_ = fwdSolvePL.get(Tolerance_name,defaultFwdTolerance_);
00243   // Adjoint Solve parameters
00244   if( !paramList_->getPtr<Teuchos::ParameterList>(AdjointSolve_name) ) {
00245     // If adjoint solve sublist is not set, then use the forward solve parameters
00246     paramList_->sublist(AdjointSolve_name).setParameters(fwdSolvePL);
00247   }
00248   Teuchos::ParameterList
00249     &adjSolvePL = paramList_->sublist(AdjointSolve_name);
00250   defaultAdjMaxIterations_ = adjSolvePL.get(MaxIterations_name,defaultAdjMaxIterations_);
00251   defaultAdjTolerance_ = adjSolvePL.get(Tolerance_name,defaultAdjTolerance_);
00252   //
00253   if(precFactory_.get()) {
00254     // Only reset the PF's PL if the sublist exists or the PF does ot already
00255     // have a PL.  We don't want to overwrite an externally set PL for the PF
00256     // if we don't have a nested sublist defined here!
00257     const bool nestedPFSublistExists = paramList_->isSublist(precFactoryName_);
00258     const bool alreadyHasSublist = !is_null(precFactory_->getParameterList());
00259     if( nestedPFSublistExists || !alreadyHasSublist ) {
00260       precFactory_->setParameterList(Teuchos::sublist(paramList_,precFactoryName_));
00261     }
00262   }
00263 }
00264 
00265 Teuchos::RefCountPtr<Teuchos::ParameterList>
00266 AztecOOLinearOpWithSolveFactory::getParameterList()
00267 {
00268   return paramList_;
00269 }
00270 
00271 Teuchos::RefCountPtr<Teuchos::ParameterList>
00272 AztecOOLinearOpWithSolveFactory::unsetParameterList()
00273 {
00274   Teuchos::RefCountPtr<Teuchos::ParameterList> _paramList = paramList_;
00275   paramList_ = Teuchos::null;
00276   return _paramList;
00277 }
00278 
00279 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00280 AztecOOLinearOpWithSolveFactory::getParameterList() const
00281 {
00282   return paramList_;
00283 }
00284 
00285 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00286 AztecOOLinearOpWithSolveFactory::getValidParameters() const
00287 {
00288   return thisValidParamList_;
00289 }
00290 
00291 // Public functions overridden from Teuchos::Describable
00292 
00293 std::string AztecOOLinearOpWithSolveFactory::description() const
00294 {
00295   std::ostringstream oss;
00296   oss << "Thyra::AztecOOLinearOpWithSolveFactory";
00297   return oss.str();
00298 }
00299 
00300 // private
00301 
00302 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00303 AztecOOLinearOpWithSolveFactory::generateAndGetValidParameters()
00304 {
00305   static Teuchos::RefCountPtr<Teuchos::ParameterList> validParamList;
00306   if(validParamList.get()==NULL) {
00307     validParamList = Teuchos::rcp(new Teuchos::ParameterList("AztecOOLinearOpWithSolveFactory"));
00308     validParamList->set(OutputEveryRhs_name,OutputEveryRhs_default);
00309     static Teuchos::RefCountPtr<const Teuchos::ParameterList>
00310       aztecParamList = getValidAztecOOParameters();
00311     Teuchos::ParameterList
00312       &fwdSolvePL = validParamList->sublist(ForwardSolve_name);
00313     fwdSolvePL.set(Tolerance_name,Tolerance_default);
00314     fwdSolvePL.set(MaxIterations_name,MaxIterations_default);
00315     fwdSolvePL.sublist(AztecOO_Settings_name).setParameters(*aztecParamList);
00316     Teuchos::ParameterList
00317       &adjSolvePL = validParamList->sublist(AdjointSolve_name);
00318     adjSolvePL.setParameters(fwdSolvePL); // Make the adjoint solve have same defaults as forward solve
00319   }
00320   return validParamList;
00321 }
00322 
00323 void AztecOOLinearOpWithSolveFactory::updateThisValidParamList()
00324 {
00325   thisValidParamList_ = Teuchos::rcp(
00326     new Teuchos::ParameterList(*generateAndGetValidParameters())
00327     );
00328   if(precFactory_.get()) {
00329     Teuchos::RefCountPtr<const Teuchos::ParameterList>
00330       precFactoryValidParamList = precFactory_->getValidParameters();
00331     if(precFactoryValidParamList.get()) {
00332       thisValidParamList_->sublist(precFactoryName_).setParameters(*precFactoryValidParamList);
00333     }
00334   }
00335 }
00336 
00337 void AztecOOLinearOpWithSolveFactory::initializeOp_impl(
00338   const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >       &fwdOpSrc
00339   ,const Teuchos::RefCountPtr<const PreconditionerBase<double> >      &prec
00340   ,const Teuchos::RefCountPtr<const LinearOpSourceBase<double> >      &approxFwdOpSrc
00341   ,const bool                                                         reusePrec
00342   ,LinearOpWithSolveBase<double>                                      *Op
00343   ) const
00344 {
00345   using Teuchos::RefCountPtr;
00346   using Teuchos::null;
00347   using Teuchos::rcp;
00348   using Teuchos::rcp_dynamic_cast;
00349   using Teuchos::rcp_const_cast;
00350   using Teuchos::set_extra_data;
00351   using Teuchos::get_optional_extra_data;
00352   typedef EpetraExt::ProductOperator PO;
00353 
00354   const Teuchos::RefCountPtr<Teuchos::FancyOStream> out       = this->getOStream();
00355   const Teuchos::EVerbosityLevel                    verbLevel = this->getVerbLevel();
00356   Teuchos::OSTab tab(out);
00357   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00358     *out << "\nEntering Thyra::AztecOOLinearOpWithSolveFactory::initializeOp_impl(...) ...\n";
00359 
00360   typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<double> > VOTSPF;
00361   VOTSPF precFactoryOutputTempState(precFactory_,out,verbLevel);
00362 
00363 #ifdef TEUCHOS_DEBUG
00364   TEST_FOR_EXCEPT(Op==NULL);
00365   TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00366   TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
00367 #endif
00368   Teuchos::RefCountPtr<const LinearOpBase<double> >
00369     tmpFwdOp = fwdOpSrc->getOp(),
00370     tmpApproxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
00371   Teuchos::RefCountPtr<const LinearOpBase<double> > fwdOp;
00372   Teuchos::RefCountPtr<const LinearOpBase<double> > approxFwdOp;
00373 
00374   // 
00375   // Determine whether the operators are EpetraLinearOp objects. If so, we're
00376   // good to go.  If not, we need to wrap it.
00377   //
00378 
00379   if ( dynamic_cast<const EpetraLinearOpBase*>(tmpFwdOp.get())!=0 )
00380   {
00381     fwdOp = tmpFwdOp;
00382     approxFwdOp = tmpApproxFwdOp;
00383   }
00384   else
00385   {
00386     fwdOp = makeEpetraWrapper(ConstLinearOperator<double>(tmpFwdOp));
00387     if (tmpApproxFwdOp.get() && dynamic_cast<const EpetraLinearOpBase*>(&*tmpApproxFwdOp.get()) )
00388       approxFwdOp = makeEpetraWrapper(ConstLinearOperator<double>(tmpApproxFwdOp));
00389   }
00390 
00391   //
00392   // Get the AztecOOLinearOpWithSolve object
00393   //
00394   AztecOOLinearOpWithSolve
00395     *aztecOp = &Teuchos::dyn_cast<AztecOOLinearOpWithSolve>(*Op);
00396 
00397   //
00398   // Unwrap and get the forward operator or matrix
00399   //
00400   Teuchos::RefCountPtr<const Epetra_Operator> epetra_epetraFwdOp;
00401   ETransp                                     epetra_epetraFwdOpTransp;
00402   EApplyEpetraOpAs                            epetra_epetraFwdOpApplyAs;
00403   EAdjointEpetraOp                            epetra_epetraFwdOpAdjointSupport;
00404   double                                      epetra_epetraFwdOpScalar;
00405   epetraFwdOpViewExtractor_->getEpetraOpView(
00406     fwdOp
00407     ,&epetra_epetraFwdOp,&epetra_epetraFwdOpTransp,&epetra_epetraFwdOpApplyAs
00408     ,&epetra_epetraFwdOpAdjointSupport,&epetra_epetraFwdOpScalar
00409     );
00410   TEST_FOR_EXCEPTION(
00411     epetra_epetraFwdOp.get()==NULL, std::logic_error
00412     ,"Error, The input fwdOp object must be fully initialized before calling this function!"
00413     );
00414   //
00415   // Get the preconditioner object to use
00416   //
00417   Teuchos::RefCountPtr<PreconditionerBase<double> >        myPrec;
00418   Teuchos::RefCountPtr<const PreconditionerBase<double> >  precUsed;
00419   if(prec.get()) {
00420     // We will be used the passed in external preconditioner
00421     precUsed = prec;
00422   }
00423   else if(precFactory_.get() ) {
00424     // We will be creating our own preconditioner using a preconditioner factory
00425     myPrec =
00426       ( !aztecOp->isExternalPrec()
00427         ? Teuchos::rcp_const_cast<PreconditionerBase<double> >(aztecOp->extract_prec())
00428         : Teuchos::null
00429         );
00430     if(myPrec.get()) {
00431       // ToDo: Get the forward operator and validate that it is the same
00432       // operator that is used here!
00433     }
00434     else {
00435       myPrec = precFactory_->createPrec();
00436     }
00437     precFactory_->initializePrec(fwdOpSrc,&*myPrec);
00438     precUsed = myPrec;
00439   }
00440   //
00441   // Unwrap and get the preconditioner operator
00442   //
00443   RefCountPtr<const LinearOpBase<double> > rightPrecOp;
00444   if(precUsed.get()) {
00445     RefCountPtr<const LinearOpBase<double> > unspecified = precUsed->getUnspecifiedPrecOp();
00446     RefCountPtr<const LinearOpBase<double> > left        = precUsed->getLeftPrecOp();
00447     RefCountPtr<const LinearOpBase<double> > right       = precUsed->getRightPrecOp();
00448     TEST_FOR_EXCEPTION(
00449       !( left.get() || right.get() || unspecified.get() ), std::logic_error
00450       ,"Error, at least one preconditoner linear operator objects must be set!"
00451       );
00452     if(unspecified.get()) {
00453       rightPrecOp = unspecified;
00454     }
00455     else {
00456       // Set a left, right or split preconditioner
00457       TEST_FOR_EXCEPTION(
00458         left.get(),std::logic_error
00459         ,"Error, we can not currently handle a left preconditioner with the Thyra/AztecOO adapters!"
00460         );
00461       rightPrecOp = right;
00462     }
00463   }
00464   double                                       wrappedPrecOpScalar = 0.0;
00465   ETransp                                      wrappedPrecOpTransp = NOTRANS;
00466   RefCountPtr<const LinearOpBase<double> >     wrappedPrecOp = null;
00467   RefCountPtr<const EpetraLinearOpBase>        epetraPrecOp;
00468   Teuchos::RefCountPtr<const Epetra_Operator>  epetra_epetraPrecOp;
00469   ETransp                                      epetra_epetraPrecOpTransp;
00470   EApplyEpetraOpAs                             epetra_epetraPrecOpApplyAs;
00471   EAdjointEpetraOp                             epetra_epetraPrecOpAdjointSupport;
00472   ETransp                                      overall_epetra_epetraPrecOpTransp;
00473   if(rightPrecOp.get()) {
00474     RefCountPtr<const LinearOpBase<double> > tmpWrappedPrecOp;
00475     unwrap(rightPrecOp,&wrappedPrecOpScalar,&wrappedPrecOpTransp,&tmpWrappedPrecOp);
00476     if( dynamic_cast<const EpetraLinearOpBase*>(&*tmpWrappedPrecOp) ) {
00477       wrappedPrecOp = tmpWrappedPrecOp;
00478     }
00479     else {
00480       wrappedPrecOp = makeEpetraWrapper(ConstLinearOperator<double>(tmpWrappedPrecOp));
00481     }
00482     epetraPrecOp = rcp_dynamic_cast<const EpetraLinearOpBase>(wrappedPrecOp,true);
00483     epetraPrecOp->getEpetraOpView(&epetra_epetraPrecOp,&epetra_epetraPrecOpTransp,&epetra_epetraPrecOpApplyAs,&epetra_epetraPrecOpAdjointSupport);
00484     TEST_FOR_EXCEPTION(
00485       epetra_epetraPrecOp.get()==NULL,std::logic_error
00486       ,"Error, The input prec object and its embedded precondtioner operator must be fully initialized before calling this function!"
00487       );
00488     set_extra_data(epetraPrecOp,AOOLOWSF_epetraPrecOp_str,&epetra_epetraPrecOp,Teuchos::POST_DESTROY,false);
00489     overall_epetra_epetraPrecOpTransp
00490       = trans_trans(real_trans(wrappedPrecOpTransp),real_trans(epetra_epetraPrecOpTransp));
00491   }
00492   //
00493   // Unwrap and get the approximate forward operator to be used to generate a preconditioner
00494   //
00495   if(approxFwdOp.get()) {
00496     // Note, here we just use the same members data that would be set for an
00497     // extenral preconditioner operator since it is not getting used.
00498     unwrap(approxFwdOp,&wrappedPrecOpScalar,&wrappedPrecOpTransp,&wrappedPrecOp);
00499     epetraPrecOp = rcp_dynamic_cast<const EpetraLinearOpBase>(wrappedPrecOp,true);
00500     epetraPrecOp->getEpetraOpView(&epetra_epetraPrecOp,&epetra_epetraPrecOpTransp,&epetra_epetraPrecOpApplyAs,&epetra_epetraPrecOpAdjointSupport);
00501     TEST_FOR_EXCEPTION(
00502       epetra_epetraPrecOp.get()==NULL,std::logic_error
00503       ,"Error, The input approxFwdOp object must be fully initialized before calling this function!"
00504       );
00505     set_extra_data(epetraPrecOp,AOOLOWSF_epetraPrecOp_str,&epetra_epetraPrecOp,Teuchos::POST_DESTROY,false);
00506     overall_epetra_epetraPrecOpTransp
00507       = trans_trans(real_trans(wrappedPrecOpTransp),real_trans(epetra_epetraPrecOpTransp));
00508   }
00509   //
00510   // Determine if the forward and preconditioner operators are a row matrices or not
00511   //
00512   RefCountPtr<const Epetra_RowMatrix>
00513     rowmatrix_epetraFwdOp  = rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_epetraFwdOp),
00514     rowmatrix_epetraPrecOp = rcp_dynamic_cast<const Epetra_RowMatrix>(epetra_epetraPrecOp);
00515   //
00516   // Determine the type of preconditoner
00517   //
00518   this->supportsPreconditionerInputType(PRECONDITIONER_INPUT_TYPE_AS_MATRIX); // Updates useAztecPrec_, input value does not matter
00519   enum ELocalPrecType { PT_NONE, PT_AZTEC_FROM_OP, PT_AZTEC_FROM_APPROX_FWD_MATRIX, PT_FROM_PREC_OP };
00520   ELocalPrecType localPrecType;
00521   if( precUsed.get()==NULL && approxFwdOp.get()==NULL && !useAztecPrec_ ) {
00522     // No preconditioning at all!
00523     localPrecType = PT_NONE;
00524   }
00525   else if( precUsed.get()==NULL && approxFwdOp.get()==NULL && useAztecPrec_ ) {
00526     // We are using the forward matrix for the preconditioner using aztec preconditioners
00527     localPrecType = PT_AZTEC_FROM_OP;
00528   }
00529   else if( approxFwdOp.get() && useAztecPrec_ ) {
00530     // The preconditioner comes from the input as a matrix and we are using aztec preconditioners
00531     localPrecType = PT_AZTEC_FROM_APPROX_FWD_MATRIX;
00532   }
00533   else if( precUsed.get() ) {
00534     // The preconditioner comes as an external operator so let's use it as such
00535     localPrecType = PT_FROM_PREC_OP;
00536   }
00537   //
00538   // Determine if aztecOp already contains solvers and if we need to reinitialize or not
00539   //
00540   RefCountPtr<AztecOO> aztecFwdSolver, aztecAdjSolver;
00541   bool startingOver;
00542   if(1){
00543     // Let's assume that fwdOp, prec and/or approxFwdOp are compatible with
00544     // the already created AztecOO objects.  If they are not, then the client
00545     // should have created a new LOWSB object from scratch!
00546     Teuchos::RefCountPtr<const LinearOpBase<double> >        old_fwdOp;
00547     Teuchos::RefCountPtr<const LinearOpSourceBase<double> >  old_fwdOpSrc;
00548     Teuchos::RefCountPtr<const PreconditionerBase<double> >  old_prec;
00549     bool                                                     old_isExternalPrec;
00550     Teuchos::RefCountPtr<const LinearOpSourceBase<double> >  old_approxFwdOpSrc;
00551     Teuchos::RefCountPtr<AztecOO>                            old_aztecFwdSolver;
00552     Teuchos::RefCountPtr<AztecOO>                            old_aztecAdjSolver;
00553     double                                                   old_aztecSolverScalar;
00554     aztecOp->uninitialize(
00555       &old_fwdOp
00556       ,&old_fwdOpSrc
00557       ,&old_prec
00558       ,&old_isExternalPrec
00559       ,&old_approxFwdOpSrc
00560       ,&old_aztecFwdSolver
00561       ,NULL
00562       ,&old_aztecAdjSolver
00563       ,NULL
00564       ,&old_aztecSolverScalar
00565       );
00566     if( old_aztecFwdSolver.get()==NULL ) {
00567       // This has never been initialized before
00568       startingOver = true;
00569     }
00570     else {
00571       // Let's assume that fwdOp, prec and/or approxFwdOp are compatible with
00572       // the already created AztecOO objects.  If they are not, then the
00573       // client should have created a new LOWSB object from scratch!
00574       aztecFwdSolver = old_aztecFwdSolver;
00575       aztecAdjSolver = old_aztecAdjSolver;
00576       startingOver = false;
00577       // We must wipe out the old preconditoner if we are not reusing the preconditioner
00578       bool *constructedAztecPreconditioner = NULL;
00579       if(
00580         !reusePrec
00581         && ( constructedAztecPreconditioner = get_optional_extra_data<bool>(aztecFwdSolver,"AOOLOWSF::constructedAztecPreconditoner") )
00582         && *constructedAztecPreconditioner
00583         )
00584       {
00585         aztecFwdSolver->DestroyPreconditioner();
00586         *constructedAztecPreconditioner = false;
00587       }
00588       // We must see if we set an external preconditioner but will not do so
00589       // again in which case we must blow away AztecOO and start over again!
00590       bool *setPreconditionerOperator = NULL;
00591       if(
00592         localPrecType != PT_FROM_PREC_OP
00593         && ( setPreconditionerOperator = get_optional_extra_data<bool>(aztecFwdSolver,"AOOLOWSF::setPreconditonerOperator") )
00594         && *setPreconditionerOperator
00595         )
00596       {
00597         // We must start over again since there is no way to unset an external preconditioner!
00598         startingOver = true;
00599       }
00600     }
00601   }
00602   //
00603   // Create the AztecOO solvers if we are starting over
00604   //
00605   startingOver = true; // ToDo: Remove this and figure out why this is not working!
00606   if(startingOver) {
00607     // Forward solver
00608     aztecFwdSolver = rcp(new AztecOO());
00609     //aztecFwdSolver->SetAztecOption(AZ_output,AZ_none); // Don't mess up output
00610     //aztecFwdSolver->SetAztecOption(AZ_conv,AZ_rhs);    // Specified by this interface (may change)
00611     aztecFwdSolver->SetAztecOption(AZ_diagnostics,AZ_none); // This was turned off in NOX?
00612     aztecFwdSolver->SetAztecOption(AZ_keep_info,1);
00613     // Adjoint solver (if supported)
00614     if(
00615       epetra_epetraFwdOpAdjointSupport==EPETRA_OP_ADJOINT_SUPPORTED
00616       && localPrecType!=PT_AZTEC_FROM_OP && localPrecType!=PT_AZTEC_FROM_APPROX_FWD_MATRIX
00617       //&& (epetra_epetraPrecOp.get()==NULL ||epetra_epetraPrecOpAdjointSupport==EPETRA_OP_ADJOINT_SUPPORTED)
00618       )
00619     {
00620       aztecAdjSolver = rcp(new AztecOO());
00621       //aztecAdjSolver->SetAztecOption(AZ_output,AZ_none);
00622       //aztecAdjSolver->SetAztecOption(AZ_conv,AZ_rhs);
00623       aztecAdjSolver->SetAztecOption(AZ_diagnostics,AZ_none);
00624       //aztecAdjSolver->SetAztecOption(AZ_keep_info,1);
00625     }
00626   }
00627   //
00628   // Set the options on the AztecOO solvers
00629   //
00630   if( startingOver ) {
00631     if(paramList_.get())
00632       setAztecOOParameters(&paramList_->sublist(ForwardSolve_name).sublist(AztecOO_Settings_name),&*aztecFwdSolver);
00633     if(aztecAdjSolver.get() && paramList_.get())
00634       setAztecOOParameters(&paramList_->sublist(AdjointSolve_name).sublist(AztecOO_Settings_name),&*aztecAdjSolver);
00635   }
00636   //
00637   // Process the forward operator
00638   //
00639   RefCountPtr<const Epetra_Operator>
00640     aztec_epetra_epetraFwdOp,
00641     aztec_epetra_epetraAdjOp;
00642   // Forward solve
00643   RefCountPtr<const Epetra_Operator>
00644     epetraOps[]
00645     = { epetra_epetraFwdOp };
00646   Teuchos::ETransp
00647     epetraOpsTransp[]
00648     = { epetra_epetraFwdOpTransp==NOTRANS ? Teuchos::NO_TRANS : Teuchos::TRANS };
00649   PO::EApplyMode
00650     epetraOpsApplyMode[]
00651     = { epetra_epetraFwdOpApplyAs==EPETRA_OP_APPLY_APPLY ? PO::APPLY_MODE_APPLY : PO::APPLY_MODE_APPLY_INVERSE };
00652   if( epetraOpsTransp[0] == Teuchos::NO_TRANS && epetraOpsApplyMode[0] == PO::APPLY_MODE_APPLY )
00653     aztec_epetra_epetraFwdOp = epetra_epetraFwdOp;
00654   else
00655     aztec_epetra_epetraFwdOp = rcp(new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
00656   if( startingOver || aztec_epetra_epetraFwdOp.get() != aztecFwdSolver->GetUserOperator() ) {
00657     // Here we will be careful not to reset the forward operator in fears that
00658     // it will blow out the internally created stuff.
00659     aztecFwdSolver->SetUserOperator(const_cast<Epetra_Operator*>(&*aztec_epetra_epetraFwdOp));
00660     set_extra_data(aztec_epetra_epetraFwdOp,AOOLOWSF_aztec_epetra_epetraFwdOp_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00661   }
00662   // Adjoint solve
00663   if( aztecAdjSolver.get() ) {
00664     epetraOpsTransp[0] = ( epetra_epetraFwdOpTransp==NOTRANS ? Teuchos::TRANS : Teuchos::NO_TRANS );
00665     if( epetraOpsTransp[0] == Teuchos::NO_TRANS && epetraOpsApplyMode[0] == PO::APPLY_MODE_APPLY )
00666       aztec_epetra_epetraAdjOp = epetra_epetraFwdOp;
00667     else
00668       aztec_epetra_epetraAdjOp = rcp(new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
00669     aztecAdjSolver->SetUserOperator(const_cast<Epetra_Operator*>(&*aztec_epetra_epetraAdjOp));
00670     set_extra_data(aztec_epetra_epetraAdjOp,AOOLOWSF_aztec_epetra_epetraAdjOp_str,&aztecAdjSolver,Teuchos::POST_DESTROY,false);
00671   }
00672   //
00673   // Process the preconditioner
00674   //
00675   RefCountPtr<const Epetra_Operator>
00676     aztec_fwd_epetra_epetraPrecOp,
00677     aztec_adj_epetra_epetraPrecOp;
00678   bool setAztecPreconditioner = false;
00679   switch(localPrecType) {
00680     case PT_NONE: {
00681       //
00682       // No preconditioning at all!
00683       //
00684       break;
00685     }
00686     case PT_AZTEC_FROM_OP: {
00687       //
00688       // We are using the forward matrix for the preconditioner using aztec preconditioners
00689       //
00690       if( startingOver || !reusePrec ) {
00691         TEST_FOR_EXCEPTION(
00692           rowmatrix_epetraFwdOp.get()==NULL, std::logic_error
00693           ,"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): Error, There is no preconditioner given by client, but the client "
00694           "passed in an Epetra_Operator for the forward operator of type \'" <<typeid(*epetra_epetraFwdOp).name()<<"\' that does not "
00695           "support the Epetra_RowMatrix interface!"
00696           );
00697         TEST_FOR_EXCEPTION(
00698           epetra_epetraFwdOpTransp!=NOTRANS, std::logic_error
00699           ,"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): Error, There is no preconditioner given by client and the client "
00700           "passed in an Epetra_RowMatrix for the forward operator but the overall transpose is not NOTRANS and therefore we can can just "
00701           "hand this over to aztec without making a copy which is not supported here!"
00702           );
00703         aztecFwdSolver->SetPrecMatrix(const_cast<Epetra_RowMatrix*>(&*rowmatrix_epetraFwdOp));
00704         set_extra_data(rowmatrix_epetraFwdOp,AOOLOWSF_rowmatrix_epetraFwdOp_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00705       }
00706       setAztecPreconditioner = true;
00707       break;
00708     }
00709     case PT_AZTEC_FROM_APPROX_FWD_MATRIX: {
00710       //
00711       // The preconditioner comes from the input as a matrix and we are using aztec preconditioners
00712       //
00713       if( startingOver || !reusePrec ) {
00714         TEST_FOR_EXCEPTION(
00715           rowmatrix_epetraPrecOp.get()==NULL, std::logic_error
00716           ,"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): The client "
00717           "passed in an Epetra_Operator for the preconditioner matrix of type \'" <<typeid(*epetra_epetraPrecOp).name()<<"\' that does not "
00718           "support the Epetra_RowMatrix interface!"
00719           );
00720         TEST_FOR_EXCEPTION(
00721           overall_epetra_epetraPrecOpTransp!=NOTRANS, std::logic_error
00722           ,"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): Error, The client "
00723           "passed in an Epetra_RowMatrix for the preconditoner matrix but the overall transpose is not NOTRANS and therefore we can can just "
00724           "hand this over to aztec without making a copy which is not supported here!"
00725           );
00726         aztecFwdSolver->SetPrecMatrix(const_cast<Epetra_RowMatrix*>(&*rowmatrix_epetraPrecOp));
00727         set_extra_data(rowmatrix_epetraPrecOp,AOOLOWSF_rowmatrix_epetraPrecOp_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00728       }
00729       setAztecPreconditioner = true;
00730       break;
00731     }
00732     case PT_FROM_PREC_OP: {
00733       //
00734       // The preconditioner comes as an operator so let's use it as such
00735       //
00736       // Forawrd solve
00737       RefCountPtr<const Epetra_Operator>
00738         epetraOps[]
00739         = { epetra_epetraPrecOp };
00740       Teuchos::ETransp
00741         epetraOpsTransp[]
00742         = { overall_epetra_epetraPrecOpTransp==NOTRANS ? Teuchos::NO_TRANS : Teuchos::TRANS };
00743       PO::EApplyMode
00744         epetraOpsApplyMode[] // Here we must toggle the apply mode since aztecoo applies the preconditioner using ApplyInverse(...)
00745         = { epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY ? PO::APPLY_MODE_APPLY_INVERSE : PO::APPLY_MODE_APPLY };
00746       if( epetraOpsTransp[0] == Teuchos::NO_TRANS && epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY_INVERSE )
00747         aztec_fwd_epetra_epetraPrecOp = epetra_epetraPrecOp;
00748       else
00749         aztec_fwd_epetra_epetraPrecOp = rcp(new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
00750       aztecFwdSolver->SetPrecOperator(const_cast<Epetra_Operator*>(&*aztec_fwd_epetra_epetraPrecOp));
00751       set_extra_data(aztec_fwd_epetra_epetraPrecOp,AOOLOWSF_aztec_fwd_epetra_epetraPrecOp_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00752       // Adjoint solve
00753       if( aztecAdjSolver.get() && epetra_epetraPrecOpAdjointSupport == EPETRA_OP_ADJOINT_SUPPORTED ) {
00754         epetraOpsTransp[0] = ( overall_epetra_epetraPrecOpTransp==NOTRANS ? Teuchos::TRANS : Teuchos::NO_TRANS );
00755         if( epetraOpsTransp[0] == Teuchos::NO_TRANS && epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY_INVERSE )
00756           aztec_adj_epetra_epetraPrecOp = epetra_epetraPrecOp;
00757         else
00758           aztec_adj_epetra_epetraPrecOp = rcp(new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
00759         aztecAdjSolver->SetPrecOperator(const_cast<Epetra_Operator*>(&*aztec_adj_epetra_epetraPrecOp));
00760         set_extra_data(aztec_adj_epetra_epetraPrecOp,AOOLOWSF_aztec_adj_epetra_epetraPrecOp_str,&aztecAdjSolver,Teuchos::POST_DESTROY,false);
00761         set_extra_data<bool>(true,AOOLOWSF_setPrecondtionerOperator_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00762       }
00763       break;
00764     }
00765     default:
00766       TEST_FOR_EXCEPT(true);
00767   }
00768   //
00769   // Initialize the aztec preconditioner
00770   //
00771   if(setAztecPreconditioner) {
00772     if( startingOver || !reusePrec ) {
00773       double condNumEst = -1.0;
00774       TEST_FOR_EXCEPT(0!=aztecFwdSolver->ConstructPreconditioner(condNumEst));
00775       //aztecFwdSolver->SetAztecOption(AZ_pre_calc, AZ_calc);
00776       set_extra_data<bool>(true,AOOLOWSF_constructedAztecPreconditoner_str,&aztecFwdSolver,Teuchos::POST_DESTROY,false);
00777     }
00778     else {
00779       //aztecFwdSolver->SetAztecOption(AZ_pre_calc, AZ_reuse);
00780     }
00781   }
00782   //
00783   // Initialize the AztecOOLinearOpWithSolve object and set its options
00784   //
00785   if(aztecAdjSolver.get()) {
00786     aztecOp->initialize(
00787       fwdOp,fwdOpSrc,precUsed,prec.get()!=NULL,approxFwdOpSrc
00788       ,aztecFwdSolver,true,aztecAdjSolver,true,epetra_epetraFwdOpScalar
00789       );
00790   }
00791   else {
00792     aztecOp->initialize(
00793       fwdOp,fwdOpSrc,precUsed,prec.get()!=NULL,approxFwdOpSrc
00794       ,aztecFwdSolver,true,null,false,epetra_epetraFwdOpScalar
00795       );
00796   }
00797   aztecOp->fwdDefaultMaxIterations(defaultFwdMaxIterations_);
00798   aztecOp->fwdDefaultTol(defaultFwdTolerance_);
00799   aztecOp->adjDefaultMaxIterations(defaultAdjMaxIterations_);
00800   aztecOp->adjDefaultTol(defaultAdjTolerance_);
00801   aztecOp->outputEveryRhs(outputEveryRhs_);
00802   aztecOp->setOStream(this->getOStream());
00803   aztecOp->setVerbLevel(this->getVerbLevel());
00804 #ifdef TEUCHOS_DEBUG
00805   if(paramList_.get())
00806     paramList_->validateParameters(*this->getValidParameters());
00807 #endif
00808   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00809     *out << "\nLeaving Thyra::AztecOOLinearOpWithSolveFactory::initializeOp_impl(...) ...\n";
00810 }
00811 
00812 } // namespace Thyra
00813 
00814 #endif // __sun

Generated on Thu Sep 18 12:29:57 2008 for Aztecoo/Thyra Linear Solver Adapter Software by doxygen 1.3.9.1