Reference documentation for deal.II version 9.4.1
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
trilinos_solver.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2020 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
17
18#ifdef DEAL_II_WITH_TRILINOS
19
21
25
26# include <AztecOO_StatusTest.h>
27# include <AztecOO_StatusTestCombo.h>
28# include <AztecOO_StatusTestMaxIters.h>
29# include <AztecOO_StatusTestResNorm.h>
30# include <AztecOO_StatusType.h>
31
32# include <cmath>
33# include <limits>
34# include <memory>
35
37
38namespace TrilinosWrappers
39{
41 const bool output_solver_details,
42 const unsigned int gmres_restart_parameter)
43 : output_solver_details(output_solver_details)
44 , gmres_restart_parameter(gmres_restart_parameter)
45 {}
46
47
48
51 , solver_control(cn)
52 , additional_data(data)
53 {}
54
55
56
58 SolverControl & cn,
59 const AdditionalData & data)
60 : solver_name(solver_name)
61 , solver_control(cn)
62 , additional_data(data)
63 {}
64
65
66
69 {
70 return solver_control;
71 }
72
73
74
75 void
77 MPI::Vector & x,
78 const MPI::Vector & b,
79 const PreconditionBase &preconditioner)
80 {
81 // We need an Epetra_LinearProblem object to let the AztecOO solver know
82 // about the matrix and vectors.
83 linear_problem = std::make_unique<Epetra_LinearProblem>(
84 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
85 &x.trilinos_vector(),
86 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
87
88 do_solve(preconditioner);
89 }
90
91
92
93 // Note: "A" is set as a constant reference so that all patterns for ::solve
94 // can be used by the inverse_operator of LinearOperator
95 void
97 MPI::Vector & x,
98 const MPI::Vector & b,
99 const PreconditionBase &preconditioner)
100 {
101 // We need an Epetra_LinearProblem object to let the AztecOO solver know
102 // about the matrix and vectors.
104 std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
105 &x.trilinos_vector(),
106 const_cast<Epetra_MultiVector *>(
107 &b.trilinos_vector()));
108
109 do_solve(preconditioner);
110 }
111
112
113
114 // Note: "A" is set as a constant reference so that all patterns for ::solve
115 // can be used by the inverse_operator of LinearOperator
116 void
118 MPI::Vector & x,
119 const MPI::Vector & b,
120 const Epetra_Operator &preconditioner)
121 {
122 // We need an Epetra_LinearProblem object to let the AztecOO solver know
123 // about the matrix and vectors.
125 std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
126 &x.trilinos_vector(),
127 const_cast<Epetra_MultiVector *>(
128 &b.trilinos_vector()));
129
130 do_solve(preconditioner);
131 }
132
133
134
135 // Note: "A" is set as a constant reference so that all patterns for ::solve
136 // can be used by the inverse_operator of LinearOperator
137 void
139 Epetra_MultiVector & x,
140 const Epetra_MultiVector &b,
141 const PreconditionBase & preconditioner)
142 {
143 // We need an Epetra_LinearProblem object to let the AztecOO solver know
144 // about the matrix and vectors.
146 std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
147 &x,
148 const_cast<Epetra_MultiVector *>(
149 &b));
150
151 do_solve(preconditioner);
152 }
153
154
155
156 // Note: "A" is set as a constant reference so that all patterns for ::solve
157 // can be used by the inverse_operator of LinearOperator
158 void
160 Epetra_MultiVector & x,
161 const Epetra_MultiVector &b,
162 const Epetra_Operator & preconditioner)
163 {
164 // We need an Epetra_LinearProblem object to let the AztecOO solver know
165 // about the matrix and vectors.
167 std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
168 &x,
169 const_cast<Epetra_MultiVector *>(
170 &b));
171
172 do_solve(preconditioner);
173 }
174
175
176
177 void
180 const ::Vector<double> &b,
181 const PreconditionBase & preconditioner)
182 {
183 // In case we call the solver with deal.II vectors, we create views of the
184 // vectors in Epetra format.
185 Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
186 Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
187 Assert(A.local_range().second == A.m(),
188 ExcMessage("Can only work in serial when using deal.II vectors."));
189 Assert(A.trilinos_matrix().Filled(),
190 ExcMessage("Matrix is not compressed. Call compress() method."));
191
192 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
193 Epetra_Vector ep_b(View,
194 A.trilinos_matrix().RangeMap(),
195 const_cast<double *>(b.begin()));
196
197 // We need an Epetra_LinearProblem object to let the AztecOO solver know
198 // about the matrix and vectors.
199 linear_problem = std::make_unique<Epetra_LinearProblem>(
200 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
201
202 do_solve(preconditioner);
203 }
204
205
206
207 void
210 const ::Vector<double> &b,
211 const PreconditionBase & preconditioner)
212 {
213 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
214 Epetra_Vector ep_b(View,
215 A.OperatorRangeMap(),
216 const_cast<double *>(b.begin()));
217
218 // We need an Epetra_LinearProblem object to let the AztecOO solver know
219 // about the matrix and vectors.
220 linear_problem = std::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
221
222 do_solve(preconditioner);
223 }
224
225
226
227 void
230 const ::LinearAlgebra::distributed::Vector<double> &b,
231 const PreconditionBase &preconditioner)
232 {
233 // In case we call the solver with deal.II vectors, we create views of the
234 // vectors in Epetra format.
236 A.trilinos_matrix().DomainMap().NumMyElements());
237 AssertDimension(b.local_size(),
238 A.trilinos_matrix().RangeMap().NumMyElements());
239
240 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
241 Epetra_Vector ep_b(View,
242 A.trilinos_matrix().RangeMap(),
243 const_cast<double *>(b.begin()));
244
245 // We need an Epetra_LinearProblem object to let the AztecOO solver know
246 // about the matrix and vectors.
247 linear_problem = std::make_unique<Epetra_LinearProblem>(
248 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
249
250 do_solve(preconditioner);
251 }
252
253
254
255 void
258 const ::LinearAlgebra::distributed::Vector<double> &b,
259 const PreconditionBase &preconditioner)
260 {
261 AssertDimension(x.local_size(), A.OperatorDomainMap().NumMyElements());
262 AssertDimension(b.local_size(), A.OperatorRangeMap().NumMyElements());
263
264 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
265 Epetra_Vector ep_b(View,
266 A.OperatorRangeMap(),
267 const_cast<double *>(b.begin()));
268
269 // We need an Epetra_LinearProblem object to let the AztecOO solver know
270 // about the matrix and vectors.
271 linear_problem = std::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
272
273 do_solve(preconditioner);
274 }
275
276
277 namespace internal
278 {
279 namespace
280 {
281 double
282 compute_residual(const Epetra_MultiVector *const residual_vector)
283 {
284 Assert(residual_vector->NumVectors() == 1,
285 ExcMessage("Residual multivector holds more than one vector"));
286 TrilinosScalar res_l2_norm = 0.0;
287 const int ierr = residual_vector->Norm2(&res_l2_norm);
288 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
289 return res_l2_norm;
290 }
291
292 class TrilinosReductionControl : public AztecOO_StatusTest
293 {
294 public:
295 TrilinosReductionControl(const int max_steps,
296 const double tolerance,
297 const double reduction,
298 const Epetra_LinearProblem &linear_problem);
299
300 virtual ~TrilinosReductionControl() override = default;
301
302 virtual bool
303 ResidualVectorRequired() const override
304 {
305 return status_test_collection->ResidualVectorRequired();
306 }
307
308 virtual AztecOO_StatusType
309 CheckStatus(int CurrentIter,
310 Epetra_MultiVector *CurrentResVector,
311 double CurrentResNormEst,
312 bool SolutionUpdated) override
313 {
314 // Note: CurrentResNormEst is set to -1.0 if no estimate of the
315 // residual value is available
317 (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
318 CurrentResNormEst);
319 if (CurrentIter == 0)
321
322 return status_test_collection->CheckStatus(CurrentIter,
323 CurrentResVector,
324 CurrentResNormEst,
325 SolutionUpdated);
326 }
327
328 virtual AztecOO_StatusType
329 GetStatus() const override
330 {
331 return status_test_collection->GetStatus();
332 }
333
334 virtual std::ostream &
335 Print(std::ostream &stream, int indent = 0) const override
336 {
337 return status_test_collection->Print(stream, indent);
338 }
339
340 double
341 get_initial_residual() const
342 {
343 return initial_residual;
344 }
345
346 double
347 get_current_residual() const
348 {
349 return current_residual;
350 }
351
352 private:
355 std::unique_ptr<AztecOO_StatusTestCombo> status_test_collection;
356 std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
357 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
358 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
359 };
360
361
362 TrilinosReductionControl::TrilinosReductionControl(
363 const int max_steps,
364 const double tolerance,
365 const double reduction,
366 const Epetra_LinearProblem &linear_problem)
367 : initial_residual(std::numeric_limits<double>::max())
368 , current_residual(std::numeric_limits<double>::max())
369 // Consider linear problem converged if any of the collection of
370 // criterion are met
371 , status_test_collection(std::make_unique<AztecOO_StatusTestCombo>(
372 AztecOO_StatusTestCombo::OR))
373 {
374 // Maximum number of iterations
375 Assert(max_steps >= 0, ExcInternalError());
377 std::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
379
380 Assert(linear_problem.GetRHS()->NumVectors() == 1,
381 ExcMessage("RHS multivector holds more than one vector"));
382
383 // Residual norm is below some absolute value
384 status_test_abs_tol = std::make_unique<AztecOO_StatusTestResNorm>(
385 *linear_problem.GetOperator(),
386 *(linear_problem.GetLHS()->operator()(0)),
387 *(linear_problem.GetRHS()->operator()(0)),
388 tolerance);
389 status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
390 AztecOO_StatusTestResNorm::TwoNorm);
391 status_test_abs_tol->DefineScaleForm(
392 AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
394
395 // Residual norm, scaled by some initial value, is below some threshold
396 status_test_rel_tol = std::make_unique<AztecOO_StatusTestResNorm>(
397 *linear_problem.GetOperator(),
398 *(linear_problem.GetLHS()->operator()(0)),
399 *(linear_problem.GetRHS()->operator()(0)),
400 reduction);
401 status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
402 AztecOO_StatusTestResNorm::TwoNorm);
403 status_test_rel_tol->DefineScaleForm(
404 AztecOO_StatusTestResNorm::NormOfInitRes,
405 AztecOO_StatusTestResNorm::TwoNorm);
407 }
408
409 } // namespace
410 } // namespace internal
411
412
413 template <typename Preconditioner>
414 void
415 SolverBase::do_solve(const Preconditioner &preconditioner)
416 {
417 int ierr;
418
419 // Next we can allocate the AztecOO solver...
420 solver.SetProblem(*linear_problem);
421
422 // ... and we can specify the solver to be used.
423 switch (solver_name)
424 {
425 case cg:
426 solver.SetAztecOption(AZ_solver, AZ_cg);
427 break;
428 case cgs:
429 solver.SetAztecOption(AZ_solver, AZ_cgs);
430 break;
431 case gmres:
432 solver.SetAztecOption(AZ_solver, AZ_gmres);
433 solver.SetAztecOption(AZ_kspace,
434 additional_data.gmres_restart_parameter);
435 break;
436 case bicgstab:
437 solver.SetAztecOption(AZ_solver, AZ_bicgstab);
438 break;
439 case tfqmr:
440 solver.SetAztecOption(AZ_solver, AZ_tfqmr);
441 break;
442 default:
443 Assert(false, ExcNotImplemented());
444 }
445
446 // Set the preconditioner
447 set_preconditioner(solver, preconditioner);
448
449 // ... set some options, ...
450 solver.SetAztecOption(AZ_output,
451 additional_data.output_solver_details ? AZ_all :
452 AZ_none);
453 solver.SetAztecOption(AZ_conv, AZ_noscaled);
454
455 // By default, the Trilinos solver chooses convergence criterion based on
456 // the number of iterations made and an absolute tolerance.
457 // This implies that the use of the standard Trilinos convergence test
458 // actually coincides with ::IterationNumberControl because the
459 // solver, unless explicitly told otherwise, will Iterate() until a number
460 // of max_steps() are taken or an absolute tolerance() is attained.
461 // It is therefore suitable for use with both SolverControl or
462 // IterationNumberControl. The final check at the end will determine whether
463 // failure to converge to the defined residual norm constitutes failure
464 // (SolverControl) or is alright (IterationNumberControl).
465 // In the case that the SolverControl wants to perform ReductionControl,
466 // then we have to do a little extra something by prescribing a custom
467 // status test.
468 if (!status_test)
469 {
470 if (const ReductionControl *const reduction_control =
471 dynamic_cast<const ReductionControl *const>(&solver_control))
472 {
473 status_test = std::make_unique<internal::TrilinosReductionControl>(
474 reduction_control->max_steps(),
475 reduction_control->tolerance(),
476 reduction_control->reduction(),
477 *linear_problem);
478 solver.SetStatusTest(status_test.get());
479 }
480 }
481
482 // ... and then solve!
483 ierr =
484 solver.Iterate(solver_control.max_steps(), solver_control.tolerance());
485
486 // report errors in more detail than just by checking whether the return
487 // status is zero or greater. the error strings are taken from the
488 // implementation of the AztecOO::Iterate function
489 switch (ierr)
490 {
491 case -1:
492 AssertThrow(false,
493 ExcMessage("AztecOO::Iterate error code -1: "
494 "option not implemented"));
495 break;
496 case -2:
497 AssertThrow(false,
498 ExcMessage("AztecOO::Iterate error code -2: "
499 "numerical breakdown"));
500 break;
501 case -3:
502 AssertThrow(false,
503 ExcMessage("AztecOO::Iterate error code -3: "
504 "loss of precision"));
505 break;
506 case -4:
507 AssertThrow(false,
508 ExcMessage("AztecOO::Iterate error code -4: "
509 "GMRES Hessenberg ill-conditioned"));
510 break;
511 default:
512 AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
513 }
514
515 // Finally, let the deal.II SolverControl object know what has
516 // happened. If the solve succeeded, the status of the solver control will
517 // turn into SolverControl::success.
518 // If the residual is not computed/stored by the solver, as can happen for
519 // certain choices of solver or if a custom status test is set, then the
520 // result returned by TrueResidual() is equal to -1. In this case we must
521 // compute it ourself.
522 if (const internal::TrilinosReductionControl
523 *const reduction_control_status =
524 dynamic_cast<const internal::TrilinosReductionControl *const>(
525 status_test.get()))
526 {
527 Assert(dynamic_cast<const ReductionControl *const>(&solver_control),
529
530 // Check to see if solver converged in one step
531 // This can happen if the matrix is diagonal and a non-trivial
532 // preconditioner is used.
533 if (solver.NumIters() > 0)
534 {
535 // For ReductionControl, we must first register the initial residual
536 // value. This is the basis from which it will determine whether the
537 // current residual corresponds to a converged state.
538 solver_control.check(
539 0, reduction_control_status->get_initial_residual());
540 solver_control.check(
541 solver.NumIters(),
542 reduction_control_status->get_current_residual());
543 }
544 else
545 solver_control.check(
546 solver.NumIters(),
547 reduction_control_status->get_current_residual());
548 }
549 else
550 {
551 Assert(solver.TrueResidual() >= 0.0, ExcInternalError());
552 solver_control.check(solver.NumIters(), solver.TrueResidual());
553 }
554
555 if (solver_control.last_check() != SolverControl::success)
556 AssertThrow(false,
557 SolverControl::NoConvergence(solver_control.last_step(),
558 solver_control.last_value()));
559 }
560
561
562
563 template <>
564 void
565 SolverBase::set_preconditioner(AztecOO & solver,
566 const PreconditionBase &preconditioner)
567 {
568 // Introduce the preconditioner, if the identity preconditioner is used,
569 // the precondioner is set to none, ...
570 if (preconditioner.preconditioner.strong_count() != 0)
571 {
572 const int ierr = solver.SetPrecOperator(
573 const_cast<Epetra_Operator *>(preconditioner.preconditioner.get()));
574 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
575 }
576 else
577 solver.SetAztecOption(AZ_precond, AZ_none);
578 }
579
580
581 template <>
582 void
583 SolverBase::set_preconditioner(AztecOO & solver,
584 const Epetra_Operator &preconditioner)
585 {
586 const int ierr =
587 solver.SetPrecOperator(const_cast<Epetra_Operator *>(&preconditioner));
588 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
589 }
590
591
592 /* ---------------------- SolverCG ------------------------ */
593
594 SolverCG::AdditionalData::AdditionalData(const bool output_solver_details)
595 : SolverBase::AdditionalData(output_solver_details)
596 {}
597
598
599
601 : SolverBase(cn, data)
602 , additional_data(data)
603 {
604 solver_name = cg;
605 }
606
607
608 /* ---------------------- SolverGMRES ------------------------ */
609
611 const bool output_solver_details,
612 const unsigned int restart_parameter)
613 : SolverBase::AdditionalData(output_solver_details, restart_parameter)
614 {}
615
616
617
619 : SolverBase(cn, data)
620 , additional_data(data)
621 {
623 }
624
625
626 /* ---------------------- SolverBicgstab ------------------------ */
627
629 const bool output_solver_details)
630 : SolverBase::AdditionalData(output_solver_details)
631 {}
632
633
634
636 : SolverBase(cn, data)
637 , additional_data(data)
638 {
640 }
641
642
643 /* ---------------------- SolverCGS ------------------------ */
644
645 SolverCGS::AdditionalData::AdditionalData(const bool output_solver_details)
646 : SolverBase::AdditionalData(output_solver_details)
647 {}
648
649
650
652 : SolverBase(cn, data)
653 , additional_data(data)
654 {
656 }
657
658
659 /* ---------------------- SolverTFQMR ------------------------ */
660
661 SolverTFQMR::AdditionalData::AdditionalData(const bool output_solver_details)
662 : SolverBase::AdditionalData(output_solver_details)
663 {}
664
665
666
668 : SolverBase(cn, data)
669 , additional_data(data)
670 {
672 }
673
674
675
676 /* ---------------------- SolverDirect ------------------------ */
677
678 SolverDirect::AdditionalData::AdditionalData(const bool output_solver_details,
679 const std::string &solver_type)
680 : output_solver_details(output_solver_details)
681 , solver_type(solver_type)
682 {}
683
684
685
687 : solver_control(cn)
688 , additional_data(data.output_solver_details, data.solver_type)
689 {}
690
691
692
695 {
696 return solver_control;
697 }
698
699
700
701 void
703 {
704 // We need an Epetra_LinearProblem object to let the Amesos solver know
705 // about the matrix and vectors.
706 linear_problem = std::make_unique<Epetra_LinearProblem>();
707
708 // Assign the matrix operator to the Epetra_LinearProblem object
709 linear_problem->SetOperator(
710 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()));
711
712 // Fetch return value of Amesos Solver functions
713 int ierr;
714
715 // First set whether we want to print the solver information to screen or
716 // not.
717 ConditionalOStream verbose_cout(std::cout,
719
720 // Next allocate the Amesos solver, this is done in two steps, first we
721 // create a solver Factory and generate with that the concrete Amesos
722 // solver, if possible.
723 Amesos Factory;
724
725 AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
727 std::string("You tried to select the solver type <") +
729 "> but this solver is not supported by Trilinos either "
730 "because it does not exist, or because Trilinos was not "
731 "configured for its use."));
732
733 solver.reset(
734 Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
735
736 verbose_cout << "Starting symbolic factorization" << std::endl;
737 ierr = solver->SymbolicFactorization();
738 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
739
740 verbose_cout << "Starting numeric factorization" << std::endl;
741 ierr = solver->NumericFactorization();
742 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
743 }
744
745
746 void
748 {
749 // Assign the empty LHS vector to the Epetra_LinearProblem object
750 linear_problem->SetLHS(&x.trilinos_vector());
751
752 // Assign the RHS vector to the Epetra_LinearProblem object
753 linear_problem->SetRHS(
754 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
755
756 // First set whether we want to print the solver information to screen or
757 // not.
758 ConditionalOStream verbose_cout(std::cout,
760
761
762 verbose_cout << "Starting solve" << std::endl;
763
764 // Fetch return value of Amesos Solver functions
765 int ierr = solver->Solve();
766 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
767
768 // Finally, force the SolverControl object to report convergence
769 solver_control.check(0, 0);
770 }
771
772
773
774 void
777 const ::LinearAlgebra::distributed::Vector<double> &b)
778 {
779 Epetra_Vector ep_x(View,
780 linear_problem->GetOperator()->OperatorDomainMap(),
781 x.begin());
782 Epetra_Vector ep_b(View,
783 linear_problem->GetOperator()->OperatorRangeMap(),
784 const_cast<double *>(b.begin()));
785
786 // Assign the empty LHS vector to the Epetra_LinearProblem object
787 linear_problem->SetLHS(&ep_x);
788
789 // Assign the RHS vector to the Epetra_LinearProblem object
790 linear_problem->SetRHS(&ep_b);
791
792 // First set whether we want to print the solver information to screen or
793 // not.
794 ConditionalOStream verbose_cout(std::cout,
796
797 verbose_cout << "Starting solve" << std::endl;
798
799 // Fetch return value of Amesos Solver functions
800 int ierr = solver->Solve();
801 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
802
803 // Finally, force the SolverControl object to report convergence
804 solver_control.check(0, 0);
805 }
806
807
808
809 void
811 {
812 // Fetch return value of Amesos Solver functions
813 int ierr;
814
815 // First set whether we want to print the solver information to screen or
816 // not.
817 ConditionalOStream verbose_cout(std::cout,
819
820 // Next allocate the Amesos solver, this is done in two steps, first we
821 // create a solver Factory and generate with that the concrete Amesos
822 // solver, if possible.
823 Amesos Factory;
824
825 AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
827 std::string("You tried to select the solver type <") +
829 "> but this solver is not supported by Trilinos either "
830 "because it does not exist, or because Trilinos was not "
831 "configured for its use."));
832
833 solver.reset(
834 Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
835
836 verbose_cout << "Starting symbolic factorization" << std::endl;
837 ierr = solver->SymbolicFactorization();
838 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
839
840 verbose_cout << "Starting numeric factorization" << std::endl;
841 ierr = solver->NumericFactorization();
842 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
843
844 verbose_cout << "Starting solve" << std::endl;
845 ierr = solver->Solve();
846 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
847
848 // Finally, let the deal.II SolverControl object know what has
849 // happened. If the solve succeeded, the status of the solver control will
850 // turn into SolverControl::success.
851 solver_control.check(0, 0);
852
854 AssertThrow(false,
857 }
858
859
860 void
862 MPI::Vector & x,
863 const MPI::Vector & b)
864 {
865 // We need an Epetra_LinearProblem object to let the Amesos solver know
866 // about the matrix and vectors.
867 linear_problem = std::make_unique<Epetra_LinearProblem>(
868 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
869 &x.trilinos_vector(),
870 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
871
872 do_solve();
873 }
874
875
876
877 void
880 const ::Vector<double> &b)
881 {
882 // In case we call the solver with deal.II vectors, we create views of the
883 // vectors in Epetra format.
884 Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
885 Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
886 Assert(A.local_range().second == A.m(),
887 ExcMessage("Can only work in serial when using deal.II vectors."));
888 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
889 Epetra_Vector ep_b(View,
890 A.trilinos_matrix().RangeMap(),
891 const_cast<double *>(b.begin()));
892
893 // We need an Epetra_LinearProblem object to let the Amesos solver know
894 // about the matrix and vectors.
895 linear_problem = std::make_unique<Epetra_LinearProblem>(
896 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
897
898 do_solve();
899 }
900
901
902
903 void
905 const SparseMatrix & A,
907 const ::LinearAlgebra::distributed::Vector<double> &b)
908 {
910 A.trilinos_matrix().DomainMap().NumMyElements());
911 AssertDimension(b.local_size(),
912 A.trilinos_matrix().RangeMap().NumMyElements());
913 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
914 Epetra_Vector ep_b(View,
915 A.trilinos_matrix().RangeMap(),
916 const_cast<double *>(b.begin()));
917
918 // We need an Epetra_LinearProblem object to let the Amesos solver know
919 // about the matrix and vectors.
920 linear_problem = std::make_unique<Epetra_LinearProblem>(
921 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
922
923 do_solve();
924 }
925} // namespace TrilinosWrappers
926
927
928// explicit instantiations
929// TODO: put these instantiations into generic file
930namespace TrilinosWrappers
931{
932 template void
933 SolverBase::do_solve(const PreconditionBase &preconditioner);
934
935 template void
936 SolverBase::do_solve(const Epetra_Operator &preconditioner);
937} // namespace TrilinosWrappers
938
940
941#endif // DEAL_II_WITH_PETSC
State last_check() const
unsigned int last_step() const
double last_value() const
virtual State check(const unsigned int step, const double check_value)
@ success
Stop iteration, goal reached.
SolverControl & control() const
void solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b, const PreconditionBase &preconditioner)
const AdditionalData additional_data
void do_solve(const Preconditioner &preconditioner)
std::unique_ptr< Epetra_LinearProblem > linear_problem
SolverBase(SolverControl &cn, const AdditionalData &data=AdditionalData())
enum TrilinosWrappers::SolverBase::SolverName solver_name
const AdditionalData additional_data
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
void initialize(const SparseMatrix &A)
std::unique_ptr< Amesos_BaseSolver > solver
void solve(MPI::Vector &x, const MPI::Vector &b)
std::unique_ptr< Epetra_LinearProblem > linear_problem
SolverControl & control() const
const AdditionalData additional_data
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
const Epetra_CrsMatrix & trilinos_matrix() const
std::pair< size_type, size_type > local_range() const
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
__global__ void reduction(Number *result, const Number *v, const size_type N)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
const Epetra_MultiVector & trilinos_vector() const
Teuchos::RCP< Epetra_Operator > preconditioner
size_type size() const
iterator begin()
STL namespace.
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
AdditionalData(const bool output_solver_details=false)
AdditionalData(const bool output_solver_details=false)
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
AdditionalData(const bool output_solver_details=false, const unsigned int restart_parameter=30)
AdditionalData(const bool output_solver_details=false)
double initial_residual
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_abs_tol
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_rel_tol
std::unique_ptr< AztecOO_StatusTestMaxIters > status_test_max_steps
std::unique_ptr< AztecOO_StatusTestCombo > status_test_collection
double current_residual
double TrilinosScalar
Definition: types.h:163