Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
45#include "Thyra_EpetraThyraWrappers.hpp"
46#include "Thyra_MultiVectorStdOps.hpp"
47#include "Epetra_MultiVector.h"
48#include "Teuchos_TimeMonitor.hpp"
49
50
51namespace 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
106Teuchos::RCP<const LinearOpSourceBase<double> >
108{
109 Teuchos::RCP<const LinearOpSourceBase<double> >
110 _fwdOpSrc = fwdOpSrc_;
111 fwdOpSrc_ = Teuchos::null;
112 return _fwdOpSrc;
113}
114
115
117 Teuchos::RCP<const LinearOpBase<double> > *fwdOp,
118 Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
119 Teuchos::RCP<Epetra_LinearProblem> *epetraLP,
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
133 fwdOp_ = Teuchos::null;
134 fwdOpSrc_ = Teuchos::null;
135 epetraLP_ = Teuchos::null;
136 amesosSolver_ = Teuchos::null;
137 amesosSolverTransp_ = NOTRANS;
139
140}
141
142
143// Overridden from LinearOpBase
144
145
146Teuchos::RCP< const VectorSpaceBase<double> >
148{
149 return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
150}
151
152
153Teuchos::RCP< const VectorSpaceBase<double> >
155{
156 return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
157}
158
159
160Teuchos::RCP<const LinearOpBase<double> >
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;
173 oss << Teuchos::Describable::description();
174 if(!is_null(amesosSolver_)) {
175 oss
176 << "{fwdOp="<<fwdOp_->description()
177 << ",amesosSolver="<<typeName(*amesosSolver_)<<"}";
178 }
179 return oss.str();
180}
181
182
184 Teuchos::FancyOStream &out,
185 const Teuchos::EVerbosityLevel verbLevel
186 ) const
187{
188 using Teuchos::OSTab;
189 using Teuchos::typeName;
190 using Teuchos::describe;
191 switch(verbLevel) {
192 case Teuchos::VERB_DEFAULT:
193 case Teuchos::VERB_LOW:
194 out << this->description() << std::endl;
195 break;
196 case Teuchos::VERB_MEDIUM:
197 case Teuchos::VERB_HIGH:
198 case Teuchos::VERB_EXTREME:
199 {
200 out
201 << Teuchos::Describable::description() << "{"
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
225bool 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
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
274SolveStatus<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 //
312 Teuchos::RCP<const Epetra_MultiVector>
313 epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
314 Teuchos::RCP<Epetra_MultiVector>
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
SolveStatus< double > solveImpl(const EOpTransp M_trans, const MultiVectorBase< double > &B, const Ptr< MultiVectorBase< double > > &X, const Ptr< const SolveCriteria< double > > solveCriteria) const
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 LinearOpSourceBase< double > > fwdOpSrc_
Teuchos::RCP< const LinearOpBase< double > > clone() const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
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.
Teuchos::RCP< const VectorSpaceBase< double > > domain() 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.
virtual bool opSupportedImpl(EOpTransp M_trans) const
virtual bool solveSupportsImpl(EOpTransp M_trans) const
Teuchos::RCP< const LinearOpBase< double > > fwdOp_
Teuchos::RCP< const VectorSpaceBase< double > > range() const
Teuchos::RCP< const LinearOpSourceBase< double > > extract_fwdOpSrc()
Extract the LinearOpSourceBase<double> object so that it can be modified.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Teuchos::RCP< Epetra_LinearProblem > epetraLP_
Teuchos::RCP< Amesos_BaseSolver > amesosSolver_
AmesosLinearOpWithSolve()
Construct to uninitialized.