Reference documentation for deal.II version GIT 5983d193e2 2023-05-27 16:50: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 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
17 
18 #ifdef DEAL_II_WITH_TRILINOS
19 
21 
25 
26 # include <AztecOO_StatusTest.h>
27 # include <AztecOO_StatusTestCombo.h>
28 # include <AztecOO_StatusTestMaxIters.h>
29 # include <AztecOO_StatusTestResNorm.h>
30 # include <AztecOO_StatusType.h>
31 
32 # include <cmath>
33 # include <limits>
34 # include <memory>
35 
37 
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.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),
528  ExcInternalError());
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 
645  SolverControl &
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()),
678  ExcMessage(
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()),
778  ExcMessage(
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
831  ::Vector<double> & x,
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
882 namespace 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())
Definition: vector.h:109
size_type size() const
iterator begin()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
__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:1614
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1787
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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