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