Thyra_AmesosLinearOpWithSolve.cpp

00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //                Amesos: Direct Sparse Solver Package
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 */
00030 
00031 #ifndef SUN_CXX
00032 
00033 #include "Thyra_AmesosLinearOpWithSolve.hpp"
00034 #include "Thyra_EpetraThyraWrappers.hpp"
00035 #include "Epetra_MultiVector.h"
00036 #include "Teuchos_TimeMonitor.hpp"
00037 
00038 
00039 namespace Thyra {
00040 
00041 
00042 // Constructors/initializers/accessors
00043 
00044 
00045 AmesosLinearOpWithSolve::AmesosLinearOpWithSolve()
00046 {}
00047 
00048 
00049 AmesosLinearOpWithSolve::AmesosLinearOpWithSolve(
00050   const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
00051   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00052   const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
00053   const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
00054   const EOpTransp amesosSolverTransp,
00055   const double amesosSolverScalar
00056   )
00057 {
00058   this->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,
00059     amesosSolverTransp,amesosSolverScalar);
00060 }
00061 
00062 
00063 void AmesosLinearOpWithSolve::initialize(
00064   const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
00065   const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
00066   const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
00067   const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
00068   const EOpTransp amesosSolverTransp,
00069   const double amesosSolverScalar
00070   )
00071 {
00072 #ifdef TEUCHOS_DEBUG
00073   TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00074   TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00075   TEST_FOR_EXCEPT(epetraLP.get()==NULL);
00076   TEST_FOR_EXCEPT(amesosSolver.get()==NULL);
00077   TEST_FOR_EXCEPT(epetraLP->GetLHS()!=NULL);
00078   TEST_FOR_EXCEPT(epetraLP->GetRHS()!=NULL);
00079 #endif
00080   fwdOp_ = fwdOp;
00081   fwdOpSrc_ = fwdOpSrc;
00082   epetraLP_ = epetraLP;
00083   amesosSolver_ = amesosSolver;
00084   amesosSolverTransp_ = amesosSolverTransp;
00085   amesosSolverScalar_ = amesosSolverScalar;
00086   const std::string fwdOpLabel = fwdOp_->getObjectLabel();
00087   if(fwdOpLabel.length())
00088     this->setObjectLabel( "lows("+fwdOpLabel+")" );
00089 }
00090 
00091 
00092 Teuchos::RCP<const LinearOpSourceBase<double> >
00093 AmesosLinearOpWithSolve::extract_fwdOpSrc()
00094 {
00095   Teuchos::RCP<const LinearOpSourceBase<double> >
00096     _fwdOpSrc = fwdOpSrc_;
00097   fwdOpSrc_ = Teuchos::null;
00098   return _fwdOpSrc;
00099 }
00100 
00101 
00102 void AmesosLinearOpWithSolve::uninitialize(
00103   Teuchos::RCP<const LinearOpBase<double> > *fwdOp,
00104   Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
00105   Teuchos::RCP<Epetra_LinearProblem> *epetraLP,
00106   Teuchos::RCP<Amesos_BaseSolver> *amesosSolver,
00107   EOpTransp *amesosSolverTransp,
00108   double *amesosSolverScalar
00109   )
00110 {
00111 
00112   if(fwdOp) *fwdOp = fwdOp_;
00113   if(fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00114   if(epetraLP) *epetraLP = epetraLP_;
00115   if(amesosSolver) *amesosSolver = amesosSolver_;
00116   if(amesosSolverTransp) *amesosSolverTransp = amesosSolverTransp_;
00117   if(amesosSolverScalar) *amesosSolverScalar = amesosSolverScalar_;
00118 
00119   fwdOp_ = Teuchos::null;
00120   fwdOpSrc_ = Teuchos::null;
00121   epetraLP_ = Teuchos::null;
00122   amesosSolver_ = Teuchos::null;
00123   amesosSolverTransp_ = NOTRANS;
00124   amesosSolverScalar_ = 0.0;
00125 
00126 }
00127 
00128 
00129 // Overridden from LinearOpBase
00130 
00131 
00132 Teuchos::RCP< const VectorSpaceBase<double> >
00133 AmesosLinearOpWithSolve::range() const
00134 {
00135   return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
00136 }
00137 
00138 
00139 Teuchos::RCP< const VectorSpaceBase<double> >
00140 AmesosLinearOpWithSolve::domain() const
00141 {
00142   return  ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
00143 }
00144 
00145 
00146 Teuchos::RCP<const LinearOpBase<double> >
00147 AmesosLinearOpWithSolve::clone() const
00148 {
00149   return Teuchos::null; // Not supported yet but could be
00150 }
00151 
00152 
00153 // Overridden from Teuchos::Describable
00154 
00155 
00156 std::string AmesosLinearOpWithSolve::description() const
00157 {
00158   std::ostringstream oss;
00159   oss << Teuchos::Describable::description();
00160   if(!is_null(amesosSolver_)) {
00161     oss
00162       << "{fwdOp="<<fwdOp_->description()
00163       << ",amesosSolver="<<typeName(*amesosSolver_)<<"}";
00164   }
00165   return oss.str();
00166 }
00167 
00168 
00169 void AmesosLinearOpWithSolve::describe(
00170   Teuchos::FancyOStream &out,
00171   const Teuchos::EVerbosityLevel verbLevel
00172   ) const
00173 {
00174   using Teuchos::OSTab;
00175   using Teuchos::typeName;
00176   using Teuchos::describe;
00177   switch(verbLevel) {
00178     case Teuchos::VERB_DEFAULT:
00179     case Teuchos::VERB_LOW:
00180       out << this->description() << std::endl;
00181       break;
00182     case Teuchos::VERB_MEDIUM:
00183     case Teuchos::VERB_HIGH:
00184     case Teuchos::VERB_EXTREME:
00185     {
00186       out
00187         << Teuchos::Describable::description() << "{"
00188         << "rangeDim=" << this->range()->dim()
00189         << ",domainDim="<< this->domain()->dim() << "}\n";
00190       OSTab tab(out);
00191       if(!is_null(fwdOp_)) {
00192         out << "fwdOp = " << describe(*fwdOp_,verbLevel);
00193       }
00194       if(!is_null(amesosSolver_)) {
00195         out << "amesosSolver=" << typeName(*amesosSolver_) << "\n";
00196       }
00197       break;
00198     }
00199     default:
00200       TEST_FOR_EXCEPT(true); // Should never get here!
00201   }
00202 }
00203 
00204 
00205 // protected
00206 
00207 
00208 // Overridden from SingleScalarLinearOpBase
00209 
00210 
00211 bool AmesosLinearOpWithSolve::opSupported(EOpTransp M_trans) const
00212 {
00213   return ::Thyra::opSupported(*fwdOp_,M_trans);
00214 }
00215 
00216 
00217 void AmesosLinearOpWithSolve::apply(
00218   const EOpTransp M_trans,
00219   const MultiVectorBase<double> &X,
00220   MultiVectorBase<double> *Y,
00221   const double alpha,
00222   const double beta
00223   ) const
00224 {
00225   Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
00226 }
00227 
00228 
00229 // Overridden from SingleScalarLinearOpWithSolveBase
00230 
00231 
00232 bool AmesosLinearOpWithSolve::solveSupportsTrans(EOpTransp M_trans) const
00233 {
00234   return true; // ToDo: Determine if the solver supports adjoints or not!
00235 }
00236 
00237 
00238 bool AmesosLinearOpWithSolve::solveSupportsSolveMeasureType(
00239   EOpTransp M_trans, const SolveMeasureType& solveMeasureType
00240   ) const
00241 {
00242   return true; // I am a direct solver so I should be able to do it all!
00243 }
00244 
00245 
00246 // Overridden from SingleRhsLinearOpWithSolveBase
00247 
00248 
00249 void AmesosLinearOpWithSolve::solve(
00250   const EOpTransp M_trans,
00251   const MultiVectorBase<double> &B,
00252   MultiVectorBase<double> *X,
00253   const int numBlocks,
00254   const BlockSolveCriteria<double> blockSolveCriteria[],
00255   SolveStatus<double> blockSolveStatus[]
00256   ) const
00257 {
00258   using Teuchos::OSTab;
00259   typedef SolveCriteria<double> SC;
00260   typedef SolveStatus<double> SS;
00261 #ifdef TEUCHOS_DEBUG
00262   TEST_FOR_EXCEPT(X==NULL);
00263   TEST_FOR_EXCEPT(blockSolveCriteria==NULL && blockSolveStatus!=NULL);
00264 #endif
00265   Teuchos::Time totalTimer("");
00266   totalTimer.start(true);
00267   TEUCHOS_FUNC_TIME_MONITOR("AmesosLOWS");
00268   //
00269   Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
00270   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00271   OSTab tab = this->getOSTab();
00272   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00273     *out << "\nSolving block system using Amesos solver "
00274          << typeName(*amesosSolver_) << " ...\n\n";
00275   //
00276   // Get the op(...) range and domain maps
00277   //
00278   const EOpTransp amesosOpTransp = real_trans(trans_trans(amesosSolverTransp_,M_trans));
00279   const Epetra_Operator *amesosOp = epetraLP_->GetOperator();
00280   const Epetra_Map
00281     &opRangeMap  = ( amesosOpTransp == NOTRANS
00282       ? amesosOp->OperatorRangeMap()  : amesosOp->OperatorDomainMap() ),
00283     &opDomainMap = ( amesosOpTransp == NOTRANS
00284       ? amesosOp->OperatorDomainMap() : amesosOp->OperatorRangeMap()  );
00285   //
00286   // Get Epetra_MultiVector views of B and X
00287   //
00288   Teuchos::RCP<const Epetra_MultiVector>
00289     epetra_B = get_Epetra_MultiVector(opRangeMap,Teuchos::rcp(&B,false));
00290   Teuchos::RCP<Epetra_MultiVector>
00291     epetra_X = get_Epetra_MultiVector(opDomainMap,Teuchos::rcp(X,false));
00292   //
00293   // Set B and X in the linear problem
00294   //
00295   epetraLP_->SetLHS(&*epetra_X);
00296   epetraLP_->SetRHS(const_cast<Epetra_MultiVector*>(&*epetra_B));
00297   // Above should be okay but cross your fingers!
00298   //
00299   // Solve the linear system
00300   //
00301   const bool oldUseTranspose = amesosSolver_->UseTranspose();
00302   amesosSolver_->SetUseTranspose(amesosOpTransp==TRANS);
00303   const int err = amesosSolver_->Solve();
00304   TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
00305     "Error, the function Solve() on the amesos solver of type\n"
00306     "\'"<<typeName(*amesosSolver_)<<"\' failed with error code "<<err<<"!"
00307     );
00308   amesosSolver_->SetUseTranspose(oldUseTranspose);
00309   //
00310   // Unset B and X
00311   //
00312   epetraLP_->SetLHS(NULL);
00313   epetraLP_->SetRHS(NULL);
00314   epetra_X = Teuchos::null;
00315   epetra_B = Teuchos::null;
00316   //
00317   // Scale X if needed
00318   //
00319   if(amesosSolverScalar_!=1.0)
00320     Thyra::scale(1.0/amesosSolverScalar_,X);
00321   //
00322   // Set the solve status if requested
00323   //
00324   if(numBlocks && blockSolveStatus) {
00325     for( int i = 0; i < numBlocks; ++i ) {
00326       blockSolveStatus[i].solveStatus = SOLVE_STATUS_CONVERGED;
00327       blockSolveStatus[i].achievedTol = SS::unknownTolerance();
00328       blockSolveStatus[i].message
00329         = std::string("Solver ")+typeName(*amesosSolver_)+std::string(" converged!");
00330     }
00331   }
00332   //
00333   // Report the overall time
00334   //
00335   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00336     *out
00337       << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00338 }
00339 
00340 
00341 } // end namespace Thyra
00342 
00343 #endif // SUN_CXX

Generated on Wed Jul 22 12:56:30 2009 for Stratimikos by  doxygen 1.5.8