Stratimikos Version of the Day
Loading...
Searching...
No Matches
Thyra_AztecOOLinearOpWithSolve.cpp
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Thyra_AztecOOLinearOpWithSolve.hpp"
43#include "Thyra_LinearOpWithSolveHelpers.hpp"
44#include "Thyra_EpetraThyraWrappers.hpp"
45#include "Thyra_EpetraOperatorWrapper.hpp"
46#include "Teuchos_BLAS_types.hpp"
47#include "Teuchos_TimeMonitor.hpp"
48#include "Teuchos_Time.hpp"
49#include "Teuchos_implicit_cast.hpp"
50
51
52namespace {
53
54#if 0
55// unused implementation detail internal to this file
56inline
57Teuchos::ETransp convert( Thyra::EOpTransp trans_in )
58{
59 Teuchos::ETransp trans_out;
60 switch(trans_in) {
61 case Thyra::NOTRANS:
62 trans_out = Teuchos::NO_TRANS;
63 break;
64 case Thyra::TRANS:
65 trans_out = Teuchos::TRANS;
66 break;
67 default:
68 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
69 }
70 return trans_out;
71}
72#endif // 0
73
74
75// This class sets some solve instance specific state and then sets it back to
76// the default state on destruction. But using the destructor to unset the
77// state we can be sure that the state is rest correctly even if an exception
78// is thrown.
79class SetAztecSolveState {
80public:
81 SetAztecSolveState(
82 const Teuchos::RCP<AztecOO> &aztecSolver,
83 const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
84 const Teuchos::EVerbosityLevel verbLevel,
85 const Thyra::SolveMeasureType &solveMeasureType
86 );
87 ~SetAztecSolveState();
88private:
89 Teuchos::RCP<AztecOO> aztecSolver_;
90 Teuchos::RCP<Teuchos::FancyOStream> fancyOStream_;
91 Teuchos::EVerbosityLevel verbLevel_;
92 int outputFrequency_;
93 int convergenceTest_;
94 SetAztecSolveState(); // Not defined and not to be called!
95};
96
97
98SetAztecSolveState::SetAztecSolveState(
99 const Teuchos::RCP<AztecOO> &aztecSolver,
100 const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
101 const Teuchos::EVerbosityLevel verbLevel,
102 const Thyra::SolveMeasureType &solveMeasureType
103 )
104 :aztecSolver_(aztecSolver.assert_not_null())
105 ,outputFrequency_(0)
106{
107
108 // Output state
109 verbLevel_ = verbLevel;
110 if ( Teuchos::VERB_NONE != verbLevel_ ) {
111 if (!is_null(fancyOStream)) {
112 // AztecOO puts in two tabs before it prints anything. Therefore,
113 // there is not much that we can do to improve the layout of the
114 // indentation so just leave it!
115 fancyOStream_= Teuchos::tab(
116 fancyOStream.create_weak(),
117 0, // Don't indent since AztecOO already puts in two tabs (not spaces!)
118 Teuchos::implicit_cast<std::string>("AZTECOO")
119 );
120 aztecSolver_->SetOutputStream(*fancyOStream_);
121 aztecSolver_->SetErrorStream(*fancyOStream_);
122 // Note, above we can not save the current output and error streams
123 // since AztecOO does not define functions to get them. In the
124 // future, AztecOO should define these functions if we are to avoid
125 // treading on each others print statements. However, since the
126 // AztecOO object is most likely owned by these Thyra wrappers, this
127 // should not be a problem.
128 }
129 }
130 else {
131 outputFrequency_ = aztecSolver_->GetAllAztecOptions()[AZ_output];
132 aztecSolver_->SetAztecOption(AZ_output,0);
133 }
134
135 // Convergence test
136 convergenceTest_ = aztecSolver_->GetAztecOption(AZ_conv);
137 if (solveMeasureType.useDefault())
138 {
139 // Just use the default solve measure type already set!
140 }
141 else if (
142 solveMeasureType(
143 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
144 Thyra::SOLVE_MEASURE_NORM_RHS
145 )
146 )
147 {
148 aztecSolver_->SetAztecOption(AZ_conv,AZ_rhs);
149 }
150 else if (
151 solveMeasureType(
152 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
153 Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL
154 )
155 )
156 {
157 aztecSolver_->SetAztecOption(AZ_conv,AZ_r0);
158 }
159 else {
160 TEUCHOS_TEST_FOR_EXCEPT("Invalid solve measure type, you should never get here!");
161 }
162
163}
164
165
166SetAztecSolveState::~SetAztecSolveState()
167{
168
169 // Output state
170 if ( Teuchos::VERB_NONE != verbLevel_ ) {
171 if (!is_null(fancyOStream_)) {
172 aztecSolver_->SetOutputStream(std::cout);
173 aztecSolver_->SetErrorStream(std::cerr);
174 *fancyOStream_ << "\n";
175 }
176 }
177 else {
178 aztecSolver_->SetAztecOption(AZ_output,outputFrequency_);
179 }
180
181 // Convergence test
182 aztecSolver_->SetAztecOption(AZ_conv,convergenceTest_);
183
184}
185
186
187} // namespace
188
189
190namespace Thyra {
191
192
193// Constructors/initializers/accessors
194
195
197 const int fwdDefaultMaxIterations_in
198 ,const double fwdDefaultTol_in
199 ,const int adjDefaultMaxIterations_in
200 ,const double adjDefaultTol_in
201 ,const bool outputEveryRhs_in
202 )
203 :fwdDefaultMaxIterations_(fwdDefaultMaxIterations_in)
204 ,fwdDefaultTol_(fwdDefaultTol_in)
205 ,adjDefaultMaxIterations_(adjDefaultMaxIterations_in)
206 ,adjDefaultTol_(adjDefaultTol_in)
207 ,outputEveryRhs_(outputEveryRhs_in)
208 ,isExternalPrec_(false)
209 ,allowInexactFwdSolve_(false)
210 ,allowInexactAdjSolve_(false)
211 ,aztecSolverScalar_(0.0)
212{}
213
214
216 const RCP<const LinearOpBase<double> > &fwdOp
217 ,const RCP<const LinearOpSourceBase<double> > &fwdOpSrc
218 ,const RCP<const PreconditionerBase<double> > &prec
219 ,const bool isExternalPrec_in
220 ,const RCP<const LinearOpSourceBase<double> > &approxFwdOpSrc
221 ,const RCP<AztecOO> &aztecFwdSolver
222 ,const bool allowInexactFwdSolve
223 ,const RCP<AztecOO> &aztecAdjSolver
224 ,const bool allowInexactAdjSolve
225 ,const double aztecSolverScalar
226 )
227{
228#ifdef TEUCHOS_DEBUG
229 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
230 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
231 TEUCHOS_TEST_FOR_EXCEPT(aztecFwdSolver.get()==NULL);
232#endif
233 fwdOp_ = fwdOp;
234 fwdOpSrc_ = fwdOpSrc;
235 isExternalPrec_ = isExternalPrec_in;
236 prec_ = prec;
237 approxFwdOpSrc_ = approxFwdOpSrc;
238 aztecFwdSolver_ = aztecFwdSolver;
239 allowInexactFwdSolve_ = allowInexactFwdSolve;
240 aztecAdjSolver_ = aztecAdjSolver;
241 allowInexactAdjSolve_ = allowInexactAdjSolve;
242 aztecSolverScalar_ = aztecSolverScalar;
243 const std::string fwdOpLabel = fwdOp_->getObjectLabel();
244 if (fwdOpLabel.length())
245 this->setObjectLabel( "lows("+fwdOpLabel+")" );
246}
247
248
249RCP<const LinearOpSourceBase<double> >
251{
252 RCP<const LinearOpSourceBase<double> >
253 _fwdOpSrc = fwdOpSrc_;
254 fwdOpSrc_ = Teuchos::null;
255 return _fwdOpSrc;
256}
257
258
259RCP<const PreconditionerBase<double> >
261{
262 RCP<const PreconditionerBase<double> >
263 _prec = prec_;
264 prec_ = Teuchos::null;
265 return _prec;
266}
267
268
270{
271 return isExternalPrec_;
272}
273
274
275RCP<const LinearOpSourceBase<double> >
277{
278 RCP<const LinearOpSourceBase<double> >
279 _approxFwdOpSrc = approxFwdOpSrc_;
280 approxFwdOpSrc_ = Teuchos::null;
281 return _approxFwdOpSrc;
282}
283
284
286 RCP<const LinearOpBase<double> > *fwdOp,
287 RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
288 RCP<const PreconditionerBase<double> > *prec,
289 bool *isExternalPrec_inout,
290 RCP<const LinearOpSourceBase<double> > *approxFwdOpSrc,
291 RCP<AztecOO> *aztecFwdSolver,
292 bool *allowInexactFwdSolve,
293 RCP<AztecOO> *aztecAdjSolver,
294 bool *allowInexactAdjSolve,
295 double *aztecSolverScalar
296 )
297{
298 if (fwdOp) *fwdOp = fwdOp_;
299 if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
300 if (prec) *prec = prec_;
301 if (isExternalPrec_inout) *isExternalPrec_inout = isExternalPrec_;
302 if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
303 if (aztecFwdSolver) *aztecFwdSolver = aztecFwdSolver_;
304 if (allowInexactFwdSolve) *allowInexactFwdSolve = allowInexactFwdSolve_;
305 if (aztecAdjSolver) *aztecAdjSolver = aztecAdjSolver_;
306 if (allowInexactAdjSolve) *allowInexactAdjSolve = allowInexactAdjSolve_;
307 if (aztecSolverScalar) *aztecSolverScalar = aztecSolverScalar_;
308
309 fwdOp_ = Teuchos::null;
310 fwdOpSrc_ = Teuchos::null;
311 prec_ = Teuchos::null;
312 isExternalPrec_ = false; // Just to make unique
313 approxFwdOpSrc_ = Teuchos::null;
314 aztecFwdSolver_ = Teuchos::null;
315 allowInexactFwdSolve_ = false;
316 aztecAdjSolver_ = Teuchos::null;
317 allowInexactAdjSolve_ = false;
318 aztecSolverScalar_ = 0.0;
319}
320
321
322// Overridden from LinearOpBase
323
324
325RCP< const VectorSpaceBase<double> >
327{
328 return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
329}
330
331
332RCP< const VectorSpaceBase<double> >
334{
335 return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
336}
337
338
339RCP<const LinearOpBase<double> >
341{
342 return Teuchos::null; // Not supported yet but could be
343}
344
345
346// Overridden from Teuchos::Describable
347
348
350{
351 std::ostringstream oss;
352 oss << Teuchos::Describable::description();
353 if (fwdOp_.get()) {
354 oss << "{";
355 oss << "fwdOp="<<fwdOp_->description()<<"";
356 oss << "}";
357 }
358 return oss.str();
359}
360
361
363 Teuchos::FancyOStream &out,
364 const Teuchos::EVerbosityLevel verbLevel
365 ) const
366{
367 using Teuchos::OSTab;
368 using Teuchos::typeName;
369 using Teuchos::describe;
370 switch(verbLevel) {
371 case Teuchos::VERB_DEFAULT:
372 case Teuchos::VERB_LOW:
373 out << this->description() << std::endl;
374 break;
375 case Teuchos::VERB_MEDIUM:
376 case Teuchos::VERB_HIGH:
377 case Teuchos::VERB_EXTREME:
378 {
379 out
380 << Teuchos::Describable::description() << "{"
381 << "rangeDim=" << this->range()->dim()
382 << ",domainDim="<< this->domain()->dim() << "}\n";
383 OSTab tab(out);
384 if (!is_null(fwdOp_)) {
385 out << "fwdOp = " << describe(*fwdOp_,verbLevel);
386 }
387 if (!is_null(prec_)) {
388 out << "prec = " << describe(*prec_,verbLevel);
389 }
390 if (!is_null(aztecFwdSolver_)) {
391 if (aztecFwdSolver_->GetUserOperator())
392 out
393 << "Aztec Fwd Op = "
394 << typeName(*aztecFwdSolver_->GetUserOperator()) << "\n";
395 if (aztecFwdSolver_->GetUserMatrix())
396 out
397 << "Aztec Fwd Mat = "
398 << typeName(*aztecFwdSolver_->GetUserMatrix()) << "\n";
399 if (aztecFwdSolver_->GetPrecOperator())
400 out
401 << "Aztec Fwd Prec Op = "
402 << typeName(*aztecFwdSolver_->GetPrecOperator()) << "\n";
403 if (aztecFwdSolver_->GetPrecMatrix())
404 out
405 << "Aztec Fwd Prec Mat = "
406 << typeName(*aztecFwdSolver_->GetPrecMatrix()) << "\n";
407 }
408 if (!is_null(aztecAdjSolver_)) {
409 if (aztecAdjSolver_->GetUserOperator())
410 out
411 << "Aztec Adj Op = "
412 << typeName(*aztecAdjSolver_->GetUserOperator()) << "\n";
413 if (aztecAdjSolver_->GetUserMatrix())
414 out
415 << "Aztec Adj Mat = "
416 << typeName(*aztecAdjSolver_->GetUserMatrix()) << "\n";
417 if (aztecAdjSolver_->GetPrecOperator())
418 out
419 << "Aztec Adj Prec Op = "
420 << typeName(*aztecAdjSolver_->GetPrecOperator()) << "\n";
421 if (aztecAdjSolver_->GetPrecMatrix())
422 out
423 << "Aztec Adj Prec Mat = "
424 << typeName(*aztecAdjSolver_->GetPrecMatrix()) << "\n";
425 }
426 break;
427 }
428 default:
429 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
430 }
431}
432
433
434// protected
435
436
437// Overridden from LinearOpBase
438
439
440bool AztecOOLinearOpWithSolve::opSupportedImpl(EOpTransp M_trans) const
441{
442 return ::Thyra::opSupported(*fwdOp_,M_trans);
443}
444
445
447 const EOpTransp M_trans,
448 const MultiVectorBase<double> &X,
449 const Ptr<MultiVectorBase<double> > &Y,
450 const double alpha,
451 const double beta
452 ) const
453{
454 Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
455}
456
457
458// Overridden from LinearOpWithSolveBase
459
460
461bool
463{
464 if (real_trans(M_trans)==NOTRANS) return true;
465 return (nonnull(aztecAdjSolver_));
466}
467
468
469bool
471 EOpTransp M_trans, const SolveMeasureType& solveMeasureType
472 ) const
473{
474 if (real_trans(M_trans)==NOTRANS) {
475 if (solveMeasureType.useDefault())
476 {
477 return true;
478 }
479 else if (
480 solveMeasureType(
481 SOLVE_MEASURE_NORM_RESIDUAL,
482 SOLVE_MEASURE_NORM_RHS
483 )
484 &&
485 allowInexactFwdSolve_
486 )
487 {
488 return true;
489 }
490 else if (
491 solveMeasureType(
492 SOLVE_MEASURE_NORM_RESIDUAL,
493 SOLVE_MEASURE_NORM_INIT_RESIDUAL
494 )
495 &&
496 allowInexactFwdSolve_
497 )
498 {
499 return true;
500 }
501 }
502 else {
503 // TRANS
504 if (aztecAdjSolver_.get()==NULL)
505 {
506 return false;
507 }
508 else if (solveMeasureType.useDefault())
509 {
510 return true;
511 }
512 else if (
513 solveMeasureType(
514 SOLVE_MEASURE_NORM_RESIDUAL,
515 SOLVE_MEASURE_NORM_RHS
516 )
517 &&
518 allowInexactFwdSolve_
519 )
520 {
521 return true;
522 }
523 else if (
524 solveMeasureType(
525 SOLVE_MEASURE_NORM_RESIDUAL,
526 SOLVE_MEASURE_NORM_INIT_RESIDUAL
527 )
528 &&
529 allowInexactFwdSolve_
530 )
531 {
532 return true;
533 }
534 }
535 // If you get here then we don't support the solve measure type!
536 return false;
537}
538
539
540// Overridden from LinearOpWithSolveBase
541
542
543SolveStatus<double>
545 const EOpTransp M_trans,
546 const MultiVectorBase<double> &B,
547 const Ptr<MultiVectorBase<double> > &X,
548 const Ptr<const SolveCriteria<double> > solveCriteria
549 ) const
550{
551
552 using Teuchos::rcp;
553 using Teuchos::rcpFromRef;
554 using Teuchos::rcpFromPtr;
555 using Teuchos::OSTab;
556 typedef SolveCriteria<double> SC;
557 typedef SolveStatus<double> SS;
558
559 THYRA_FUNC_TIME_MONITOR("Stratimikos: AztecOOLOWS");
560 Teuchos::Time totalTimer(""), timer("");
561 totalTimer.start(true);
562
563 RCP<Teuchos::FancyOStream> out = this->getOStream();
564 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
565 OSTab tab = this->getOSTab();
566 if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
567 *out << "\nSolving block system using AztecOO ...\n\n";
568
569 //
570 // Validate input
571 //
572 TEUCHOS_ASSERT(this->solveSupportsImpl(M_trans));
573 SolveMeasureType solveMeasureType;
574 if (nonnull(solveCriteria)) {
575 solveMeasureType = solveCriteria->solveMeasureType;
576 assertSupportsSolveMeasureType(*this, M_trans, solveMeasureType);
577 }
578
579 //
580 // Get the transpose argument
581 //
582 const EOpTransp aztecOpTransp = real_trans(M_trans);
583
584 //
585 // Get the solver, operator, and preconditioner that we will use
586 //
587 RCP<AztecOO>
588 aztecSolver = ( aztecOpTransp == NOTRANS ? aztecFwdSolver_ : aztecAdjSolver_ );
589 const Epetra_Operator
590 *aztecOp = aztecSolver->GetUserOperator();
591
592 //
593 // Get the op(...) range and domain maps
594 //
595 const Epetra_Map
596 &opRangeMap = aztecOp->OperatorRangeMap(),
597 &opDomainMap = aztecOp->OperatorDomainMap();
598
599 //
600 // Get the convergence criteria
601 //
602 double tol = ( aztecOpTransp==NOTRANS ? fwdDefaultTol() : adjDefaultTol() );
603 int maxIterations = ( aztecOpTransp==NOTRANS
604 ? fwdDefaultMaxIterations() : adjDefaultMaxIterations() );
605 bool isDefaultSolveCriteria = true;
606 if (nonnull(solveCriteria)) {
607 if ( solveCriteria->requestedTol != SC::unspecifiedTolerance() ) {
608 tol = solveCriteria->requestedTol;
609 isDefaultSolveCriteria = false;
610 }
611 if (nonnull(solveCriteria->extraParameters)) {
612 maxIterations = solveCriteria->extraParameters->get("Maximum Iterations",maxIterations);
613 }
614 }
615
616 //
617 // Get Epetra_MultiVector views of B and X
618 //
619
620 RCP<const Epetra_MultiVector> epetra_B;
621 RCP<Epetra_MultiVector> epetra_X;
622
623 const EpetraOperatorWrapper* opWrapper =
624 dynamic_cast<const EpetraOperatorWrapper*>(aztecOp);
625
626 if (opWrapper == 0) {
627 epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
628 epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
629 }
630
631 //
632 // Use AztecOO to solve each RHS one at a time (which is all that I can do anyway)
633 //
634
635 int totalIterations = 0;
636 SolveStatus<double> solveStatus;
637 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
638 solveStatus.achievedTol = -1.0;
639
640 /* Get the number of columns in the multivector. We use Thyra
641 * functions rather than Epetra functions to do this, as we
642 * might not yet have created an Epetra multivector. - KL */
643 //const int m = epetra_B->NumVectors();
644 const int m = B.domain()->dim();
645
646 for( int j = 0; j < m; ++j ) {
647
648 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: AztecOOLOWS:SingleSolve", SingleSolve);
649
650 //
651 // Get Epetra_Vector views of B(:,j) and X(:,j)
652 // How this is done will depend on whether we have a true Epetra operator
653 // or we are wrapping a general Thyra operator in an Epetra operator.
654 //
655
656 // We need to declare epetra_x_j as non-const because when we have a phony
657 // Epetra operator we'll have to copy a thyra vector into it.
658 RCP<Epetra_Vector> epetra_b_j;
659 RCP<Epetra_Vector> epetra_x_j;
660
661 if (opWrapper == 0) {
662 epetra_b_j = rcpFromRef(*const_cast<Epetra_Vector*>((*epetra_B)(j)));
663 epetra_x_j = rcpFromRef(*(*epetra_X)(j));
664 }
665 else {
666 if (is_null(epetra_b_j)) {
667 epetra_b_j = rcp(new Epetra_Vector(opRangeMap));
668 epetra_x_j = rcp(new Epetra_Vector(opDomainMap));
669 }
670 opWrapper->copyThyraIntoEpetra(*B.col(j), *epetra_b_j);
671 opWrapper->copyThyraIntoEpetra(*X->col(j), *epetra_x_j);
672 }
673
674 //
675 // Set the RHS and LHS
676 //
677
678 aztecSolver->SetRHS(&*epetra_b_j);
679 aztecSolver->SetLHS(&*epetra_x_j);
680
681 //
682 // Solve the linear system
683 //
684 timer.start(true);
685 {
686 SetAztecSolveState
687 setAztecSolveState(aztecSolver,out,verbLevel,solveMeasureType);
688 aztecSolver->Iterate( maxIterations, tol );
689 // NOTE: We ignore the returned status but get it below
690 }
691 timer.stop();
692
693 //
694 // Scale the solution
695 // (Originally, this was at the end of the loop after all columns had been
696 // processed. It's moved here because we need to do it before copying the
697 // solution back into a Thyra vector. - KL
698 //
699 if (aztecSolverScalar_ != 1.0)
700 epetra_x_j->Scale(1.0/aztecSolverScalar_);
701
702 //
703 // If necessary, convert the solution back to a non-epetra vector
704 //
705 if (opWrapper != 0) {
706 opWrapper->copyEpetraIntoThyra(*epetra_x_j, X->col(j).ptr());
707 }
708
709 //
710 // Set the return solve status
711 //
712
713 const int iterations = aztecSolver->NumIters();
714 const double achievedTol = aztecSolver->ScaledResidual();
715 const double *AZ_status = aztecSolver->GetAztecStatus();
716 std::ostringstream oss;
717 bool converged = false;
718 if (AZ_status[AZ_why]==AZ_normal) { oss << "Aztec returned AZ_normal."; converged = true; }
719 else if (AZ_status[AZ_why]==AZ_param) oss << "Aztec returned AZ_param.";
720 else if (AZ_status[AZ_why]==AZ_breakdown) oss << "Aztec returned AZ_breakdown.";
721 else if (AZ_status[AZ_why]==AZ_loss) oss << "Aztec returned AZ_loss.";
722 else if (AZ_status[AZ_why]==AZ_ill_cond) oss << "Aztec returned AZ_ill_cond.";
723 else if (AZ_status[AZ_why]==AZ_maxits) oss << "Aztec returned AZ_maxits.";
724 else oss << "Aztec returned an unknown status?";
725 oss << " Iterations = " << iterations << ".";
726 oss << " Achieved Tolerance = " << achievedTol << ".";
727 oss << " Total time = " << timer.totalElapsedTime() << " sec.";
728 if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE) && outputEveryRhs())
729 Teuchos::OSTab(out).o() << "j="<<j<<": " << oss.str() << "\n";
730
731 solveStatus.achievedTol = TEUCHOS_MAX(solveStatus.achievedTol, achievedTol);
732 // Note, achieveTol may actually be greater than tol due to ill conditioning and roundoff!
733
734 totalIterations += iterations;
735
736 solveStatus.message = oss.str();
737 if ( isDefaultSolveCriteria ) {
738 switch(solveStatus.solveStatus) {
739 case SOLVE_STATUS_UNKNOWN:
740 // Leave overall unknown!
741 break;
742 case SOLVE_STATUS_CONVERGED:
743 solveStatus.solveStatus = ( converged ? SOLVE_STATUS_CONVERGED : SOLVE_STATUS_UNCONVERGED );
744 break;
745 case SOLVE_STATUS_UNCONVERGED:
746 // Leave overall unconverged!
747 break;
748 default:
749 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
750 }
751 }
752 }
753
754 aztecSolver->UnsetLHSRHS();
755
756 //
757 // Release the Epetra_MultiVector views of X and B
758 //
759 epetra_X = Teuchos::null;
760 epetra_B = Teuchos::null;
761
762 //
763 // Update the overall solve criteria
764 //
765 totalTimer.stop();
766 SolveStatus<double> overallSolveStatus;
767 if (isDefaultSolveCriteria) {
768 overallSolveStatus.solveStatus = SOLVE_STATUS_UNKNOWN;
769 overallSolveStatus.achievedTol = SS::unknownTolerance();
770 }
771 else {
772 overallSolveStatus.solveStatus = solveStatus.solveStatus;
773 overallSolveStatus.achievedTol = solveStatus.achievedTol;
774 }
775 std::ostringstream oss;
776 oss
777 << "AztecOO solver "
778 << ( overallSolveStatus.solveStatus==SOLVE_STATUS_CONVERGED ? "converged" : "unconverged" )
779 << " on m = "<<m<<" RHSs using " << totalIterations << " cumulative iterations"
780 << " for an average of " << (totalIterations/m) << " iterations/RHS and"
781 << " total CPU time of "<<totalTimer.totalElapsedTime()<<" sec.";
782 overallSolveStatus.message = oss.str();
783
784 // Added these statistics following what was done for Belos
785 if (overallSolveStatus.extraParameters.is_null()) {
786 overallSolveStatus.extraParameters = Teuchos::parameterList ();
787 }
788 overallSolveStatus.extraParameters->set ("AztecOO/Iteration Count",
789 totalIterations);
790 // package independent version of the same
791 overallSolveStatus.extraParameters->set ("Iteration Count",
792 totalIterations);
793 overallSolveStatus.extraParameters->set ("AztecOO/Achieved Tolerance",
794 overallSolveStatus.achievedTol);
795
796 //
797 // Report the overall time
798 //
799 if (out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
800 *out
801 << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
802
803 return overallSolveStatus;
804
805}
806
807
808} // end namespace Thyra
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
RCP< const PreconditionerBase< double > > extract_prec()
Extract the preconditioner.
RCP< const LinearOpBase< double > > clone() const
void initialize(const RCP< const LinearOpBase< double > > &fwdOp, const RCP< const LinearOpSourceBase< double > > &fwdOpSrc, const RCP< const PreconditionerBase< double > > &prec, const bool isExternalPrec, const RCP< const LinearOpSourceBase< double > > &approxFwdOpSrc, const RCP< AztecOO > &aztecFwdSolver, const bool allowInexactFwdSolve=false, const RCP< AztecOO > &aztecAdjSolver=Teuchos::null, const bool allowInexactAdjSolve=false, const double aztecSolverScalar=1.0)
Sets up this object.
virtual bool opSupportedImpl(EOpTransp M_trans) const
AztecOOLinearOpWithSolve(const int fwdDefaultMaxIterations=400, const double fwdDefaultTol=1e-6, const int adjDefaultMaxIterations=400, const double adjDefaultTol=1e-6, const bool outputEveryRhs=false)
RCP< const LinearOpSourceBase< double > > extract_fwdOpSrc()
Extract the forward LinearOpBase<double> object so that it can be modified.
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
RCP< const VectorSpaceBase< double > > range() const
bool isExternalPrec() const
Determine if the preconditioner was external or not.
RCP< const LinearOpSourceBase< double > > extract_approxFwdOpSrc()
Extract the approximate forward LinearOpBase<double> object used to build the preconditioner.
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
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< const VectorSpaceBase< double > > domain() const
void uninitialize(RCP< const LinearOpBase< double > > *fwdOp=NULL, RCP< const LinearOpSourceBase< double > > *fwdOpSrc=NULL, RCP< const PreconditionerBase< double > > *prec=NULL, bool *isExternalPrec=NULL, RCP< const LinearOpSourceBase< double > > *approxFwdOpSrc=NULL, RCP< AztecOO > *aztecFwdSolver=NULL, bool *allowInexactFwdSolve=NULL, RCP< AztecOO > *aztecAdjSolver=NULL, bool *allowInexactAdjSolve=NULL, double *aztecSolverScalar=NULL)
Uninitialize.
virtual bool solveSupportsImpl(EOpTransp M_trans) const
NOTRANS

Generated for Stratimikos by doxygen 1.9.6