Reference documentation for deal.II version 9.2.0
\(\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 
22 
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 
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_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
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
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
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
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.trilinos_matrix().DomainMap(), x.begin());
189  Epetra_Vector ep_b(View,
190  A.trilinos_matrix().RangeMap(),
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
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.
233  A.trilinos_matrix().DomainMap().NumMyElements());
234  AssertDimension(b.local_size(),
235  A.trilinos_matrix().RangeMap().NumMyElements());
236 
237  Epetra_Vector ep_x(View, A.trilinos_matrix().DomainMap(), x.begin());
238  Epetra_Vector ep_b(View,
239  A.trilinos_matrix().RangeMap(),
240  const_cast<double *>(b.begin()));
241 
242  // We need an Epetra_LinearProblem object to let the AztecOO solver know
243  // about the matrix and vectors.
244  linear_problem = std_cxx14::make_unique<Epetra_LinearProblem>(
245  const_cast<Epetra_CrsMatrix *>(&A.trilinos_matrix()), &ep_x, &ep_b);
246 
247  do_solve(preconditioner);
248  }
249 
250 
251 
252  void
255  const ::LinearAlgebra::distributed::Vector<double> &b,
256  const PreconditionBase &preconditioner)
257  {
258  AssertDimension(x.local_size(), A.OperatorDomainMap().NumMyElements());
259  AssertDimension(b.local_size(), A.OperatorRangeMap().NumMyElements());
260 
261  Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.begin());
262  Epetra_Vector ep_b(View,
263  A.OperatorRangeMap(),
264  const_cast<double *>(b.begin()));
265 
266  // We need an Epetra_LinearProblem object to let the AztecOO solver know
267  // about the matrix and vectors.
269  std_cxx14::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
270 
271  do_solve(preconditioner);
272  }
273 
274 
275  namespace internal
276  {
277  namespace
278  {
279  double
280  compute_residual(const Epetra_MultiVector *const residual_vector)
281  {
282  Assert(residual_vector->NumVectors() == 1,
283  ExcMessage("Residual multivector holds more than one vector"));
284  TrilinosScalar res_l2_norm = 0.0;
285  const int ierr = residual_vector->Norm2(&res_l2_norm);
286  AssertThrow(ierr == 0, ExcTrilinosError(ierr));
287  return res_l2_norm;
288  }
289 
290  class TrilinosReductionControl : public AztecOO_StatusTest
291  {
292  public:
293  TrilinosReductionControl(const int max_steps,
294  const double tolerance,
295  const double reduction,
296  const Epetra_LinearProblem &linear_problem);
297 
298  virtual ~TrilinosReductionControl() override = default;
299 
300  virtual bool
301  ResidualVectorRequired() const override
302  {
303  return status_test_collection->ResidualVectorRequired();
304  }
305 
306  virtual AztecOO_StatusType
307  CheckStatus(int CurrentIter,
308  Epetra_MultiVector *CurrentResVector,
309  double CurrentResNormEst,
310  bool SolutionUpdated) override
311  {
312  // Note: CurrentResNormEst is set to -1.0 if no estimate of the
313  // residual value is available
315  (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
316  CurrentResNormEst);
317  if (CurrentIter == 0)
319 
320  return status_test_collection->CheckStatus(CurrentIter,
321  CurrentResVector,
322  CurrentResNormEst,
323  SolutionUpdated);
324  }
325 
326  virtual AztecOO_StatusType
327  GetStatus() const override
328  {
329  return status_test_collection->GetStatus();
330  }
331 
332  virtual std::ostream &
333  Print(std::ostream &stream, int indent = 0) const override
334  {
335  return status_test_collection->Print(stream, indent);
336  }
337 
338  double
339  get_initial_residual() const
340  {
341  return initial_residual;
342  }
343 
344  double
345  get_current_residual() const
346  {
347  return current_residual;
348  }
349 
350  private:
353  std::unique_ptr<AztecOO_StatusTestCombo> status_test_collection;
354  std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
355  std::unique_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
356  std::unique_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
357  };
358 
359 
360  TrilinosReductionControl::TrilinosReductionControl(
361  const int max_steps,
362  const double tolerance,
363  const double reduction,
364  const Epetra_LinearProblem &linear_problem)
365  : initial_residual(std::numeric_limits<double>::max())
366  , current_residual(std::numeric_limits<double>::max())
367  // Consider linear problem converged if any of the collection of
368  // criterion are met
370  std_cxx14::make_unique<AztecOO_StatusTestCombo>(
371  AztecOO_StatusTestCombo::OR))
372  {
373  // Maximum number of iterations
374  Assert(max_steps >= 0, ExcInternalError());
376  std_cxx14::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
378 
379  Assert(linear_problem.GetRHS()->NumVectors() == 1,
380  ExcMessage("RHS multivector holds more than one vector"));
381 
382  // Residual norm is below some absolute value
383  status_test_abs_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
384  *linear_problem.GetOperator(),
385  *(linear_problem.GetLHS()->operator()(0)),
386  *(linear_problem.GetRHS()->operator()(0)),
387  tolerance);
388  status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
389  AztecOO_StatusTestResNorm::TwoNorm);
390  status_test_abs_tol->DefineScaleForm(
391  AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
393 
394  // Residual norm, scaled by some initial value, is below some threshold
395  status_test_rel_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
396  *linear_problem.GetOperator(),
397  *(linear_problem.GetLHS()->operator()(0)),
398  *(linear_problem.GetRHS()->operator()(0)),
399  reduction);
400  status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
401  AztecOO_StatusTestResNorm::TwoNorm);
402  status_test_rel_tol->DefineScaleForm(
403  AztecOO_StatusTestResNorm::NormOfInitRes,
404  AztecOO_StatusTestResNorm::TwoNorm);
406  }
407 
408  } // namespace
409  } // namespace internal
410 
411 
412  template <typename Preconditioner>
413  void
414  SolverBase::do_solve(const Preconditioner &preconditioner)
415  {
416  int ierr;
417 
418  // Next we can allocate the AztecOO solver...
419  solver.SetProblem(*linear_problem);
420 
421  // ... and we can specify the solver to be used.
422  switch (solver_name)
423  {
424  case cg:
425  solver.SetAztecOption(AZ_solver, AZ_cg);
426  break;
427  case cgs:
428  solver.SetAztecOption(AZ_solver, AZ_cgs);
429  break;
430  case gmres:
431  solver.SetAztecOption(AZ_solver, AZ_gmres);
432  solver.SetAztecOption(AZ_kspace,
433  additional_data.gmres_restart_parameter);
434  break;
435  case bicgstab:
436  solver.SetAztecOption(AZ_solver, AZ_bicgstab);
437  break;
438  case tfqmr:
439  solver.SetAztecOption(AZ_solver, AZ_tfqmr);
440  break;
441  default:
442  Assert(false, ExcNotImplemented());
443  }
444 
445  // Set the preconditioner
446  set_preconditioner(solver, preconditioner);
447 
448  // ... set some options, ...
449  solver.SetAztecOption(AZ_output,
450  additional_data.output_solver_details ? AZ_all :
451  AZ_none);
452  solver.SetAztecOption(AZ_conv, AZ_noscaled);
453 
454  // By default, the Trilinos solver chooses convergence criterion based on
455  // the number of iterations made and an absolute tolerance.
456  // This implies that the use of the standard Trilinos convergence test
457  // actually coincides with ::IterationNumberControl because the
458  // solver, unless explicitly told otherwise, will Iterate() until a number
459  // of max_steps() are taken or an absolute tolerance() is attained.
460  // It is therefore suitable for use with both SolverControl or
461  // IterationNumberControl. The final check at the end will determine whether
462  // failure to converge to the defined residual norm constitutes failure
463  // (SolverControl) or is alright (IterationNumberControl).
464  // In the case that the SolverControl wants to perform ReductionControl,
465  // then we have to do a little extra something by prescribing a custom
466  // status test.
467  if (!status_test)
468  {
469  if (const ReductionControl *const reduction_control =
470  dynamic_cast<const ReductionControl *const>(&solver_control))
471  {
472  status_test =
473  std_cxx14::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)
688  , additional_data(data.output_solver_details, data.solver_type)
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_cxx14::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 
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(
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 
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(
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_cxx14::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_cxx14::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_cxx14::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
SolverControl::check
virtual State check(const unsigned int step, const double check_value)
Definition: solver_control.cc:52
TrilinosWrappers::SolverBase::gmres
@ gmres
Definition: trilinos_solver.h:91
TrilinosWrappers::SolverTFQMR::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:552
TrilinosWrappers::MPI::Vector
Definition: trilinos_vector.h:400
TrilinosWrappers::SolverBase::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:342
SolverControl::last_check
State last_check() const
Definition: solver_control.cc:106
TrilinosWrappers::SolverBase::solve
void solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b, const PreconditionBase &preconditioner)
Definition: trilinos_solver.cc:76
SolverControl::NoConvergence
Definition: solver_control.h:96
TrilinosWrappers::SolverDirect::linear_problem
std::unique_ptr< Epetra_LinearProblem > linear_problem
Definition: trilinos_solver.h:716
TrilinosWrappers::SparseMatrix
Definition: trilinos_sparse_matrix.h:515
LinearAlgebra::distributed::Vector< double >
TrilinosWrappers::SolverGMRES::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:472
TrilinosWrappers::SolverTFQMR::SolverTFQMR
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:667
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Vector::begin
iterator begin()
TrilinosWrappers::SolverDirect::solver
std::unique_ptr< Amesos_BaseSolver > solver
Definition: trilinos_solver.h:722
TrilinosWrappers::SolverTFQMR::AdditionalData
Definition: trilinos_solver.h:530
TrilinosWrappers::SolverDirect::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:727
TrilinosWrappers
Definition: types.h:161
TrilinosWrappers::SolverDirect::AdditionalData::output_solver_details
bool output_solver_details
Definition: trilinos_solver.h:589
TrilinosWrappers::SolverCGS::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false)
Definition: trilinos_solver.cc:645
trilinos_sparse_matrix.h
TrilinosWrappers::SolverBase::SolverBase
SolverBase(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:49
initial_residual
double initial_residual
Definition: trilinos_solver.cc:351
LinearAlgebra::CUDAWrappers::kernel::reduction
__global__ void reduction(Number *result, const Number *v, const size_type N)
status_test_max_steps
std::unique_ptr< AztecOO_StatusTestMaxIters > status_test_max_steps
Definition: trilinos_solver.cc:354
SolverControl::last_value
double last_value() const
Definition: solver_control.cc:120
ReductionControl
Definition: solver_control.h:427
Vector::size
size_type size() const
TrilinosWrappers::SolverBase::solver_control
SolverControl & solver_control
Definition: trilinos_solver.h:302
TrilinosWrappers::SolverCGS::SolverCGS
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:651
TrilinosWrappers::SolverBase::tfqmr
@ tfqmr
Definition: trilinos_solver.h:99
TrilinosWrappers::SolverDirect::control
SolverControl & control() const
Definition: trilinos_solver.cc:694
TrilinosWrappers::SolverBase::bicgstab
@ bicgstab
Definition: trilinos_solver.h:95
SolverControl::last_step
unsigned int last_step() const
Definition: solver_control.cc:127
TrilinosWrappers::SolverBase::SolverName
SolverName
Definition: trilinos_solver.h:78
TrilinosWrappers::SolverBase::solver_name
enum TrilinosWrappers::SolverBase::SolverName solver_name
TrilinosWrappers::SolverDirect::initialize
void initialize(const SparseMatrix &A)
Definition: trilinos_solver.cc:702
TrilinosWrappers::SolverCG::SolverCG
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:600
conditional_ostream.h
TrilinosWrappers::SolverCG::AdditionalData
Definition: trilinos_solver.h:372
TrilinosWrappers::SolverGMRES::SolverGMRES
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:618
TrilinosWrappers::SolverDirect::AdditionalData
Definition: trilinos_solver.h:577
TrilinosWrappers::SolverBase::AdditionalData
Definition: trilinos_solver.h:106
TrilinosWrappers::SolverBase::linear_problem
std::unique_ptr< Epetra_LinearProblem > linear_problem
Definition: trilinos_solver.h:325
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
TrilinosWrappers::PreconditionBase::preconditioner
Teuchos::RCP< Epetra_Operator > preconditioner
Definition: trilinos_precondition.h:239
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
TrilinosWrappers::SolverTFQMR::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false)
Definition: trilinos_solver.cc:661
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
TrilinosWrappers::SolverBicgstab::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false)
Definition: trilinos_solver.cc:628
current_residual
double current_residual
Definition: trilinos_solver.cc:352
LinearAlgebra::distributed::Vector::local_size
size_type local_size() const
TrilinosWrappers::SolverBicgstab::AdditionalData
Definition: trilinos_solver.h:490
status_test_collection
std::unique_ptr< AztecOO_StatusTestCombo > status_test_collection
Definition: trilinos_solver.cc:353
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
TrilinosWrappers::MPI::Vector::trilinos_vector
const Epetra_MultiVector & trilinos_vector() const
TrilinosWrappers::SolverGMRES::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false, const unsigned int restart_parameter=30)
Definition: trilinos_solver.cc:610
TrilinosWrappers::SolverDirect::do_solve
void do_solve()
Definition: trilinos_solver.cc:810
TrilinosWrappers::PreconditionBase
Definition: trilinos_precondition.h:82
TrilinosWrappers::SolverDirect::solve
void solve(MPI::Vector &x, const MPI::Vector &b)
Definition: trilinos_solver.cc:747
TrilinosWrappers::SolverBase::control
SolverControl & control() const
Definition: trilinos_solver.cc:68
TrilinosWrappers::SolverBase::do_solve
void do_solve(const Preconditioner &preconditioner)
Definition: trilinos_solver.cc:414
trilinos_precondition.h
TrilinosWrappers::SolverDirect::AdditionalData::solver_type
std::string solver_type
Definition: trilinos_solver.h:609
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
TrilinosWrappers::SolverCGS::AdditionalData
Definition: trilinos_solver.h:410
TrilinosWrappers::SolverDirect::SolverDirect
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:686
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
TrilinosWrappers::SolverBase::cgs
@ cgs
Definition: trilinos_solver.h:87
LACExceptions::ExcTrilinosError
static ::ExceptionBase & ExcTrilinosError(int arg1)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
TrilinosWrappers::SolverGMRES::AdditionalData
Definition: trilinos_solver.h:448
TrilinosWrappers::SolverBase::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
Definition: trilinos_solver.cc:40
ConditionalOStream
Definition: conditional_ostream.h:81
TrilinosWrappers::SolverBicgstab::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:512
TrilinosWrappers::SolverCG::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:393
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SolverControl::success
@ success
Stop iteration, goal reached.
Definition: solver_control.h:79
TrilinosWrappers::SolverBase::cg
@ cg
Definition: trilinos_solver.h:83
internal
Definition: aligned_vector.h:369
SolverControl
Definition: solver_control.h:67
memory.h
trilinos_solver.h
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
status_test_abs_tol
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_abs_tol
Definition: trilinos_solver.cc:355
TrilinosWrappers::SolverBase
Definition: trilinos_solver.h:68
trilinos_vector.h
Epetra_Operator
std_cxx14
Definition: c++.h:129
Vector< double >
TrilinosWrappers::SolverDirect::solver_control
SolverControl & solver_control
Definition: trilinos_solver.h:709
TrilinosWrappers::SolverCGS::additional_data
const AdditionalData additional_data
Definition: trilinos_solver.h:431
TrilinosWrappers::SolverDirect::ExcTrilinosError
static ::ExceptionBase & ExcTrilinosError(int arg1)
Factory
LinearAlgebra::distributed::Vector::begin
iterator begin()
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
FEValuesViews::View
typename ::internal::FEValuesViews::ViewType< dim, spacedim, Extractor >::type View
Definition: fe_values.h:1980
TrilinosWrappers::SolverBicgstab::SolverBicgstab
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
Definition: trilinos_solver.cc:635
status_test_rel_tol
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_rel_tol
Definition: trilinos_solver.cc:356
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
TrilinosWrappers::SolverDirect::AdditionalData::AdditionalData
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
Definition: trilinos_solver.cc:678