deal.II version GIT relicensing-2879-g253433159b 2025-03-21 00:10:00+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
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
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
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.
234 AssertDimension(x.locally_owned_size(),
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 {
260 AssertDimension(x.locally_owned_size(),
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"));
288 const int ierr = residual_vector->Norm2(&res_l2_norm);
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
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) :
320 if (CurrentIter == 0)
322
323 return status_test_collection->CheckStatus(CurrentIter,
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
343 {
344 return initial_residual;
345 }
346
347 double
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)
370 // Consider linear problem converged if any of the collection of
371 // criterion are met
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);
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);
404 status_test_rel_tol->DefineScaleForm(
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 *>(&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:
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
525 dynamic_cast<const internal::TrilinosReductionControl *>(
526 status_test.get()))
527 {
528 Assert(dynamic_cast<const ReductionControl *>(&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()));
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));
590 }
591
592
593 /* ---------------------- SolverCG ------------------------ */
594
598
599
600 /* ---------------------- SolverGMRES ------------------------ */
601
605
606
607 /* ---------------------- SolverBicgstab ------------------------ */
608
612
613
614 /* ---------------------- SolverCGS ------------------------ */
615
619
620
621 /* ---------------------- SolverTFQMR ------------------------ */
622
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.
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 "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();
699
700 verbose_cout << "Starting numeric factorization" << std::endl;
701 ierr = solver->NumericFactorization();
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.
747
748
749 verbose_cout << "Starting solve" << std::endl;
750
751 // Fetch return value of Amesos Solver functions
752 int ierr = solver->Solve();
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.
783
784 verbose_cout << "Starting solve" << std::endl;
785
786 // Fetch return value of Amesos Solver functions
787 int ierr = solver->Solve();
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.
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 "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();
826
827 verbose_cout << "Starting numeric factorization" << std::endl;
828 ierr = solver->NumericFactorization();
830
831 verbose_cout << "Starting solve" << std::endl;
832 ierr = solver->Solve();
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 Epetra_MultiVector &x,
850 const Epetra_MultiVector &b)
851 {
852 // We need an Epetra_LinearProblem object to let the Amesos solver know
853 // about the matrix and vectors.
855 std::make_unique<Epetra_LinearProblem>(const_cast<Epetra_Operator *>(&A),
856 &x,
857 const_cast<Epetra_MultiVector *>(
858 &b));
859
860 do_solve();
861 }
862
863
864 void
866 FullMatrix<double> &solution,
867 const FullMatrix<double> &rhs)
868 {
871 Assert(rhs.m() == solution.m(), ExcInternalError());
872 Assert(rhs.n() == solution.n(), ExcInternalError());
873
874 const unsigned int m = rhs.m();
875 const unsigned int n = rhs.n();
876
879
880 rhs_t.copy_transposed(rhs);
881 solution_t.copy_transposed(solution);
882
883 std::vector<double *> rhs_ptrs(n);
884 std::vector<double *> sultion_ptrs(n);
885
886 for (unsigned int i = 0; i < n; ++i)
887 {
888 rhs_ptrs[i] = &rhs_t[i][0];
889 sultion_ptrs[i] = &solution_t[i][0];
890 }
891
892 const Epetra_CrsMatrix &mat = sparse_matrix.trilinos_matrix();
893
894 Epetra_MultiVector trilinos_dst(View,
895 mat.OperatorRangeMap(),
896 sultion_ptrs.data(),
897 sultion_ptrs.size());
898 Epetra_MultiVector trilinos_src(View,
899 mat.OperatorDomainMap(),
900 rhs_ptrs.data(),
901 rhs_ptrs.size());
902
903 this->initialize(sparse_matrix);
904 this->solve(mat, trilinos_dst, trilinos_src);
905
906 solution.copy_transposed(solution_t);
907 }
908
909
910 void
912 MPI::Vector &x,
913 const MPI::Vector &b)
914 {
915 // We need an Epetra_LinearProblem object to let the Amesos solver know
916 // about the matrix and vectors.
917 linear_problem = std::make_unique<Epetra_LinearProblem>(
918 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
919 &x.trilinos_vector(),
920 const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
921
922 do_solve();
923 }
924
925
926
927 void
930 const ::Vector<double> &b)
931 {
932 // In case we call the solver with deal.II vectors, we create views of the
933 // vectors in Epetra format.
934 Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
935 Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
936 Assert(A.local_range().second == A.m(),
937 ExcMessage("Can only work in serial when using deal.II vectors."));
938 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
939 Epetra_Vector ep_b(View,
940 A.trilinos_matrix().RangeMap(),
941 const_cast<double *>(b.begin()));
942
943 // We need an Epetra_LinearProblem object to let the Amesos solver know
944 // about the matrix and vectors.
945 linear_problem = std::make_unique<Epetra_LinearProblem>(
946 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
947
948 do_solve();
949 }
950
951
952
953 void
955 const SparseMatrix &A,
957 const ::LinearAlgebra::distributed::Vector<double> &b)
958 {
959 AssertDimension(x.locally_owned_size(),
960 A.trilinos_matrix().DomainMap().NumMyElements());
961 AssertDimension(b.locally_owned_size(),
962 A.trilinos_matrix().RangeMap().NumMyElements());
963 Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
964 Epetra_Vector ep_b(View,
965 A.trilinos_matrix().RangeMap(),
966 const_cast<double *>(b.begin()));
967
968 // We need an Epetra_LinearProblem object to let the Amesos solver know
969 // about the matrix and vectors.
970 linear_problem = std::make_unique<Epetra_LinearProblem>(
971 const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
972
973 do_solve();
974 }
975} // namespace TrilinosWrappers
976
977
978// explicit instantiations
979// TODO: put these instantiations into generic file
980namespace TrilinosWrappers
981{
982 template void
983 SolverBase::do_solve(const PreconditionBase &preconditioner);
984
985 template void
986 SolverBase::do_solve(const Epetra_Operator &preconditioner);
987} // namespace TrilinosWrappers
988
990
991#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.
friend class Tensor
Definition tensor.h:865
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())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
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)
std::vector< index_type > data
Definition mpi.cc:746
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:190