Reference documentation for deal.II version Git eab6000ff4 2020-05-30 14:50:38 -0400
\(\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);
393  status_test_collection->AddStatusTest(*status_test_abs_tol);
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 
594  SolverCG::AdditionalData::AdditionalData(const bool output_solver_details)
595  : SolverBase::AdditionalData(output_solver_details)
596  {}
597 
598 
599 
601  : SolverBase(cn, data)
602  , additional_data(data)
603  {
604  solver_name = cg;
605  }
606 
607 
608  /* ---------------------- SolverGMRES ------------------------ */
609 
611  const bool output_solver_details,
612  const unsigned int restart_parameter)
613  : SolverBase::AdditionalData(output_solver_details, restart_parameter)
614  {}
615 
616 
617 
619  : SolverBase(cn, data)
620  , additional_data(data)
621  {
622  solver_name = gmres;
623  }
624 
625 
626  /* ---------------------- SolverBicgstab ------------------------ */
627 
629  const bool output_solver_details)
630  : SolverBase::AdditionalData(output_solver_details)
631  {}
632 
633 
634 
636  : SolverBase(cn, data)
637  , additional_data(data)
638  {
640  }
641 
642 
643  /* ---------------------- SolverCGS ------------------------ */
644 
645  SolverCGS::AdditionalData::AdditionalData(const bool output_solver_details)
646  : SolverBase::AdditionalData(output_solver_details)
647  {}
648 
649 
650 
652  : SolverBase(cn, data)
653  , additional_data(data)
654  {
655  solver_name = cgs;
656  }
657 
658 
659  /* ---------------------- SolverTFQMR ------------------------ */
660 
661  SolverTFQMR::AdditionalData::AdditionalData(const bool output_solver_details)
662  : SolverBase::AdditionalData(output_solver_details)
663  {}
664 
665 
666 
668  : SolverBase(cn, data)
669  , additional_data(data)
670  {
671  solver_name = tfqmr;
672  }
673 
674 
675 
676  /* ---------------------- SolverDirect ------------------------ */
677 
678  SolverDirect::AdditionalData::AdditionalData(const bool output_solver_details,
679  const std::string &solver_type)
680  : output_solver_details(output_solver_details)
681  , solver_type(solver_type)
682  {}
683 
684 
685 
687  : solver_control(cn)
689  {}
690 
691 
692 
693  SolverControl &
695  {
696  return solver_control;
697  }
698 
699 
700 
701  void
703  {
704  // We need an Epetra_LinearProblem object to let the Amesos solver know
705  // about the matrix and vectors.
706  linear_problem = std::make_unique<Epetra_LinearProblem>();
707 
708  // Assign the matrix operator to the Epetra_LinearProblem object
709  linear_problem->SetOperator(
710  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()));
711 
712  // Fetch return value of Amesos Solver functions
713  int ierr;
714 
715  // First set whether we want to print the solver information to screen or
716  // not.
717  ConditionalOStream verbose_cout(std::cout,
719 
720  // Next allocate the Amesos solver, this is done in two steps, first we
721  // create a solver Factory and and generate with that the concrete Amesos
722  // solver, if possible.
723  Amesos Factory;
724 
725  AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
726  ExcMessage(
727  std::string("You tried to select the solver type <") +
729  "> but this solver is not supported by Trilinos either "
730  "because it does not exist, or because Trilinos was not "
731  "configured for its use."));
732 
733  solver.reset(
734  Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
735 
736  verbose_cout << "Starting symbolic factorization" << std::endl;
737  ierr = solver->SymbolicFactorization();
738  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
739 
740  verbose_cout << "Starting numeric factorization" << std::endl;
741  ierr = solver->NumericFactorization();
742  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
743  }
744 
745 
746  void
748  {
749  // Assign the empty LHS vector to the Epetra_LinearProblem object
750  linear_problem->SetLHS(&x.trilinos_vector());
751 
752  // Assign the RHS vector to the Epetra_LinearProblem object
753  linear_problem->SetRHS(
754  const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
755 
756  // First set whether we want to print the solver information to screen or
757  // not.
758  ConditionalOStream verbose_cout(std::cout,
760 
761 
762  verbose_cout << "Starting solve" << std::endl;
763 
764  // Fetch return value of Amesos Solver functions
765  int ierr = solver->Solve();
766  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
767 
768  // Finally, force the SolverControl object to report convergence
769  solver_control.check(0, 0);
770  }
771 
772 
773 
774  void
777  const ::LinearAlgebra::distributed::Vector<double> &b)
778  {
779  Epetra_Vector ep_x(View,
780  linear_problem->GetOperator()->OperatorDomainMap(),
781  x.begin());
782  Epetra_Vector ep_b(View,
783  linear_problem->GetOperator()->OperatorRangeMap(),
784  const_cast<double *>(b.begin()));
785 
786  // Assign the empty LHS vector to the Epetra_LinearProblem object
787  linear_problem->SetLHS(&ep_x);
788 
789  // Assign the RHS vector to the Epetra_LinearProblem object
790  linear_problem->SetRHS(&ep_b);
791 
792  // First set whether we want to print the solver information to screen or
793  // not.
794  ConditionalOStream verbose_cout(std::cout,
796 
797  verbose_cout << "Starting solve" << std::endl;
798 
799  // Fetch return value of Amesos Solver functions
800  int ierr = solver->Solve();
801  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
802 
803  // Finally, force the SolverControl object to report convergence
804  solver_control.check(0, 0);
805  }
806 
807 
808 
809  void
811  {
812  // Fetch return value of Amesos Solver functions
813  int ierr;
814 
815  // First set whether we want to print the solver information to screen or
816  // not.
817  ConditionalOStream verbose_cout(std::cout,
819 
820  // Next allocate the Amesos solver, this is done in two steps, first we
821  // create a solver Factory and and generate with that the concrete Amesos
822  // solver, if possible.
823  Amesos Factory;
824 
825  AssertThrow(Factory.Query(additional_data.solver_type.c_str()),
826  ExcMessage(
827  std::string("You tried to select the solver type <") +
829  "> but this solver is not supported by Trilinos either "
830  "because it does not exist, or because Trilinos was not "
831  "configured for its use."));
832 
833  solver.reset(
834  Factory.Create(additional_data.solver_type.c_str(), *linear_problem));
835 
836  verbose_cout << "Starting symbolic factorization" << std::endl;
837  ierr = solver->SymbolicFactorization();
838  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
839 
840  verbose_cout << "Starting numeric factorization" << std::endl;
841  ierr = solver->NumericFactorization();
842  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
843 
844  verbose_cout << "Starting solve" << std::endl;
845  ierr = solver->Solve();
846  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
847 
848  // Finally, let the deal.II SolverControl object know what has
849  // happened. If the solve succeeded, the status of the solver control will
850  // turn into SolverControl::success.
851  solver_control.check(0, 0);
852 
854  AssertThrow(false,
857  }
858 
859 
860  void
862  MPI::Vector & x,
863  const MPI::Vector & b)
864  {
865  // We need an Epetra_LinearProblem object to let the Amesos solver know
866  // about the matrix and vectors.
867  linear_problem = std::make_unique<Epetra_LinearProblem>(
868  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()),
869  &x.trilinos_vector(),
870  const_cast<Epetra_MultiVector *>(&b.trilinos_vector()));
871 
872  do_solve();
873  }
874 
875 
876 
877  void
879  ::Vector<double> & x,
880  const ::Vector<double> &b)
881  {
882  // In case we call the solver with deal.II vectors, we create views of the
883  // vectors in Epetra format.
884  Assert(x.size() == A.n(), ExcDimensionMismatch(x.size(), A.n()));
885  Assert(b.size() == A.m(), ExcDimensionMismatch(b.size(), A.m()));
886  Assert(A.local_range().second == A.m(),
887  ExcMessage("Can only work in serial when using deal.II vectors."));
888  Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
889  Epetra_Vector ep_b(View,
890  A.trilinos_matrix().RangeMap(),
891  const_cast<double *>(b.begin()));
892 
893  // We need an Epetra_LinearProblem object to let the Amesos solver know
894  // about the matrix and vectors.
895  linear_problem = std::make_unique<Epetra_LinearProblem>(
896  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
897 
898  do_solve();
899  }
900 
901 
902 
903  void
905  const SparseMatrix & A,
907  const ::LinearAlgebra::distributed::Vector<double> &b)
908  {
910  A.trilinos_matrix().DomainMap().NumMyElements());
911  AssertDimension(b.local_size(),
912  A.trilinos_matrix().RangeMap().NumMyElements());
913  Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
914  Epetra_Vector ep_b(View,
915  A.trilinos_matrix().RangeMap(),
916  const_cast<double *>(b.begin()));
917 
918  // We need an Epetra_LinearProblem object to let the Amesos solver know
919  // about the matrix and vectors.
920  linear_problem = std::make_unique<Epetra_LinearProblem>(
921  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
922 
923  do_solve();
924  }
925 } // namespace TrilinosWrappers
926 
927 
928 // explicit instantiations
929 // TODO: put these instantiations into generic file
930 namespace TrilinosWrappers
931 {
932  template void
933  SolverBase::do_solve(const PreconditionBase &preconditioner);
934 
935  template void
936  SolverBase::do_solve(const Epetra_Operator &preconditioner);
937 } // namespace TrilinosWrappers
938 
940 
941 #endif // DEAL_II_WITH_PETSC
static ::ExceptionBase & ExcTrilinosError(int arg1)
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_abs_tol
__global__ void reduction(Number *result, const Number *v, const size_type N)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
virtual State check(const unsigned int step, const double check_value)
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
Definition: fe_values.h:1980
SolverBase(SolverControl &cn, const AdditionalData &data=AdditionalData())
std::unique_ptr< Amesos_BaseSolver > solver
const AdditionalData additional_data
const AdditionalData additional_data
AdditionalData(const bool output_solver_details=false, const unsigned int restart_parameter=30)
const AdditionalData additional_data
STL namespace.
void initialize(const SparseMatrix &A)
static ::ExceptionBase & ExcTrilinosError(int arg1)
void solve(MPI::Vector &x, const MPI::Vector &b)
void solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b, const PreconditionBase &preconditioner)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1520
const AdditionalData additional_data
std::unique_ptr< AztecOO_StatusTestMaxIters > status_test_max_steps
static ::ExceptionBase & ExcMessage(std::string arg1)
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
double initial_residual
Stop iteration, goal reached.
#define Assert(cond, exc)
Definition: exceptions.h:1408
AdditionalData(const bool output_solver_details=false)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
std::unique_ptr< Epetra_LinearProblem > linear_problem
SolverControl & control() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:361
const Epetra_CrsMatrix & trilinos_matrix() const
std::unique_ptr< Epetra_LinearProblem > linear_problem
iterator begin()
double last_value() const
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData(const bool output_solver_details=false)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
unsigned int last_step() const
static const char A
const AdditionalData additional_data
enum TrilinosWrappers::SolverBase::SolverName solver_name
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
std::unique_ptr< AztecOO_StatusTestCombo > status_test_collection
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
State last_check() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:360
const AdditionalData additional_data
double current_residual
size_type size() const
std::pair< size_type, size_type > local_range() const
AdditionalData(const bool output_solver_details=false)
SolverControl & control() const
void do_solve(const Preconditioner &preconditioner)
static ::ExceptionBase & ExcNotImplemented()
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
const Epetra_MultiVector & trilinos_vector() const
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_rel_tol
T max(const T &t, const MPI_Comm &mpi_communicator)
Teuchos::RCP< Epetra_Operator > preconditioner
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTrilinosError(int arg1)