Reference documentation for deal.II version 9.5.0
\(\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 - 2023 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
595 : SolverBase(cg, cn, data)
596 {}
597
598
599 /* ---------------------- SolverGMRES ------------------------ */
600
602 : SolverBase(gmres, cn, data)
603 {}
604
605
606 /* ---------------------- SolverBicgstab ------------------------ */
607
609 : SolverBase(bicgstab, cn, data)
610 {}
611
612
613 /* ---------------------- SolverCGS ------------------------ */
614
616 : SolverBase(cgs, cn, data)
617 {}
618
619
620 /* ---------------------- SolverTFQMR ------------------------ */
621
623 : SolverBase(tfqmr, cn, data)
624 {}
625
626
627
628 /* ---------------------- SolverDirect ------------------------ */
629
630 SolverDirect::AdditionalData::AdditionalData(const bool output_solver_details,
631 const std::string &solver_type)
632 : output_solver_details(output_solver_details)
633 , solver_type(solver_type)
634 {}
635
636
637
639 : solver_control(cn)
640 , additional_data(data.output_solver_details, data.solver_type)
641 {}
642
643
644
647 {
648 return solver_control;
649 }
650
651
652
653 void
655 {
656 // We need an Epetra_LinearProblem object to let the Amesos solver know
657 // about the matrix and vectors.
658 linear_problem = std::make_unique<Epetra_LinearProblem>();
659
660 // Assign the matrix operator to the Epetra_LinearProblem object
661 linear_problem->SetOperator(
662 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()));
663
664 // Fetch return value of Amesos Solver functions
665 int ierr;
666
667 // First set whether we want to print the solver information to screen or
668 // not.
669 ConditionalOStream verbose_cout(std::cout,
671
672 // Next allocate the Amesos solver, this is done in two steps, first we
673 // create a solver Factory and generate with that the concrete Amesos
674 // solver, if possible.
675 Amesos Factory;
676
677 AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
679 std::string("You tried to select the solver type <") +
681 "> but this solver is not supported by Trilinos either "
682 "because it does not exist, or because Trilinos was not "
683 "configured for its use."));
684
685 solver.reset(
686 Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
687
688 verbose_cout << "Starting symbolic factorization" << std::endl;
689 ierr = solver->SymbolicFactorization();
690 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
691
692 verbose_cout << "Starting numeric factorization" << std::endl;
693 ierr = solver->NumericFactorization();
694 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
695 }
696
697
698 void
700 {
701 // Assign the empty LHS vector to the Epetra_LinearProblem object
702 linear_problem->SetLHS(&x.trilinos_vector());
703
704 // Assign the RHS vector to the Epetra_LinearProblem object
705 linear_problem->SetRHS(
706 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
707
708 // First set whether we want to print the solver information to screen or
709 // not.
710 ConditionalOStream verbose_cout(std::cout,
712
713
714 verbose_cout << "Starting solve" << std::endl;
715
716 // Fetch return value of Amesos Solver functions
717 int ierr = solver->Solve();
718 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
719
720 // Finally, force the SolverControl object to report convergence
721 solver_control.check(0, 0);
722 }
723
724
725
726 void
729 const ::LinearAlgebra::distributed::Vector<double> &b)
730 {
731 Epetra_Vector ep_x(View,
732 linear_problem->GetOperator()->OperatorDomainMap(),
733 x.begin());
734 Epetra_Vector ep_b(View,
735 linear_problem->GetOperator()->OperatorRangeMap(),
736 const_cast<double *>(b.begin()));
737
738 // Assign the empty LHS vector to the Epetra_LinearProblem object
739 linear_problem->SetLHS(&ep_x);
740
741 // Assign the RHS vector to the Epetra_LinearProblem object
742 linear_problem->SetRHS(&ep_b);
743
744 // First set whether we want to print the solver information to screen or
745 // not.
746 ConditionalOStream verbose_cout(std::cout,
748
749 verbose_cout << "Starting solve" << std::endl;
750
751 // Fetch return value of Amesos Solver functions
752 int ierr = solver->Solve();
753 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
754
755 // Finally, force the SolverControl object to report convergence
756 solver_control.check(0, 0);
757 }
758
759
760
761 void
763 {
764 // Fetch return value of Amesos Solver functions
765 int ierr;
766
767 // First set whether we want to print the solver information to screen or
768 // not.
769 ConditionalOStream verbose_cout(std::cout,
771
772 // Next allocate the Amesos solver, this is done in two steps, first we
773 // create a solver Factory and generate with that the concrete Amesos
774 // solver, if possible.
775 Amesos Factory;
776
777 AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
779 std::string("You tried to select the solver type <") +
781 "> but this solver is not supported by Trilinos either "
782 "because it does not exist, or because Trilinos was not "
783 "configured for its use."));
784
785 solver.reset(
786 Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
787
788 verbose_cout << "Starting symbolic factorization" << std::endl;
789 ierr = solver->SymbolicFactorization();
790 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
791
792 verbose_cout << "Starting numeric factorization" << std::endl;
793 ierr = solver->NumericFactorization();
794 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
795
796 verbose_cout << "Starting solve" << std::endl;
797 ierr = solver->Solve();
798 AssertThrow(ierr == 0, ExcTrilinosError(ierr));
799
800 // Finally, let the deal.II SolverControl object know what has
801 // happened. If the solve succeeded, the status of the solver control will
802 // turn into SolverControl::success.
803 solver_control.check(0, 0);
804
806 AssertThrow(false,
809 }
810
811
812 void
814 MPI::Vector & x,
815 const MPI::Vector & b)
816 {
817 // We need an Epetra_LinearProblem object to let the Amesos solver know
818 // about the matrix and vectors.
819 linear_problem = std::make_unique<Epetra_LinearProblem>(
820 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
821 &x.trilinos_vector(),
822 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
823
824 do_solve();
825 }
826
827
828
829 void
832 const ::Vector<double> &b)
833 {
834 // In case we call the solver with deal.II vectors, we create views of the
835 // vectors in Epetra format.
836 Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
837 Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
838 Assert(A.local_range().second == A.m(),
839 ExcMessage("Can only work in serial when using deal.II vectors."));
840 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
841 Epetra_Vector ep_b(View,
842 A.trilinos_matrix().RangeMap(),
843 const_cast<double *>(b.begin()));
844
845 // We need an Epetra_LinearProblem object to let the Amesos solver know
846 // about the matrix and vectors.
847 linear_problem = std::make_unique<Epetra_LinearProblem>(
848 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
849
850 do_solve();
851 }
852
853
854
855 void
857 const SparseMatrix & A,
859 const ::LinearAlgebra::distributed::Vector<double> &b)
860 {
862 A.trilinos_matrix().DomainMap().NumMyElements());
863 AssertDimension(b.local_size(),
864 A.trilinos_matrix().RangeMap().NumMyElements());
865 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
866 Epetra_Vector ep_b(View,
867 A.trilinos_matrix().RangeMap(),
868 const_cast<double *>(b.begin()));
869
870 // We need an Epetra_LinearProblem object to let the Amesos solver know
871 // about the matrix and vectors.
872 linear_problem = std::make_unique<Epetra_LinearProblem>(
873 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
874
875 do_solve();
876 }
877} // namespace TrilinosWrappers
878
879
880// explicit instantiations
881// TODO: put these instantiations into generic file
882namespace TrilinosWrappers
883{
884 template void
885 SolverBase::do_solve(const PreconditionBase &preconditioner);
886
887 template void
888 SolverBase::do_solve(const Epetra_Operator &preconditioner);
889} // namespace TrilinosWrappers
890
892
893#endif // DEAL_II_WITH_PETSC
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
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.
const Epetra_MultiVector & trilinos_vector() const
Teuchos::RCP< Epetra_Operator > preconditioner
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
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverCGS(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
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
const Epetra_CrsMatrix & trilinos_matrix() const
std::pair< size_type, size_type > local_range() const
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
__global__ void reduction(Number *result, const Number *v, const size_type N)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
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)
STL namespace.
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
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:175