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