Reference documentation for deal.II version GIT 3779fa9eb4 2023-09-28 13:00: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\}}\)
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 
38 namespace 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
179  ::Vector<double> &x,
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
209  ::Vector<double> &x,
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.locally_owned_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  {
262  A.OperatorDomainMap().NumMyElements());
263  AssertDimension(b.locally_owned_size(),
264  A.OperatorRangeMap().NumMyElements());
265 
266  Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
267  Epetra_Vector ep_b(View,
268  A.OperatorRangeMap(),
269  const_cast<double *>(b.begin()));
270 
271  // We need an Epetra_LinearProblem object to let the AztecOO solver know
272  // about the matrix and vectors.
273  linear_problem = std::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
274 
275  do_solve(preconditioner);
276  }
277 
278 
279  namespace internal
280  {
281  namespace
282  {
283  double
284  compute_residual(const Epetra_MultiVector *const residual_vector)
285  {
286  Assert(residual_vector->NumVectors() == 1,
287  ExcMessage("Residual multivector holds more than one vector"));
288  TrilinosScalar res_l2_norm = 0.0;
289  const int ierr = residual_vector->Norm2(&res_l2_norm);
290  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
291  return res_l2_norm;
292  }
293 
294  class TrilinosReductionControl : public AztecOO_StatusTest
295  {
296  public:
297  TrilinosReductionControl(const int max_steps,
298  const double tolerance,
299  const double reduction,
300  const Epetra_LinearProblem &linear_problem);
301 
302  virtual ~TrilinosReductionControl() override = default;
303 
304  virtual bool
305  ResidualVectorRequired() const override
306  {
307  return status_test_collection->ResidualVectorRequired();
308  }
309 
310  virtual AztecOO_StatusType
311  CheckStatus(int CurrentIter,
312  Epetra_MultiVector *CurrentResVector,
313  double CurrentResNormEst,
314  bool SolutionUpdated) override
315  {
316  // Note: CurrentResNormEst is set to -1.0 if no estimate of the
317  // residual value is available
319  (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
320  CurrentResNormEst);
321  if (CurrentIter == 0)
323 
324  return status_test_collection->CheckStatus(CurrentIter,
325  CurrentResVector,
326  CurrentResNormEst,
327  SolutionUpdated);
328  }
329 
330  virtual AztecOO_StatusType
331  GetStatus() const override
332  {
333  return status_test_collection->GetStatus();
334  }
335 
336  virtual std::ostream &
337  Print(std::ostream &stream, int indent = 0) const override
338  {
339  return status_test_collection->Print(stream, indent);
340  }
341 
342  double
343  get_initial_residual() const
344  {
345  return initial_residual;
346  }
347 
348  double
349  get_current_residual() const
350  {
351  return current_residual;
352  }
353 
354  private:
357  std::unique_ptr<AztecOO_StatusTestCombo> status_test_collection;
358  std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
359  std::unique_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
360  std::unique_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
361  };
362 
363 
364  TrilinosReductionControl::TrilinosReductionControl(
365  const int max_steps,
366  const double tolerance,
367  const double reduction,
368  const Epetra_LinearProblem &linear_problem)
369  : initial_residual(std::numeric_limits<double>::max())
370  , current_residual(std::numeric_limits<double>::max())
371  // Consider linear problem converged if any of the collection of
372  // criterion are met
373  , status_test_collection(std::make_unique<AztecOO_StatusTestCombo>(
374  AztecOO_StatusTestCombo::OR))
375  {
376  // Maximum number of iterations
377  Assert(max_steps >= 0, ExcInternalError());
379  std::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
381 
382  Assert(linear_problem.GetRHS()->NumVectors() == 1,
383  ExcMessage("RHS multivector holds more than one vector"));
384 
385  // Residual norm is below some absolute value
386  status_test_abs_tol = std::make_unique<AztecOO_StatusTestResNorm>(
387  *linear_problem.GetOperator(),
388  *(linear_problem.GetLHS()->operator()(0)),
389  *(linear_problem.GetRHS()->operator()(0)),
390  tolerance);
391  status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
392  AztecOO_StatusTestResNorm::TwoNorm);
393  status_test_abs_tol->DefineScaleForm(
394  AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
396 
397  // Residual norm, scaled by some initial value, is below some threshold
398  status_test_rel_tol = std::make_unique<AztecOO_StatusTestResNorm>(
399  *linear_problem.GetOperator(),
400  *(linear_problem.GetLHS()->operator()(0)),
401  *(linear_problem.GetRHS()->operator()(0)),
402  reduction);
403  status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
404  AztecOO_StatusTestResNorm::TwoNorm);
405  status_test_rel_tol->DefineScaleForm(
406  AztecOO_StatusTestResNorm::NormOfInitRes,
407  AztecOO_StatusTestResNorm::TwoNorm);
409  }
410 
411  } // namespace
412  } // namespace internal
413 
414 
415  template <typename Preconditioner>
416  void
417  SolverBase::do_solve(const Preconditioner &preconditioner)
418  {
419  int ierr;
420 
421  // Next we can allocate the AztecOO solver...
422  solver.SetProblem(*linear_problem);
423 
424  // ... and we can specify the solver to be used.
425  switch (solver_name)
426  {
427  case cg:
428  solver.SetAztecOption(AZ_solver, AZ_cg);
429  break;
430  case cgs:
431  solver.SetAztecOption(AZ_solver, AZ_cgs);
432  break;
433  case gmres:
434  solver.SetAztecOption(AZ_solver, AZ_gmres);
435  solver.SetAztecOption(AZ_kspace,
436  additional_data.gmres_restart_parameter);
437  break;
438  case bicgstab:
439  solver.SetAztecOption(AZ_solver, AZ_bicgstab);
440  break;
441  case tfqmr:
442  solver.SetAztecOption(AZ_solver, AZ_tfqmr);
443  break;
444  default:
445  Assert(false, ExcNotImplemented());
446  }
447 
448  // Set the preconditioner
449  set_preconditioner(solver, preconditioner);
450 
451  // ... set some options, ...
452  solver.SetAztecOption(AZ_output,
453  additional_data.output_solver_details ? AZ_all :
454  AZ_none);
455  solver.SetAztecOption(AZ_conv, AZ_noscaled);
456 
457  // By default, the Trilinos solver chooses convergence criterion based on
458  // the number of iterations made and an absolute tolerance.
459  // This implies that the use of the standard Trilinos convergence test
460  // actually coincides with ::IterationNumberControl because the
461  // solver, unless explicitly told otherwise, will Iterate() until a number
462  // of max_steps() are taken or an absolute tolerance() is attained.
463  // It is therefore suitable for use with both SolverControl or
464  // IterationNumberControl. The final check at the end will determine whether
465  // failure to converge to the defined residual norm constitutes failure
466  // (SolverControl) or is alright (IterationNumberControl).
467  // In the case that the SolverControl wants to perform ReductionControl,
468  // then we have to do a little extra something by prescribing a custom
469  // status test.
470  if (!status_test)
471  {
472  if (const ReductionControl *const reduction_control =
473  dynamic_cast<const ReductionControl *const>(&solver_control))
474  {
475  status_test = std::make_unique<internal::TrilinosReductionControl>(
476  reduction_control->max_steps(),
477  reduction_control->tolerance(),
478  reduction_control->reduction(),
479  *linear_problem);
480  solver.SetStatusTest(status_test.get());
481  }
482  }
483 
484  // ... and then solve!
485  ierr =
486  solver.Iterate(solver_control.max_steps(), solver_control.tolerance());
487 
488  // report errors in more detail than just by checking whether the return
489  // status is zero or greater. the error strings are taken from the
490  // implementation of the AztecOO::Iterate function
491  switch (ierr)
492  {
493  case -1:
494  AssertThrow(false,
495  ExcMessage("AztecOO::Iterate error code -1: "
496  "option not implemented"));
497  break;
498  case -2:
499  AssertThrow(false,
500  ExcMessage("AztecOO::Iterate error code -2: "
501  "numerical breakdown"));
502  break;
503  case -3:
504  AssertThrow(false,
505  ExcMessage("AztecOO::Iterate error code -3: "
506  "loss of precision"));
507  break;
508  case -4:
509  AssertThrow(false,
510  ExcMessage("AztecOO::Iterate error code -4: "
511  "GMRES Hessenberg ill-conditioned"));
512  break;
513  default:
514  AssertThrow(ierr >= 0, ExcTrilinosError(ierr));
515  }
516 
517  // Finally, let the deal.II SolverControl object know what has
518  // happened. If the solve succeeded, the status of the solver control will
519  // turn into SolverControl::success.
520  // If the residual is not computed/stored by the solver, as can happen for
521  // certain choices of solver or if a custom status test is set, then the
522  // result returned by TrueResidual() is equal to -1. In this case we must
523  // compute it ourself.
524  if (const internal::TrilinosReductionControl
525  *const reduction_control_status =
526  dynamic_cast<const internal::TrilinosReductionControl *const>(
527  status_test.get()))
528  {
529  Assert(dynamic_cast<const ReductionControl *const>(&solver_control),
530  ExcInternalError());
531 
532  // Check to see if solver converged in one step
533  // This can happen if the matrix is diagonal and a non-trivial
534  // preconditioner is used.
535  if (solver.NumIters() > 0)
536  {
537  // For ReductionControl, we must first register the initial residual
538  // value. This is the basis from which it will determine whether the
539  // current residual corresponds to a converged state.
540  solver_control.check(
541  0, reduction_control_status->get_initial_residual());
542  solver_control.check(
543  solver.NumIters(),
544  reduction_control_status->get_current_residual());
545  }
546  else
547  solver_control.check(
548  solver.NumIters(),
549  reduction_control_status->get_current_residual());
550  }
551  else
552  {
553  Assert(solver.TrueResidual() >= 0.0, ExcInternalError());
554  solver_control.check(solver.NumIters(), solver.TrueResidual());
555  }
556 
557  if (solver_control.last_check() != SolverControl::success)
558  AssertThrow(false,
559  SolverControl::NoConvergence(solver_control.last_step(),
560  solver_control.last_value()));
561  }
562 
563 
564 
565  template <>
566  void
567  SolverBase::set_preconditioner(AztecOO &solver,
568  const PreconditionBase &preconditioner)
569  {
570  // Introduce the preconditioner, if the identity preconditioner is used,
571  // the precondioner is set to none, ...
572  if (preconditioner.preconditioner.strong_count() != 0)
573  {
574  const int ierr = solver.SetPrecOperator(
575  const_cast<Epetra_Operator *>(preconditioner.preconditioner.get()));
576  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
577  }
578  else
579  solver.SetAztecOption(AZ_precond, AZ_none);
580  }
581 
582 
583  template <>
584  void
585  SolverBase::set_preconditioner(AztecOO &solver,
586  const Epetra_Operator &preconditioner)
587  {
588  const int ierr =
589  solver.SetPrecOperator(const_cast<Epetra_Operator *>(&preconditioner));
590  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
591  }
592 
593 
594  /* ---------------------- SolverCG ------------------------ */
595 
597  : SolverBase(cg, cn, data)
598  {}
599 
600 
601  /* ---------------------- SolverGMRES ------------------------ */
602 
604  : SolverBase(gmres, cn, data)
605  {}
606 
607 
608  /* ---------------------- SolverBicgstab ------------------------ */
609 
611  : SolverBase(bicgstab, cn, data)
612  {}
613 
614 
615  /* ---------------------- SolverCGS ------------------------ */
616 
618  : SolverBase(cgs, cn, data)
619  {}
620 
621 
622  /* ---------------------- SolverTFQMR ------------------------ */
623 
625  : SolverBase(tfqmr, cn, data)
626  {}
627 
628 
629 
630  /* ---------------------- SolverDirect ------------------------ */
631 
632  SolverDirect::AdditionalData::AdditionalData(const bool output_solver_details,
633  const std::string &solver_type)
634  : output_solver_details(output_solver_details)
635  , solver_type(solver_type)
636  {}
637 
638 
639 
641  : solver_control(cn)
642  , additional_data(data.output_solver_details, data.solver_type)
643  {}
644 
645 
646 
647  SolverControl &
649  {
650  return solver_control;
651  }
652 
653 
654 
655  void
657  {
658  // We need an Epetra_LinearProblem object to let the Amesos solver know
659  // about the matrix and vectors.
660  linear_problem = std::make_unique<Epetra_LinearProblem>();
661 
662  // Assign the matrix operator to the Epetra_LinearProblem object
663  linear_problem->SetOperator(
664  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()));
665 
666  // Fetch return value of Amesos Solver functions
667  int ierr;
668 
669  // First set whether we want to print the solver information to screen or
670  // not.
671  ConditionalOStream verbose_cout(std::cout,
673 
674  // Next allocate the Amesos solver, this is done in two steps, first we
675  // create a solver Factory and generate with that the concrete Amesos
676  // solver, if possible.
677  Amesos Factory;
678 
679  AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
680  ExcMessage(
681  std::string("You tried to select the solver type <") +
683  "> but this solver is not supported by Trilinos either "
684  "because it does not exist, or because Trilinos was not "
685  "configured for its use."));
686 
687  solver.reset(
688  Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
689 
690  verbose_cout << "Starting symbolic factorization" << std::endl;
691  ierr = solver->SymbolicFactorization();
692  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
693 
694  verbose_cout << "Starting numeric factorization" << std::endl;
695  ierr = solver->NumericFactorization();
696  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
697  }
698 
699 
700  void
702  {
703  // Assign the empty LHS vector to the Epetra_LinearProblem object
704  linear_problem->SetLHS(&x.trilinos_vector());
705 
706  // Assign the RHS vector to the Epetra_LinearProblem object
707  linear_problem->SetRHS(
708  const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
709 
710  // First set whether we want to print the solver information to screen or
711  // not.
712  ConditionalOStream verbose_cout(std::cout,
714 
715 
716  verbose_cout << "Starting solve" << std::endl;
717 
718  // Fetch return value of Amesos Solver functions
719  int ierr = solver->Solve();
720  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
721 
722  // Finally, force the SolverControl object to report convergence
723  solver_control.check(0, 0);
724  }
725 
726 
727 
728  void
731  const ::LinearAlgebra::distributed::Vector<double> &b)
732  {
733  Epetra_Vector ep_x(View,
734  linear_problem->GetOperator()->OperatorDomainMap(),
735  x.begin());
736  Epetra_Vector ep_b(View,
737  linear_problem->GetOperator()->OperatorRangeMap(),
738  const_cast<double *>(b.begin()));
739 
740  // Assign the empty LHS vector to the Epetra_LinearProblem object
741  linear_problem->SetLHS(&ep_x);
742 
743  // Assign the RHS vector to the Epetra_LinearProblem object
744  linear_problem->SetRHS(&ep_b);
745 
746  // First set whether we want to print the solver information to screen or
747  // not.
748  ConditionalOStream verbose_cout(std::cout,
750 
751  verbose_cout << "Starting solve" << std::endl;
752 
753  // Fetch return value of Amesos Solver functions
754  int ierr = solver->Solve();
755  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
756 
757  // Finally, force the SolverControl object to report convergence
758  solver_control.check(0, 0);
759  }
760 
761 
762 
763  void
765  {
766  // Fetch return value of Amesos Solver functions
767  int ierr;
768 
769  // First set whether we want to print the solver information to screen or
770  // not.
771  ConditionalOStream verbose_cout(std::cout,
773 
774  // Next allocate the Amesos solver, this is done in two steps, first we
775  // create a solver Factory and generate with that the concrete Amesos
776  // solver, if possible.
777  Amesos Factory;
778 
779  AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
780  ExcMessage(
781  std::string("You tried to select the solver type <") +
783  "> but this solver is not supported by Trilinos either "
784  "because it does not exist, or because Trilinos was not "
785  "configured for its use."));
786 
787  solver.reset(
788  Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
789 
790  verbose_cout << "Starting symbolic factorization" << std::endl;
791  ierr = solver->SymbolicFactorization();
792  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
793 
794  verbose_cout << "Starting numeric factorization" << std::endl;
795  ierr = solver->NumericFactorization();
796  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
797 
798  verbose_cout << "Starting solve" << std::endl;
799  ierr = solver->Solve();
800  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
801 
802  // Finally, let the deal.II SolverControl object know what has
803  // happened. If the solve succeeded, the status of the solver control will
804  // turn into SolverControl::success.
805  solver_control.check(0, 0);
806 
808  AssertThrow(false,
811  }
812 
813 
814  void
816  MPI::Vector &x,
817  const MPI::Vector &b)
818  {
819  // We need an Epetra_LinearProblem object to let the Amesos solver know
820  // about the matrix and vectors.
821  linear_problem = std::make_unique<Epetra_LinearProblem>(
822  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
823  &x.trilinos_vector(),
824  const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
825 
826  do_solve();
827  }
828 
829 
830 
831  void
833  ::Vector<double> &x,
834  const ::Vector<double> &b)
835  {
836  // In case we call the solver with deal.II vectors, we create views of the
837  // vectors in Epetra format.
838  Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
839  Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
840  Assert(A.local_range().second == A.m(),
841  ExcMessage("Can only work in serial when using deal.II vectors."));
842  Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
843  Epetra_Vector ep_b(View,
844  A.trilinos_matrix().RangeMap(),
845  const_cast<double *>(b.begin()));
846 
847  // We need an Epetra_LinearProblem object to let the Amesos solver know
848  // about the matrix and vectors.
849  linear_problem = std::make_unique<Epetra_LinearProblem>(
850  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
851 
852  do_solve();
853  }
854 
855 
856 
857  void
859  const SparseMatrix &A,
861  const ::LinearAlgebra::distributed::Vector<double> &b)
862  {
864  A.trilinos_matrix().DomainMap().NumMyElements());
865  AssertDimension(b.locally_owned_size(),
866  A.trilinos_matrix().RangeMap().NumMyElements());
867  Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
868  Epetra_Vector ep_b(View,
869  A.trilinos_matrix().RangeMap(),
870  const_cast<double *>(b.begin()));
871 
872  // We need an Epetra_LinearProblem object to let the Amesos solver know
873  // about the matrix and vectors.
874  linear_problem = std::make_unique<Epetra_LinearProblem>(
875  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
876 
877  do_solve();
878  }
879 } // namespace TrilinosWrappers
880 
881 
882 // explicit instantiations
883 // TODO: put these instantiations into generic file
884 namespace TrilinosWrappers
885 {
886  template void
887  SolverBase::do_solve(const PreconditionBase &preconditioner);
888 
889  template void
890  SolverBase::do_solve(const Epetra_Operator &preconditioner);
891 } // namespace TrilinosWrappers
892 
894 
895 #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())
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())
Definition: vector.h:110
virtual size_type size() const override
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
__global__ void reduction(Number *result, const Number *v, const size_type N)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
static const char A
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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