Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_AmesosLinearOpWithSolve.cpp
Go to the documentation of this file.
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Stratimikos: Thyra-based strategies for linear solvers
6 // Copyright (2006) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
46 #include "Thyra_MultiVectorStdOps.hpp"
47 #include "Epetra_MultiVector.h"
48 #include "Teuchos_TimeMonitor.hpp"
49 
50 
51 namespace Thyra {
52 
53 
54 // Constructors/initializers/accessors
55 
56 
58  amesosSolverTransp_(Thyra::NOTRANS),
59  amesosSolverScalar_(1.0)
60 {}
61 
62 
64  const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
65  const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
66  const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
67  const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
68  const EOpTransp amesosSolverTransp,
69  const double amesosSolverScalar
70  )
71 {
72  this->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,
73  amesosSolverTransp,amesosSolverScalar);
74 }
75 
76 
78  const Teuchos::RCP<const LinearOpBase<double> > &fwdOp,
79  const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc,
80  const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
81  const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
82  const EOpTransp amesosSolverTransp,
83  const double amesosSolverScalar
84  )
85 {
86 #ifdef TEUCHOS_DEBUG
87  TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
88  TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
89  TEUCHOS_TEST_FOR_EXCEPT(epetraLP.get()==NULL);
90  TEUCHOS_TEST_FOR_EXCEPT(amesosSolver.get()==NULL);
91  TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetLHS()!=NULL);
92  TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetRHS()!=NULL);
93 #endif
94  fwdOp_ = fwdOp;
95  fwdOpSrc_ = fwdOpSrc;
96  epetraLP_ = epetraLP;
97  amesosSolver_ = amesosSolver;
98  amesosSolverTransp_ = amesosSolverTransp;
99  amesosSolverScalar_ = amesosSolverScalar;
100  const std::string fwdOpLabel = fwdOp_->getObjectLabel();
101  if(fwdOpLabel.length())
102  this->setObjectLabel( "lows("+fwdOpLabel+")" );
103 }
104 
105 
108 {
110  _fwdOpSrc = fwdOpSrc_;
112  return _fwdOpSrc;
113 }
114 
115 
117  Teuchos::RCP<const LinearOpBase<double> > *fwdOp,
118  Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
120  Teuchos::RCP<Amesos_BaseSolver> *amesosSolver,
121  EOpTransp *amesosSolverTransp,
122  double *amesosSolverScalar
123  )
124 {
125 
126  if(fwdOp) *fwdOp = fwdOp_;
127  if(fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
128  if(epetraLP) *epetraLP = epetraLP_;
129  if(amesosSolver) *amesosSolver = amesosSolver_;
130  if(amesosSolverTransp) *amesosSolverTransp = amesosSolverTransp_;
131  if(amesosSolverScalar) *amesosSolverScalar = amesosSolverScalar_;
132 
138  amesosSolverScalar_ = 0.0;
139 
140 }
141 
142 
143 // Overridden from LinearOpBase
144 
145 
148 {
149  return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
150 }
151 
152 
155 {
156  return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
157 }
158 
159 
162 {
163  return Teuchos::null; // Not supported yet but could be
164 }
165 
166 
167 // Overridden from Teuchos::Describable
168 
169 
171 {
172  std::ostringstream oss;
174  if(!is_null(amesosSolver_)) {
175  oss
176  << "{fwdOp="<<fwdOp_->description()
177  << ",amesosSolver="<<typeName(*amesosSolver_)<<"}";
178  }
179  return oss.str();
180 }
181 
182 
185  const Teuchos::EVerbosityLevel verbLevel
186  ) const
187 {
188  using Teuchos::OSTab;
189  using Teuchos::typeName;
190  using Teuchos::describe;
191  switch(verbLevel) {
193  case Teuchos::VERB_LOW:
194  out << this->description() << std::endl;
195  break;
197  case Teuchos::VERB_HIGH:
199  {
200  out
202  << "rangeDim=" << this->range()->dim()
203  << ",domainDim="<< this->domain()->dim() << "}\n";
204  OSTab tab(out);
205  if(!is_null(fwdOp_)) {
206  out << "fwdOp = " << describe(*fwdOp_,verbLevel);
207  }
208  if(!is_null(amesosSolver_)) {
209  out << "amesosSolver=" << typeName(*amesosSolver_) << "\n";
210  }
211  break;
212  }
213  default:
214  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
215  }
216 }
217 
218 
219 // protected
220 
221 
222 // Overridden from LinearOpBase
223 
224 
225 bool AmesosLinearOpWithSolve::opSupportedImpl(EOpTransp M_trans) const
226 {
227  return ::Thyra::opSupported(*fwdOp_,M_trans);
228 }
229 
230 
232  const EOpTransp M_trans,
233  const MultiVectorBase<double> &X,
234  const Ptr<MultiVectorBase<double> > &Y,
235  const double alpha,
236  const double beta
237  ) const
238 {
239  Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
240 }
241 
242 
243 // Overridden from LinearOpWithSolveBase
244 
245 
246 bool AmesosLinearOpWithSolve::solveSupportsImpl(EOpTransp M_trans) const
247 {
248  if (Thyra::real_trans(M_trans) == Thyra::NOTRANS) {
249  // Assume every amesos solver supports a basic forward solve!
250  return true;
251  }
252  // Query the amesos solver to see if it supports the transpose operation.
253  // NOTE: Amesos_BaseSolver makes you change the state of the object to
254  // determine if the object supports an adjoint solver. This is a bad design
255  // but I have no control over that. This is why you see this hacked
256  // oldUseTranspose variable and logic. NOTE: This function meets the basic
257  // guarantee but if setUseTransplse(...) throws, then the state of
258  // UseTranspose() may be different.
259  const bool oldUseTranspose = amesosSolver_->UseTranspose();
260  const bool supportsAdjoint = (amesosSolver_->SetUseTranspose(true) == 0);
261  amesosSolver_->SetUseTranspose(oldUseTranspose);
262  return supportsAdjoint;
263 }
264 
265 
267  EOpTransp /* M_trans */, const SolveMeasureType& /* solveMeasureType */
268  ) const
269 {
270  return true; // I am a direct solver so I should be able to do it all!
271 }
272 
273 
274 SolveStatus<double>
276  const EOpTransp M_trans,
277  const MultiVectorBase<double> &B,
278  const Ptr<MultiVectorBase<double> > &X,
279  const Ptr<const SolveCriteria<double> > /* solveCriteria */
280  ) const
281 {
282  using Teuchos::rcpFromPtr;
283  using Teuchos::rcpFromRef;
284  using Teuchos::OSTab;
285 
286  Teuchos::Time totalTimer("");
287  totalTimer.start(true);
288 
289  THYRA_FUNC_TIME_MONITOR("Stratimikos: AmesosLOWS");
290 
291  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
292  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
293  OSTab tab = this->getOSTab();
294  if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
295  *out << "\nSolving block system using Amesos solver "
296  << typeName(*amesosSolver_) << " ...\n\n";
297 
298  //
299  // Get the op(...) range and domain maps
300  //
301  const EOpTransp amesosOpTransp = real_trans(trans_trans(amesosSolverTransp_,M_trans));
302  const Epetra_Operator *amesosOp = epetraLP_->GetOperator();
303  const Epetra_Map
304  &opRangeMap = ( amesosOpTransp == NOTRANS
305  ? amesosOp->OperatorRangeMap() : amesosOp->OperatorDomainMap() ),
306  &opDomainMap = ( amesosOpTransp == NOTRANS
307  ? amesosOp->OperatorDomainMap() : amesosOp->OperatorRangeMap() );
308 
309  //
310  // Get Epetra_MultiVector views of B and X
311  //
313  epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
315  epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
316 
317  //
318  // Set B and X in the linear problem
319  //
320  epetraLP_->SetLHS(&*epetra_X);
321  epetraLP_->SetRHS(const_cast<Epetra_MultiVector*>(&*epetra_B));
322  // Above should be okay but cross your fingers!
323 
324  //
325  // Solve the linear system
326  //
327  const bool oldUseTranspose = amesosSolver_->UseTranspose();
328  amesosSolver_->SetUseTranspose(amesosOpTransp==TRANS);
329  const int err = amesosSolver_->Solve();
330  TEUCHOS_TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
331  "Error, the function Solve() on the amesos solver of type\n"
332  "\'"<<typeName(*amesosSolver_)<<"\' failed with error code "<<err<<"!"
333  );
334  amesosSolver_->SetUseTranspose(oldUseTranspose);
335 
336  //
337  // Unset B and X
338  //
339  epetraLP_->SetLHS(NULL);
340  epetraLP_->SetRHS(NULL);
341  epetra_X = Teuchos::null;
342  epetra_B = Teuchos::null;
343 
344  //
345  // Scale X if needed
346  //
347  if(amesosSolverScalar_!=1.0)
348  Thyra::scale(1.0/amesosSolverScalar_, X);
349 
350  //
351  // Set the solve status if requested
352  //
353  SolveStatus<double> solveStatus;
354  solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
355  solveStatus.achievedTol = SolveStatus<double>::unknownTolerance();
356  solveStatus.message =
357  std::string("Solver ")+typeName(*amesosSolver_)+std::string(" converged!");
358 
359  //
360  // Report the overall time
361  //
362  if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
363  *out
364  << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
365 
366  return solveStatus;
367 
368 }
369 
370 
371 } // end namespace Thyra
Teuchos::RCP< Epetra_LinearProblem > epetraLP_
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Teuchos::RCP< const LinearOpSourceBase< double > > extract_fwdOpSrc()
Extract the LinearOpSourceBase&lt;double&gt; object so that it can be modified.
std::string typeName(const T &t)
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::RCP< Amesos_BaseSolver > amesosSolver_
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual bool opSupportedImpl(EOpTransp M_trans) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T * get() const
basic_OSTab< char > OSTab
virtual bool solveSupportsImpl(EOpTransp M_trans) const
void start(bool reset=false)
SolveStatus< double > solveImpl(const EOpTransp M_trans, const MultiVectorBase< double > &B, const Ptr< MultiVectorBase< double > > &X, const Ptr< const SolveCriteria< double > > solveCriteria) const
Teuchos::RCP< const VectorSpaceBase< double > > range() const
Teuchos::RCP< const LinearOpSourceBase< double > > fwdOpSrc_
virtual std::string description() const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void initialize(const Teuchos::RCP< const LinearOpBase< double > > &fwdOp, const Teuchos::RCP< const LinearOpSourceBase< double > > &fwdOpSrc, const Teuchos::RCP< Epetra_LinearProblem > &epetraLP, const Teuchos::RCP< Amesos_BaseSolver > &amesosSolver, const EOpTransp amesosSolverTransp, const double amesosSolverScalar)
First initialization.
void uninitialize(Teuchos::RCP< const LinearOpBase< double > > *fwdOp=NULL, Teuchos::RCP< const LinearOpSourceBase< double > > *fwdOpSrc=NULL, Teuchos::RCP< Epetra_LinearProblem > *epetraLP=NULL, Teuchos::RCP< Amesos_BaseSolver > *amesosSolver=NULL, EOpTransp *amesosSolverTransp=NULL, double *amesosSolverScalar=NULL)
Uninitialize.
virtual void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
Teuchos::RCP< const LinearOpBase< double > > fwdOp_
Teuchos::RCP< const LinearOpBase< double > > clone() const
double totalElapsedTime(bool readCurrentTime=false) const
Teuchos::RCP< const VectorSpaceBase< double > > domain() const
AmesosLinearOpWithSolve()
Construct to uninitialized.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)