]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Cleanup of SolverGMRES and SolverFGMRES implementations
authorMartin Kronbichler <martin.kronbichler@rub.de>
Mon, 11 Mar 2024 14:50:56 +0000 (15:50 +0100)
committerMartin Kronbichler <martin.kronbichler@rub.de>
Thu, 14 Mar 2024 09:38:13 +0000 (10:38 +0100)
26 files changed:
examples/step-12/step-12.cc
examples/step-22/doc/results.dox
include/deal.II/lac/solver_gmres.h
tests/examples/step-56.output
tests/lac/gmres_reorthogonalize_01.cc
tests/lac/gmres_reorthogonalize_02.cc
tests/lac/gmres_reorthogonalize_03.cc
tests/lac/gmres_reorthogonalize_04.cc
tests/lac/gmres_reorthogonalize_05.cc
tests/lac/solver.cc
tests/lac/solver.output
tests/lac/solver_gmres_01.cc
tests/lac/solver_gmres_02.cc
tests/lac/solver_monitor.cc
tests/lac/solver_selector_01.output
tests/lac/solver_selector_03.output
tests/lac/solver_signals.cc
tests/lac/solver_signals.with_lapack=true.output
tests/matrix_free/stokes_computation.cc
tests/matrix_free/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output
tests/multigrid-global-coarsening/stokes_computation.with_lapack=true.with_p4est=true.mpirun=4.output
tests/multigrid/step-39-03.output
tests/non_matching/step-70.with_p4est=true.with_petsc=true.mpirun=2.output
tests/trilinos/deal_solver_03.cc
tests/trilinos/deal_solver_03.output
tests/trilinos_tpetra/deal_solver_03.cc

index 511acf5dabf973aae7a33945de25bf2b0cd8364f..6c7de2cbde549d1df22b98010dcd7ca0f5eca365 100644 (file)
@@ -481,7 +481,7 @@ namespace Step12
     SolverControl solver_control(1000, 1e-6 * right_hand_side.l2_norm());
 
     SolverGMRES<Vector<double>>::AdditionalData additional_data;
-    additional_data.max_n_tmp_vectors = 100;
+    additional_data.max_basis_size = 100;
     SolverGMRES<Vector<double>> solver(solver_control, additional_data);
 
     // Here we create the preconditioner,
index 32f7be0f72659575ca07a3eb2f4442438de00634..30473dc2257e26373877292f45166bd19805f7c0 100644 (file)
@@ -544,7 +544,7 @@ $k=100$ temporary vectors:
                                     1e-6*system_rhs.l2_norm());
       GrowingVectorMemory<BlockVector<double> > vector_memory;
       SolverGMRES<BlockVector<double> >::AdditionalData gmres_data;
-      gmres_data.max_n_tmp_vectors = 100;
+      gmres_data.max_basis_size = 100;
 
       SolverGMRES<BlockVector<double> > gmres(solver_control, vector_memory,
                                               gmres_data);
index f722f144b4c53ca80222b02dbb161768e9e6f9df..78cf154529159344a5a19704efb0eb6aece5d707 100644 (file)
@@ -25,7 +25,6 @@
 
 #include <deal.II/lac/block_vector_base.h>
 #include <deal.II/lac/full_matrix.h>
-#include <deal.II/lac/householder.h>
 #include <deal.II/lac/lapack_full_matrix.h>
 #include <deal.II/lac/orthogonalization.h>
 #include <deal.II/lac/solver.h>
@@ -129,11 +128,14 @@ namespace internal
  * Implementation of the Restarted Preconditioned Direct Generalized Minimal
  * Residual Method. The stopping criterion is the norm of the residual.
  *
- * The AdditionalData structure contains the number of temporary vectors used.
- * The size of the Arnoldi basis is this number minus three. Additionally, it
- * allows you to choose between right or left preconditioning. The default is
- * left preconditioning. Finally it includes a flag indicating whether or not
- * the default residual is used as stopping criterion.
+ * The AdditionalData structure contains the size of the Arnoldi basis used
+ * for orthogonalization. It is related to the number of temporary vectors
+ * used, which is the basis size plus two. Additionally, it allows you to
+ * choose between right or left preconditioning. The default is left
+ * preconditioning. Furthermore, it includes a flag indicating whether or not
+ * the default residual is used as stopping criterion and an option for the
+ * orthogonalization algorithm, see LinearAlgebra::OrthogonalizationStrategy
+ * for available options.
  *
  *
  * <h3>Left versus right preconditioning</h3>
@@ -142,7 +144,7 @@ namespace internal
  * preconditioning. As expected, this switches between solving for the systems
  * <i>P<sup>-1</sup>A</i> and <i>AP<sup>-1</sup></i>, respectively.
  *
- * A second consequence is the type of residual which is used to measure
+ * A second consequence is the type of residual, which is used to measure
  * convergence. With left preconditioning, this is the <b>preconditioned</b>
  * residual, while with right preconditioning, it is the residual of the
  * unpreconditioned system.
@@ -156,13 +158,13 @@ namespace internal
  *
  * <h3>The size of the Arnoldi basis</h3>
  *
- * The maximal basis size is controlled by AdditionalData::max_n_tmp_vectors,
- * and it is this number minus 2. If the number of iteration steps exceeds
- * this number, all basis vectors are discarded and the iteration starts anew
- * from the approximation obtained so far.
+ * The maximal basis size is controlled by AdditionalData::max_basis_size. If
+ * the number of iteration steps exceeds this number, all basis vectors are
+ * discarded and the iteration starts anew from the approximation obtained so
+ * far.
  *
- * Note that the minimizing property of GMRes only pertains to the Krylov
- * space spanned by the Arnoldi basis. Therefore, restarted GMRes is
+ * Note that the minimizing property of GMRES only pertains to the Krylov
+ * space spanned by the Arnoldi basis. Therefore, restarted GMRES is
  * <b>not</b> minimizing anymore. The choice of the basis length is a
  * trade-off between memory consumption and convergence speed, since a longer
  * basis means minimization over a larger space.
@@ -199,14 +201,15 @@ public:
   struct AdditionalData
   {
     /**
-     * Constructor. By default, set the number of temporary vectors to 30,
-     * i.e. do a restart every 28 iterations. Also set preconditioning from
-     * left, the residual of the stopping criterion to the default residual,
-     * and re-orthogonalization only if necessary. Also, the batched mode with
-     * reduced functionality to track information is disabled by default.
+     * Constructor. By default, set the size of the Arnoldi basis.  Also set
+     * preconditioning from left, the residual of the stopping criterion to
+     * the default residual, and re-orthogonalization only if necessary. Also,
+     * the batched mode with reduced functionality to track information is
+     * disabled by default. Finally, the default orthogonalization algorithm
+     * is the modified Gram-Schmidt method.
      */
     explicit AdditionalData(
-      const unsigned int max_n_tmp_vectors          = 30,
+      const unsigned int max_basis_size             = 30,
       const bool         right_preconditioning      = false,
       const bool         use_default_residual       = true,
       const bool         force_re_orthogonalization = false,
@@ -216,13 +219,27 @@ public:
           LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt);
 
     /**
-     * Maximum number of temporary vectors. This parameter controls the size
-     * of the Arnoldi basis, which for historical reasons is
-     * #max_n_tmp_vectors-2. SolverGMRES assumes that there are at least three
-     * temporary vectors, so this value must be greater than or equal to three.
+     * Maximum number of temporary vectors. Together with #max_basis_size,
+     * this parameter controls the size of the Arnoldi basis, which for
+     * historical reasons is #max_n_tmp_vectors-2. SolverGMRES assumes that
+     * there are at least three temporary vectors, so this value must be
+     * greater than or equal to three. If both this variable and
+     * #max_basis_size are set to a non-zero value, the solver uses the latter
+     * variable.
+     *
+     * @deprecated Use #max_basis_size instead.
      */
     unsigned int max_n_tmp_vectors;
 
+    /**
+     * Maximum size of the Arnoldi basis. SolverGMRES assumes that there is at
+     * least one vector in the Arnoldi basis, so this value must be greater
+     * than or equal to one. Note that whenever this variable is set to a
+     * non-zero, including the value set by the default constructor, this
+     * variable takes precedence.
+     */
+    unsigned int max_basis_size;
+
     /**
      * Flag for right preconditioning.
      *
@@ -419,17 +436,6 @@ protected:
   virtual double
   criterion();
 
-  /**
-   * Transformation of an upper Hessenberg matrix into tridiagonal structure
-   * by givens rotation of the last column
-   */
-  void
-  givens_rotation(Vector<double> &h,
-                  Vector<double> &b,
-                  Vector<double> &ci,
-                  Vector<double> &si,
-                  int             col) const;
-
   /**
    * Estimates the eigenvalues from the Hessenberg matrix, H_orig, generated
    * during the inner iterations. Uses these estimate to compute the condition
@@ -454,17 +460,12 @@ protected:
   /**
    * Auxiliary vector for orthogonalization
    */
-  Vector<double> gamma;
-
-  /**
-   * Auxiliary vector for orthogonalization
-   */
-  Vector<double> ci;
+  Vector<double> projected_rhs;
 
   /**
    * Auxiliary vector for orthogonalization
    */
-  Vector<double> si;
+  std::vector<std::pair<double, double>> givens_rotations;
 
   /**
    * Auxiliary vector for orthogonalization
@@ -572,6 +573,51 @@ private:
 
 
 #ifndef DOXYGEN
+
+template <typename VectorType>
+inline SolverGMRES<VectorType>::AdditionalData::AdditionalData(
+  const unsigned int                             max_basis_size,
+  const bool                                     right_preconditioning,
+  const bool                                     use_default_residual,
+  const bool                                     force_re_orthogonalization,
+  const bool                                     batched_mode,
+  const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
+  : max_n_tmp_vectors(0)
+  , max_basis_size(max_basis_size)
+  , right_preconditioning(right_preconditioning)
+  , use_default_residual(use_default_residual)
+  , force_re_orthogonalization(force_re_orthogonalization)
+  , batched_mode(batched_mode)
+  , orthogonalization_strategy(orthogonalization_strategy)
+{
+  Assert(max_basis_size >= 1,
+         ExcMessage("SolverGMRES needs at least one vector in the "
+                    "Arnoldi basis."));
+}
+
+
+
+template <typename VectorType>
+SolverGMRES<VectorType>::SolverGMRES(SolverControl            &cn,
+                                     VectorMemory<VectorType> &mem,
+                                     const AdditionalData     &data)
+  : SolverBase<VectorType>(cn, mem)
+  , additional_data(data)
+  , solver_control(cn)
+{}
+
+
+
+template <typename VectorType>
+SolverGMRES<VectorType>::SolverGMRES(SolverControl        &cn,
+                                     const AdditionalData &data)
+  : SolverBase<VectorType>(cn)
+  , additional_data(data)
+  , solver_control(cn)
+{}
+
+
+
 namespace internal
 {
   namespace SolverGMRESImplementation
@@ -622,111 +668,6 @@ namespace internal
 
 
 
-    // A comparator for better printing eigenvalues
-    inline bool
-    complex_less_pred(const std::complex<double> &x,
-                      const std::complex<double> &y)
-    {
-      return x.real() < y.real() ||
-             (x.real() == y.real() && x.imag() < y.imag());
-    }
-
-    // A function to solve the (upper) triangular system after Givens
-    // rotations on a matrix that has possibly unused rows and columns
-    inline void
-    solve_triangular(const unsigned int        dim,
-                     const FullMatrix<double> &H,
-                     const Vector<double>     &rhs,
-                     Vector<double>           &solution)
-    {
-      for (int i = dim - 1; i >= 0; --i)
-        {
-          double s = rhs(i);
-          for (unsigned int j = i + 1; j < dim; ++j)
-            s -= solution(j) * H(i, j);
-          solution(i) = s / H(i, i);
-          AssertIsFinite(solution(i));
-        }
-    }
-  } // namespace SolverGMRESImplementation
-} // namespace internal
-
-
-
-template <typename VectorType>
-inline SolverGMRES<VectorType>::AdditionalData::AdditionalData(
-  const unsigned int                             max_n_tmp_vectors,
-  const bool                                     right_preconditioning,
-  const bool                                     use_default_residual,
-  const bool                                     force_re_orthogonalization,
-  const bool                                     batched_mode,
-  const LinearAlgebra::OrthogonalizationStrategy orthogonalization_strategy)
-  : max_n_tmp_vectors(max_n_tmp_vectors)
-  , right_preconditioning(right_preconditioning)
-  , use_default_residual(use_default_residual)
-  , force_re_orthogonalization(force_re_orthogonalization)
-  , batched_mode(batched_mode)
-  , orthogonalization_strategy(orthogonalization_strategy)
-{
-  Assert(3 <= max_n_tmp_vectors,
-         ExcMessage("SolverGMRES needs at least three "
-                    "temporary vectors."));
-}
-
-
-
-template <typename VectorType>
-SolverGMRES<VectorType>::SolverGMRES(SolverControl            &cn,
-                                     VectorMemory<VectorType> &mem,
-                                     const AdditionalData     &data)
-  : SolverBase<VectorType>(cn, mem)
-  , additional_data(data)
-  , solver_control(cn)
-{}
-
-
-
-template <typename VectorType>
-SolverGMRES<VectorType>::SolverGMRES(SolverControl        &cn,
-                                     const AdditionalData &data)
-  : SolverBase<VectorType>(cn)
-  , additional_data(data)
-  , solver_control(cn)
-{}
-
-
-
-template <typename VectorType>
-inline void
-SolverGMRES<VectorType>::givens_rotation(Vector<double> &h,
-                                         Vector<double> &b,
-                                         Vector<double> &ci,
-                                         Vector<double> &si,
-                                         int             col) const
-{
-  for (int i = 0; i < col; ++i)
-    {
-      const double s     = si(i);
-      const double c     = ci(i);
-      const double dummy = h(i);
-      h(i)               = c * dummy + s * h(i + 1);
-      h(i + 1)           = -s * dummy + c * h(i + 1);
-    };
-
-  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
-  si(col)        = h(col + 1) * r;
-  ci(col)        = h(col) * r;
-  h(col)         = ci(col) * h(col) + si(col) * h(col + 1);
-  b(col + 1)     = -si(col) * b(col);
-  b(col) *= ci(col);
-}
-
-
-
-namespace internal
-{
-  namespace SolverGMRESImplementation
-  {
     template <typename VectorType, typename Enable = void>
     struct is_dealii_compatible_distributed_vector;
 
@@ -1225,6 +1166,65 @@ namespace internal
 
       return 0.0;
     }
+
+    // A comparator for better printing eigenvalues
+    inline bool
+    complex_less_pred(const std::complex<double> &x,
+                      const std::complex<double> &y)
+    {
+      return x.real() < y.real() ||
+             (x.real() == y.real() && x.imag() < y.imag());
+    }
+
+
+
+    // A function to compute the Givens rotation for the QR factorization of
+    // the Hessenberg matrix involved in the Arnoldi process, transforming it
+    // into an upper triangular matrix.
+    inline void
+    givens_rotation(Vector<double>                         &h,
+                    Vector<double>                         &b,
+                    std::vector<std::pair<double, double>> &rotations,
+                    const int                               col)
+    {
+      for (int i = 0; i < col; ++i)
+        {
+          const double c   = rotations[i].first;
+          const double s   = rotations[i].second;
+          const double tmp = h(i);
+          h(i)             = c * tmp + s * h(i + 1);
+          h(i + 1)         = -s * tmp + c * h(i + 1);
+        }
+
+      const double r =
+        1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
+      rotations[col].second = h(col + 1) * r;
+      rotations[col].first  = h(col) * r;
+      h(col) =
+        rotations[col].first * h(col) + rotations[col].second * h(col + 1);
+      b(col + 1) = -rotations[col].second * b(col);
+      b(col) *= rotations[col].first;
+    }
+
+
+
+    // A function to solve the (upper) triangular system after Givens
+    // rotations on a matrix that has possibly unused rows and columns
+    inline void
+    solve_triangular(const unsigned int        dim,
+                     const FullMatrix<double> &H,
+                     const Vector<double>     &rhs,
+                     Vector<double>           &solution)
+    {
+      for (int i = dim - 1; i >= 0; --i)
+        {
+          double s = rhs(i);
+          for (unsigned int j = i + 1; j < dim; ++j)
+            s -= solution(j) * H(i, j);
+          solution(i) = s / H(i, i);
+          AssertIsFinite(solution(i));
+        }
+    }
   } // namespace SolverGMRESImplementation
 } // namespace internal
 
@@ -1290,22 +1290,20 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
                                const VectorType         &b,
                                const PreconditionerType &preconditioner)
 {
-  // TODO:[?] Check, why there are two different start residuals.
-  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
-  // size
-
   std::unique_ptr<LogStream::Prefix> prefix;
   if (!additional_data.batched_mode)
     prefix = std::make_unique<LogStream::Prefix>("GMRES");
 
   // extra call to std::max to placate static analyzers: coverity rightfully
   // complains that data.max_n_tmp_vectors - 2 may overflow
-  const unsigned int n_tmp_vectors =
-    std::max(additional_data.max_n_tmp_vectors, 3u);
+  const unsigned int basis_size =
+    (additional_data.max_basis_size > 0 ?
+       additional_data.max_basis_size :
+       std::max(additional_data.max_n_tmp_vectors, 3u) - 2);
 
   // Generate an object where basis vectors are stored.
-  internal::SolverGMRESImplementation::TmpVectors<VectorType> tmp_vectors(
-    n_tmp_vectors, this->memory);
+  internal::SolverGMRESImplementation::TmpVectors<VectorType> basis_vectors(
+    basis_size + 2, this->memory);
 
   // number of the present iteration; this
   // number is not reset to zero upon a
@@ -1322,21 +1320,18 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
   // applying Givens rotations)
   FullMatrix<double> H_orig;
   if (do_eigenvalues)
-    H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
+    H_orig.reinit(basis_size + 1, basis_size);
 
   // matrix used for the orthogonalization process later
-  H.reinit(n_tmp_vectors, n_tmp_vectors - 1, /* omit_initialization */ true);
+  H.reinit(basis_size + 1, basis_size, /* omit_initialization */ true);
 
   // some additional vectors, also used in the orthogonalization
-  gamma.reinit(n_tmp_vectors);
-  ci.reinit(n_tmp_vectors - 1);
-  si.reinit(n_tmp_vectors - 1);
-  h.reinit(n_tmp_vectors - 1);
-
-  unsigned int dim = 0;
+  projected_rhs.reinit(basis_size + 1);
+  givens_rotations.resize(basis_size);
+  h.reinit(basis_size + 1);
 
   SolverControl::State iteration_state = SolverControl::iterate;
-  double               last_res        = std::numeric_limits<double>::lowest();
+  double               res             = std::numeric_limits<double>::lowest();
 
   // switch to determine whether we want a left or a right preconditioner. at
   // present, left is default, but both ways are implemented
@@ -1347,15 +1342,14 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
   // residual as stopping criterion.
   const bool use_default_residual = additional_data.use_default_residual;
 
-  // define two aliases
-  VectorType &v = tmp_vectors(0, x);
-  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
+  // define an alias
+  VectorType &p = basis_vectors(basis_size + 1, x);
 
   // Following vectors are needed when we are not using the default residuals
   // as stopping criterion
   typename VectorMemory<VectorType>::Pointer r;
   typename VectorMemory<VectorType>::Pointer x_;
-  std::unique_ptr<dealii::Vector<double>>    gamma_;
+  std::unique_ptr<dealii::Vector<double>>    gamma;
   if (!use_default_residual)
     {
       r  = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
@@ -1363,7 +1357,7 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
       r->reinit(x);
       x_->reinit(x);
 
-      gamma_ = std::make_unique<dealii::Vector<double>>(gamma.size());
+      gamma = std::make_unique<dealii::Vector<double>>(projected_rhs.size());
     }
 
   bool re_orthogonalize = additional_data.force_re_orthogonalization;
@@ -1374,42 +1368,45 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
   // restart
   do
     {
-      double rho = 0.0;
+      VectorType &v      = basis_vectors(0, x);
+      double      norm_v = 0.;
 
       if (left_precondition)
         {
           A.vmult(p, x);
           p.sadd(-1., 1., b);
           preconditioner.vmult(v, p);
-          rho = v.l2_norm();
+          norm_v = v.l2_norm();
         }
       else
         {
           A.vmult(v, x);
-          rho = dealii::internal::SolverGMRESImplementation::sadd_and_norm(v,
-                                                                           -1,
-                                                                           b,
-                                                                           1.0);
+          norm_v = dealii::internal::SolverGMRESImplementation::sadd_and_norm(
+            v, -1, b, 1.0);
         }
 
+      projected_rhs(0) = norm_v;
+      if (norm_v != 0)
+        v /= norm_v;
+
       // check the residual here as well since it may be that we got the exact
       // (or an almost exact) solution vector at the outset. if we wouldn't
       // check here, the next scaling operation would produce garbage
       if (use_default_residual)
         {
-          last_res = rho;
+          res = norm_v;
           if (additional_data.batched_mode)
-            iteration_state = solver_control.check(accumulated_iterations, rho);
+            iteration_state = solver_control.check(accumulated_iterations, res);
           else
             iteration_state =
-              this->iteration_status(accumulated_iterations, rho, x);
+              this->iteration_status(accumulated_iterations, res, x);
 
           if (iteration_state != SolverControl::iterate)
             break;
         }
       else
         {
-          deallog << "default_res=" << rho << std::endl;
+          deallog << "default_res=" << norm_v << std::endl;
 
           if (left_precondition)
             {
@@ -1419,10 +1416,9 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
           else
             preconditioner.vmult(*r, v);
 
-          double res = r->l2_norm();
-          last_res   = res;
+          res = r->l2_norm();
           if (additional_data.batched_mode)
-            iteration_state = solver_control.check(accumulated_iterations, rho);
+            iteration_state = solver_control.check(accumulated_iterations, res);
           else
             iteration_state =
               this->iteration_status(accumulated_iterations, res, x);
@@ -1431,98 +1427,91 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
             break;
         }
 
-      gamma(0) = rho;
-
-      v *= 1. / rho;
-
-      // inner iteration doing at most as many steps as there are temporary
-      // vectors. the number of steps actually been done is propagated outside
-      // through the @p dim variable
-      for (unsigned int inner_iteration = 0;
-           ((inner_iteration < n_tmp_vectors - 2) &&
-            (iteration_state == SolverControl::iterate));
+      // inner iteration doing at most as many steps as the size of the
+      // Arnoldi basis
+      unsigned int inner_iteration = 0;
+      for (; (inner_iteration < basis_size &&
+              iteration_state == SolverControl::iterate);
            ++inner_iteration)
         {
           ++accumulated_iterations;
           // yet another alias
-          VectorType &vv = tmp_vectors(inner_iteration + 1, x);
+          VectorType &vv = basis_vectors(inner_iteration + 1, x);
 
           if (left_precondition)
             {
-              A.vmult(p, tmp_vectors[inner_iteration]);
+              A.vmult(p, basis_vectors[inner_iteration]);
               preconditioner.vmult(vv, p);
             }
           else
             {
-              preconditioner.vmult(p, tmp_vectors[inner_iteration]);
+              preconditioner.vmult(p, basis_vectors[inner_iteration]);
               A.vmult(vv, p);
             }
 
-          dim = inner_iteration + 1;
-
-          const double s =
-            internal::SolverGMRESImplementation::iterated_gram_schmidt(
-              additional_data.orthogonalization_strategy,
-              tmp_vectors,
-              dim,
-              accumulated_iterations,
-              vv,
-              h,
-              re_orthogonalize,
-              re_orthogonalize_signal);
-          h(inner_iteration + 1) = s;
-
-          // s=0 is a lucky breakdown, the solver will reach convergence,
-          // but we must not divide by zero here.
-          if (s != 0)
-            vv *= 1. / s;
+          norm_v = internal::SolverGMRESImplementation::iterated_gram_schmidt(
+            additional_data.orthogonalization_strategy,
+            basis_vectors,
+            inner_iteration + 1,
+            accumulated_iterations,
+            vv,
+            h,
+            re_orthogonalize,
+            re_orthogonalize_signal);
+
+          // norm_v = 0 is a lucky breakdown, the solver will reach
+          // convergence, but we must not divide by zero here.
+          if (norm_v != 0)
+            vv /= norm_v;
+
+          h(inner_iteration + 1) = norm_v;
 
           // for eigenvalues, get the resulting coefficients from the
           // orthogonalization process
           if (do_eigenvalues)
-            for (unsigned int i = 0; i < dim + 1; ++i)
+            for (unsigned int i = 0; i < inner_iteration + 2; ++i)
               H_orig(i, inner_iteration) = h(i);
 
           //  Transformation into tridiagonal structure
-          givens_rotation(h, gamma, ci, si, inner_iteration);
+          internal::SolverGMRESImplementation::givens_rotation(h,
+                                                               projected_rhs,
+                                                               givens_rotations,
+                                                               inner_iteration);
 
           //  append vector on matrix
-          for (unsigned int i = 0; i < dim; ++i)
+          for (unsigned int i = 0; i < inner_iteration + 1; ++i)
             H(i, inner_iteration) = h(i);
 
           //  default residual
-          rho = std::fabs(gamma(dim));
+          res = std::fabs(projected_rhs(inner_iteration + 1));
 
           if (use_default_residual)
             {
-              last_res = rho;
               if (additional_data.batched_mode)
                 iteration_state =
-                  solver_control.check(accumulated_iterations, rho);
+                  solver_control.check(accumulated_iterations, res);
               else
                 iteration_state =
-                  this->iteration_status(accumulated_iterations, rho, x);
+                  this->iteration_status(accumulated_iterations, res, x);
             }
           else
             {
               if (!additional_data.batched_mode)
-                deallog << "default_res=" << rho << std::endl;
+                deallog << "default_res=" << res << std::endl;
 
-              *x_     = x;
-              *gamma_ = gamma;
-              internal::SolverGMRESImplementation::solve_triangular(dim,
-                                                                    H,
-                                                                    *gamma_,
-                                                                    h);
+              *x_    = x;
+              *gamma = projected_rhs;
+              internal::SolverGMRESImplementation::solve_triangular(
+                inner_iteration + 1, H, *gamma, h);
 
               if (left_precondition)
-                for (unsigned int i = 0; i < dim; ++i)
-                  x_->add(h(i), tmp_vectors[i]);
+                for (unsigned int i = 0; i < inner_iteration + 1; ++i)
+                  x_->add(h(i), basis_vectors[i]);
               else
                 {
                   p = 0.;
-                  for (unsigned int i = 0; i < dim; ++i)
-                    p.add(h(i), tmp_vectors[i]);
+                  for (unsigned int i = 0; i < inner_iteration + 1; ++i)
+                    p.add(h(i), basis_vectors[i]);
                   preconditioner.vmult(*r, p);
                   x_->add(1., *r);
                 };
@@ -1531,69 +1520,69 @@ SolverGMRES<VectorType>::solve(const MatrixType         &A,
               // Now *r contains the unpreconditioned residual!!
               if (left_precondition)
                 {
-                  const double res = r->l2_norm();
-                  last_res         = res;
-
+                  res = r->l2_norm();
                   iteration_state =
                     this->iteration_status(accumulated_iterations, res, x);
                 }
               else
                 {
                   preconditioner.vmult(*x_, *r);
-                  const double preconditioned_res = x_->l2_norm();
-                  last_res                        = preconditioned_res;
+                  res = x_->l2_norm();
 
                   if (additional_data.batched_mode)
                     iteration_state =
-                      solver_control.check(accumulated_iterations, rho);
+                      solver_control.check(accumulated_iterations, res);
                   else
                     iteration_state =
-                      this->iteration_status(accumulated_iterations,
-                                             preconditioned_res,
-                                             x);
+                      this->iteration_status(accumulated_iterations, res, x);
                 }
             }
         }
 
       // end of inner iteration. now calculate the solution from the temporary
       // vectors
-      internal::SolverGMRESImplementation::solve_triangular(dim, H, gamma, h);
+      internal::SolverGMRESImplementation::solve_triangular(inner_iteration,
+                                                            H,
+                                                            projected_rhs,
+                                                            h);
 
       if (do_eigenvalues)
         compute_eigs_and_cond(H_orig,
-                              dim,
+                              inner_iteration,
                               all_eigenvalues_signal,
                               all_hessenberg_signal,
                               condition_number_signal);
 
       if (left_precondition)
         dealii::internal::SolverGMRESImplementation::add(
-          x, dim, h, tmp_vectors, false);
+          x, inner_iteration, h, basis_vectors, false);
       else
         {
           dealii::internal::SolverGMRESImplementation::add(
-            p, dim, h, tmp_vectors, true);
+            p, inner_iteration, h, basis_vectors, true);
           preconditioner.vmult(v, p);
           x.add(1., v);
         }
+
+      // in the last round, print the eigenvalues from the last Arnoldi step
+      if (iteration_state != SolverControl::iterate && do_eigenvalues)
+        compute_eigs_and_cond(H_orig,
+                              inner_iteration,
+                              eigenvalues_signal,
+                              hessenberg_signal,
+                              condition_number_signal);
+
       // end of outer iteration. restart if no convergence and the number of
       // iterations is not exceeded
     }
   while (iteration_state == SolverControl::iterate);
 
-  if (do_eigenvalues)
-    compute_eigs_and_cond(H_orig,
-                          dim,
-                          eigenvalues_signal,
-                          hessenberg_signal,
-                          condition_number_signal);
-
   if (!additional_data.batched_mode && !krylov_space_signal.empty())
-    krylov_space_signal(tmp_vectors);
+    krylov_space_signal(basis_vectors);
 
   // in case of failure: throw exception
   AssertThrow(iteration_state == SolverControl::success,
-              SolverControl::NoConvergence(accumulated_iterations, last_res));
+              SolverControl::NoConvergence(accumulated_iterations, res));
 }
 
 
@@ -1721,7 +1710,7 @@ SolverFGMRES<VectorType>::solve(const MatrixType         &A,
 
   // Generate an object where basis vectors are stored.
   typename internal::SolverGMRESImplementation::TmpVectors<VectorType> v(
-    basis_size, this->memory);
+    basis_size + 1, this->memory);
   typename internal::SolverGMRESImplementation::TmpVectors<VectorType> z(
     basis_size, this->memory);
 
@@ -1731,84 +1720,85 @@ SolverFGMRES<VectorType>::solve(const MatrixType         &A,
 
   // matrix used for the orthogonalization process later
   H.reinit(basis_size + 1, basis_size);
+  std::vector<std::pair<double, double>> givens_rotations(basis_size);
 
   Vector<double> h(basis_size + 1);
 
   // Vectors for projected system
-  Vector<double> projected_rhs;
-  Vector<double> y;
+  Vector<double> projected_rhs(basis_size + 1);
+  Vector<double> y(basis_size);
 
   // Iteration starts here
   double res = std::numeric_limits<double>::lowest();
 
-  typename VectorMemory<VectorType>::Pointer aux(this->memory);
-  aux->reinit(x);
   do
     {
-      A.vmult(*aux, x);
-      aux->sadd(-1., 1., b);
+      A.vmult(v(0, x), x);
+      v[0].sadd(-1., 1., b);
 
-      double beta     = aux->l2_norm();
-      res             = beta;
+      double norm_v   = v[0].l2_norm();
+      res             = norm_v;
       iteration_state = this->iteration_status(accumulated_iterations, res, x);
       if (iteration_state == SolverControl::success)
         break;
 
       H.reinit(basis_size + 1, basis_size);
-      double a = beta;
 
-      for (unsigned int j = 0; j < basis_size; ++j)
+      projected_rhs(0)             = norm_v;
+      unsigned int inner_iteration = 0;
+      for (; (inner_iteration < basis_size &&
+              iteration_state == SolverControl::iterate);
+           ++inner_iteration)
         {
-          if (a != 0) // treat lucky breakdown
-            v(j, x).equ(1. / a, *aux);
-          else
-            v(j, x) = 0.;
-
+          // norm_v = 0 is a lucky breakdown, the solver will reach
+          // convergence, but we must not divide by zero here.
+          if (norm_v != 0)
+            v[inner_iteration] /= norm_v;
 
-          preconditioner.vmult(z(j, x), v[j]);
-          A.vmult(*aux, z[j]);
+          preconditioner.vmult(z(inner_iteration, x), v[inner_iteration]);
+          A.vmult(v(inner_iteration + 1, x), z[inner_iteration]);
 
           // Gram-Schmidt
-          bool         re_orthogonalize = false;
-          const double s =
-            internal::SolverGMRESImplementation::iterated_gram_schmidt<
-              VectorType>(additional_data.orthogonalization_strategy,
-                          v,
-                          j + 1,
-                          0,
-                          *aux,
-                          h,
-                          re_orthogonalize);
-          for (unsigned int i = 0; i <= j; ++i)
-            H(i, j) = h(i);
-          H(j + 1, j) = a = s;
+          bool re_orthogonalize = false;
+          norm_v = internal::SolverGMRESImplementation::iterated_gram_schmidt<
+            VectorType>(additional_data.orthogonalization_strategy,
+                        v,
+                        inner_iteration + 1,
+                        0,
+                        v[inner_iteration + 1],
+                        h,
+                        re_orthogonalize);
 
           // Compute projected solution
+          h(inner_iteration + 1) = norm_v;
+          internal::SolverGMRESImplementation::givens_rotation(h,
+                                                               projected_rhs,
+                                                               givens_rotations,
+                                                               inner_iteration);
+
+          //  append vector on Hessenberg matrix
+          for (unsigned int i = 0; i < inner_iteration + 1; ++i)
+            H(i, inner_iteration) = h(i);
 
-          if (j > 0)
-            {
-              H1.reinit(j + 1, j);
-              projected_rhs.reinit(j + 1);
-              y.reinit(j);
-              projected_rhs(0) = beta;
-              H1.fill(H);
-
-              // check convergence. note that the vector 'x' we pass to the
-              // criterion is not the final solution we compute if we
-              // decide to jump out of the iteration (we update 'x' again
-              // right after the current loop)
-              Householder<double> house(H1);
-              res = house.least_squares(y, projected_rhs);
-              iteration_state =
-                this->iteration_status(++accumulated_iterations, res, x);
-              if (iteration_state != SolverControl::iterate)
-                break;
-            }
+          //  default residual
+          res = std::fabs(projected_rhs(inner_iteration + 1));
+
+          // check convergence. note that the vector 'x' we pass to the
+          // criterion is not the final solution we compute if we
+          // decide to jump out of the iteration (we update 'x' again
+          // right after the current loop)
+          iteration_state =
+            this->iteration_status(++accumulated_iterations, res, x);
         }
 
-      // Update solution vector
-      for (unsigned int j = 0; j < y.size(); ++j)
-        x.add(y(j), z[j]);
+      // Solve triangular system with projected quantities and update solution
+      // vector
+      internal::SolverGMRESImplementation::solve_triangular(inner_iteration,
+                                                            H,
+                                                            projected_rhs,
+                                                            y);
+      dealii::internal::SolverGMRESImplementation::add(
+        x, inner_iteration, y, z, false);
     }
   while (iteration_state == SolverControl::iterate);
 
index 5edf0cac0872e6f07afeadcf98f9cc77f0c1b508..f72504944df8d4765b5b0c6e8df0870db1cef008 100644 (file)
@@ -7,10 +7,10 @@ Refinement cycle 0
    Assembling Multigrid...
    Solving...
        Number of FGMRES iterations: 14
-       Total number of iterations used for approximation of A inverse: 45
-       Total number of iterations used for approximation of S inverse: 15
+       Total number of iterations used for approximation of A inverse: 42
+       Total number of iterations used for approximation of S inverse: 14
 
-   Note: The mean value was adjusted by 8.38229e-15
+   Note: The mean value was adjusted by 8.58753e-15
 
    Velocity L2 Error: 0.00340778
    Pressure L2 Error: 0.017624
@@ -23,10 +23,10 @@ Refinement cycle 1
    Assembling Multigrid...
    Solving...
        Number of FGMRES iterations: 16
-       Total number of iterations used for approximation of A inverse: 52
-       Total number of iterations used for approximation of S inverse: 17
+       Total number of iterations used for approximation of A inverse: 49
+       Total number of iterations used for approximation of S inverse: 16
 
-   Note: The mean value was adjusted by -1.75663e-14
+   Note: The mean value was adjusted by -1.77281e-14
 
    Velocity L2 Error: 0.000426317
    Pressure L2 Error: 0.00414007
@@ -39,10 +39,10 @@ Refinement cycle 2
    Assembling Multigrid...
    Solving...
        Number of FGMRES iterations: 16
-       Total number of iterations used for approximation of A inverse: 52
-       Total number of iterations used for approximation of S inverse: 17
+       Total number of iterations used for approximation of A inverse: 49
+       Total number of iterations used for approximation of S inverse: 16
 
-   Note: The mean value was adjusted by -2.40243e-14
+   Note: The mean value was adjusted by -2.39441e-14
 
    Velocity L2 Error: 5.33246e-05
    Pressure L2 Error: 0.00102063
@@ -56,10 +56,10 @@ Refinement cycle 0
    Solving...   Computing preconditioner...
 
        Number of FGMRES iterations: 14
-       Total number of iterations used for approximation of A inverse: 75
-       Total number of iterations used for approximation of S inverse: 15
+       Total number of iterations used for approximation of A inverse: 70
+       Total number of iterations used for approximation of S inverse: 14
 
-   Note: The mean value was adjusted by 3.10407e-15
+   Note: The mean value was adjusted by 3.04715e-15
 
    Velocity L2 Error: 0.00340778
    Pressure L2 Error: 0.017624
@@ -72,10 +72,10 @@ Refinement cycle 1
    Solving...   Computing preconditioner...
 
        Number of FGMRES iterations: 16
-       Total number of iterations used for approximation of A inverse: 134
-       Total number of iterations used for approximation of S inverse: 17
+       Total number of iterations used for approximation of A inverse: 126
+       Total number of iterations used for approximation of S inverse: 16
 
-   Note: The mean value was adjusted by -3.36177e-15
+   Note: The mean value was adjusted by -3.40459e-15
 
    Velocity L2 Error: 0.000426317
    Pressure L2 Error: 0.00414007
@@ -88,10 +88,10 @@ Refinement cycle 2
    Solving...   Computing preconditioner...
 
        Number of FGMRES iterations: 16
-       Total number of iterations used for approximation of A inverse: 229
-       Total number of iterations used for approximation of S inverse: 17
+       Total number of iterations used for approximation of A inverse: 216
+       Total number of iterations used for approximation of S inverse: 16
 
-   Note: The mean value was adjusted by -1.16485e-15
+   Note: The mean value was adjusted by -1.33296e-15
 
    Velocity L2 Error: 5.33247e-05
    Pressure L2 Error: 0.00102063
index 077cae729883de814cd3a58a06af608f1255dbfa..78a5b995a0579bddc71e280f6932437c407190f9 100644 (file)
@@ -69,7 +69,7 @@ test(unsigned int variant, unsigned int min_convergence_steps)
 
   SolverControl control(1000, 1e2 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
-  data.max_n_tmp_vectors = 80;
+  data.max_basis_size = 80;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {
index dd44f3447656ae777fdc2eeeb81ecfbcf4d5b211..572485a7e67f9bc259488f252e0391ad5d1b842c 100644 (file)
@@ -43,7 +43,7 @@ test()
 
   SolverControl control(1000, 1e3 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
-  data.max_n_tmp_vectors = 202;
+  data.max_basis_size = 200;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {
index 399cc5c7894d564f693d53742b2b442f632f912e..bccc1a35244d671bb9c02281498942e95924ced7 100644 (file)
@@ -67,7 +67,7 @@ test(unsigned int variant)
 
   SolverControl control(1000, 1e2 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
-  data.max_n_tmp_vectors          = 80;
+  data.max_basis_size             = 80;
   data.force_re_orthogonalization = true;
 
   SolverGMRES<Vector<number>> solver(control, data);
index 18cb893b5df4b4dcdc02c0445ee989c4d260b873..fd41939e9244ce7115882c5611ba01831610149a 100644 (file)
@@ -73,7 +73,7 @@ test(const unsigned int n_expected_steps)
   // of re-orthogonalization)
   SolverControl control(1000, 1e3 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
-  data.max_n_tmp_vectors = 202;
+  data.max_basis_size = 200;
 
   SolverGMRES<Vector<number>> solver(control, data);
   auto print_re_orthogonalization = [](int accumulated_iterations) {
index 1743b1910c36d5f24ee85c38d01ca7c591a40d4f..37774657d89112e14e20c6b308c8aa565df6d52a 100644 (file)
@@ -70,7 +70,7 @@ test()
 
   SolverControl control(1000, 1e2 * std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number>>::AdditionalData data;
-  data.max_n_tmp_vectors          = 202;
+  data.max_basis_size             = 200;
   data.force_re_orthogonalization = true;
 
   SolverGMRES<Vector<number>> solver(control, data);
index 5a4abcd39b132488f3490de7044811fda9e2be17..800b4bffdbbbf9da35a1d04af5e344b5da78e5af 100644 (file)
@@ -89,9 +89,9 @@ main()
   SolverControl                 control(100, 1.e-3);
   SolverControl                 verbose_control(100, 1.e-3, true);
   SolverCG<>                    cg(control, mem);
-  SolverGMRES<>::AdditionalData data1(8);
+  SolverGMRES<>::AdditionalData data1(6);
   SolverGMRES<>                 gmres(control, mem, data1);
-  SolverGMRES<>::AdditionalData data2(8, true);
+  SolverGMRES<>::AdditionalData data2(6, true);
   SolverGMRES<>                 gmresright(control, mem, data2);
   SolverMinRes<>                minres(control, mem);
   SolverBicgstab<>              bicgstab(control, mem);
@@ -99,7 +99,7 @@ main()
   SolverQMRS<>                  qmrs(control, mem);
   SolverFIRE<>                  fire(control, mem);
 
-  SolverGMRES<>::AdditionalData data3(8);
+  SolverGMRES<>::AdditionalData data3(6);
   data3.orthogonalization_strategy =
     LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
   SolverGMRES<> gmresclassical(control, mem, data3);
index 653812ce53b78a127b6ce3ccbab61d699f543fe9..4dd322bb37a372271d48995c20a4f02460304935 100644 (file)
@@ -110,13 +110,13 @@ DEAL:no-fail:Bicgstab::Failure step 10 value 0.001961
 DEAL:no-fail::Exception: SolverControl::NoConvergence(state.last_step, state.last_residual)
 DEAL:no-fail:GMRES::Starting value 11.00
 DEAL:no-fail:GMRES::Failure step 10 value 0.8414
-DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res)
+DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res)
 DEAL:no-fail:GMRES::Starting value 11.00
 DEAL:no-fail:GMRES::Failure step 10 value 0.8414
-DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res)
+DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res)
 DEAL:no-fail:GMRES::Starting value 11.00
 DEAL:no-fail:GMRES::Failure step 10 value 0.8414
-DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, last_res)
+DEAL:no-fail::Exception: SolverControl::NoConvergence(accumulated_iterations, res)
 DEAL:no-fail:SQMR::Starting value 11.00
 DEAL:no-fail:SQMR::Failure step 10 value 0.4215
 DEAL:no-fail::Exception: SolverControl::NoConvergence(step, state.last_residual)
index e140dc5c2f60830c4194a2d54cea0e5400e3d7b2..f788c3e6588829f829be748404fb4bbc8f26a5d7 100644 (file)
@@ -105,7 +105,7 @@ main()
     deallog << "Solve with PreconditionIdentity: " << std::endl;
     SolverControl control(40, 1e-4);
     SolverGMRES<LinearAlgebra::distributed::Vector<double>>::AdditionalData
-      data3(8);
+      data3(6);
     data3.orthogonalization_strategy =
       LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<LinearAlgebra::distributed::Vector<double>> solver(control,
@@ -127,7 +127,7 @@ main()
     deallog << "Solve with PreconditionIdentity: " << std::endl;
     SolverControl control(40, 1e-4);
     SolverGMRES<LinearAlgebra::distributed::BlockVector<double>>::AdditionalData
-      data3(8);
+      data3(6);
     data3.orthogonalization_strategy =
       LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<LinearAlgebra::distributed::BlockVector<double>> solver(control,
index bcf6aee22f71a550a6579dee9fd386575241eaa6..3dd6103bc915af8b7b14b6bb1943e9727bc20b37 100644 (file)
@@ -103,7 +103,7 @@ main()
   {
     deallog << "Solve with PreconditionIdentity: " << std::endl;
     SolverControl                               control(40, 1e-4);
-    SolverGMRES<Vector<double>>::AdditionalData data3(8);
+    SolverGMRES<Vector<double>>::AdditionalData data3(6);
     data3.orthogonalization_strategy =
       LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<Vector<double>> solver(control, data3);
@@ -123,7 +123,7 @@ main()
 
     deallog << "Solve with PreconditionIdentity: " << std::endl;
     SolverControl                                    control(40, 1e-4);
-    SolverGMRES<BlockVector<double>>::AdditionalData data3(8);
+    SolverGMRES<BlockVector<double>>::AdditionalData data3(6);
     data3.orthogonalization_strategy =
       LinearAlgebra::OrthogonalizationStrategy::classical_gram_schmidt;
     SolverGMRES<BlockVector<double>> solver(control, data3);
index 4a828b0aecd338ccedc4ef14d82adc83ad939cc1..22d0c8ed46e7a6e6fc2fd89b10172344e956934d 100644 (file)
@@ -97,7 +97,7 @@ main()
 
   SolverGMRES<> gmres(control,
                       mem,
-                      SolverGMRES<>::AdditionalData(/*max_vecs=*/8));
+                      SolverGMRES<>::AdditionalData(/*basis_size=*/6));
   gmres.connect(&monitor_norm);
   gmres.connect(&monitor_mean);
 
index 8f2bd5c7499b0248b61ede61c2ef6e47503d2e84..d936b867f3e8d377efa02360a73b0faa2739f46b 100644 (file)
@@ -13,6 +13,6 @@ DEAL:cg::Convergence step 39 value 7.335e-08
 DEAL:Bicgstab::Starting value 36.00
 DEAL:Bicgstab::Convergence step 28 value 9.087e-08
 DEAL:GMRES::Starting value 34.45
-DEAL:GMRES::Convergence step 39 value 6.627e-08
+DEAL:GMRES::Convergence step 40 value 9.928e-08
 DEAL:FGMRES::Starting value 36.00
-DEAL:FGMRES::Convergence step 39 value 7.890e-08
+DEAL:FGMRES::Convergence step 40 value 9.803e-08
index a716a5ebb1c7d4357a72d4f2c843b87af84d6a73..d4c88c708aac376ef3b41031a9f475a84935d12d 100644 (file)
@@ -3,4 +3,4 @@ DEAL::Size 37 Unknowns 1296
 DEAL:cg::Starting value 36.00
 DEAL:cg::Convergence step 39 value 7.335e-08
 DEAL:GMRES::Starting value 34.45
-DEAL:GMRES::Convergence step 39 value 6.627e-08
+DEAL:GMRES::Convergence step 40 value 9.928e-08
index 48089657bf2ad5be99b9fc8374db0fb7ba54701b..761e9a5c1553b82cd214dd23c3ff1e1aeb330259 100644 (file)
@@ -154,7 +154,8 @@ main()
       solver_cg.solve(A, u, f, PreconditionIdentity());
 
       u = 0;
-      SolverGMRES<> solver_gmres(solver_control);
+      SolverGMRES<> solver_gmres(solver_control,
+                                 SolverGMRES<>::AdditionalData(28));
       // Attach all possible slots.
       solver_gmres.connect_condition_number_slot(
         std::bind(output_double_number,
index f597ef56d0922d0156cf9e64bfd0c16e147167f3..ac1b0228b1f4a2362e7ca20a6485e4a6c4d1bfc3 100644 (file)
@@ -128,36 +128,35 @@ DEAL:GMRES::Convergence step 65 value 0.0009799
 DEAL:GMRES::Eigenvalues:  (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000)
 DEAL:GMRES::Final condition number: 288.2
 DEAL:GMRES::Hessenberg matrix: 
-DEAL:GMRES::0.1982 0.8628 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0.8628 6.098 1.971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 1.971 4.008 1.721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 1.721 3.583 2.004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 2.004 4.008 1.910 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 1.910 3.325 1.950 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 1.950 4.036 1.744 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 1.744 3.951 1.992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 1.992 4.080 1.994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 1.920 4.007 2.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 2.002 3.981 2.022 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 2.022 4.028 1.908 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 1.908 4.238 1.925 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 1.925 4.011 2.042 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 2.042 3.827 2.111 0 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.111 3.742 2.046 0 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 4.070 1.879 0 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.879 4.168 1.933 0 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.933 4.013 1.886 0 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.886 3.963 1.957 0 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.957 3.818 2.033 0 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.033 4.006 2.201 0 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.201 3.944 1.862 0 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.862 4.039 1.827 0 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.827 3.704 1.844 0 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.844 4.073 2.046 0 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 3.592 1.647 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.647 3.431 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.155 0 
-DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0.1982 0.8628 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0.8628 6.098 1.971 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 1.971 4.008 1.721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 1.721 3.583 2.004 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 2.004 4.008 1.910 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 1.910 3.325 1.950 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 1.950 4.036 1.744 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 1.744 3.951 1.992 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 1.992 4.080 1.994 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 1.920 4.007 2.002 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 2.002 3.981 2.022 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 2.022 4.028 1.908 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 1.908 4.238 1.925 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 1.925 4.011 2.042 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 2.042 3.827 2.111 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.111 3.742 2.046 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 4.070 1.879 0 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.879 4.168 1.933 0 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.933 4.013 1.886 0 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.886 3.963 1.957 0 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.957 3.818 2.033 0 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.033 4.006 2.201 0 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.201 3.944 1.862 0 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.862 4.039 1.827 0 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.827 3.704 1.844 0 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.844 4.073 2.046 0 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.046 3.592 1.647 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.647 3.431 
+DEAL:GMRES::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.155 
 DEAL:GMRES::Final Eigenvalues:  (0.02681,0.000) (0.2333,0.000) (1.178,0.000) (2.161,0.000) (3.554,0.000) (4.930,0.000) (6.307,0.000) (7.170,0.000) (7.726,0.000)
 DEAL:GMRES::Final condition number: 288.2
 DEAL:GMRES::Arnoldi vectors norms: 
index 51b6d4c9080463144d1211c639bea99226fa02a1..c685e92c4dc9891bd13ec9d7cf9ea3a125f81f75 100644 (file)
@@ -1386,13 +1386,19 @@ namespace StokesClass
     solution.block(block_vel) = distributed_stokes_solution.block(block_vel);
     solution.block(block_p)   = distributed_stokes_solution.block(block_p);
 
+    LinearAlgebra::distributed::BlockVector<double> vec(
+      distributed_stokes_solution);
+    stokes_matrix.vmult(vec, distributed_stokes_solution);
+    vec -= distributed_stokes_rhs;
+
     deallog << "Solved-in "
             << (solver_control_cheap.last_step() !=
                     numbers::invalid_unsigned_int ?
                   solver_control_cheap.last_step() :
                   0)
             << " iterations, final residual: "
-            << solver_control_cheap.last_value() << std::endl;
+            << solver_control_cheap.last_value() << " and true residual "
+            << vec.l2_norm() << std::endl;
   }
 
 
index 1f2687b42cdf4a1cc4a9c6afdc3fc59062a0d17e..ba7275b72b06c9e5bcbe0025267ac16eea8122b0 100644 (file)
@@ -4,16 +4,16 @@ DEAL:2d::n_sinker: 4       max/min viscosity ratio: 1000.00
 DEAL:2d::
 DEAL:2d::Number of active cells: 1024 (on 6 levels)
 DEAL:2d::Number of degrees of freedom: 9539 (8450+1089)
-DEAL:2d::Solved-in 69 iterations, final residual: 5.58979e-08
+DEAL:2d::Solved-in 64 iterations, final residual: 5.31067e-08 and true residual 5.31067e-08
 DEAL:2d::
 DEAL:2d::Number of active cells: 4096 (on 7 levels)
 DEAL:2d::Number of degrees of freedom: 37507 (33282+4225)
-DEAL:2d::Solved-in 69 iterations, final residual: 2.58019e-08
+DEAL:2d::Solved-in 64 iterations, final residual: 2.55967e-08 and true residual 2.55967e-08
 DEAL:2d::
 DEAL:3d::Sinker problem in 3D.
 DEAL:3d::n_sinker: 4       max/min viscosity ratio: 1000.00
 DEAL:3d::
 DEAL:3d::Number of active cells: 512 (on 4 levels)
 DEAL:3d::Number of degrees of freedom: 15468 (14739+729)
-DEAL:3d::Solved-in 62 iterations, final residual: 1.99453e-08
+DEAL:3d::Solved-in 58 iterations, final residual: 2.08315e-08 and true residual 2.08315e-08
 DEAL:3d::
index 1f2687b42cdf4a1cc4a9c6afdc3fc59062a0d17e..46baf3fb41a80b4823a30857bb4971165984931c 100644 (file)
@@ -4,16 +4,16 @@ DEAL:2d::n_sinker: 4       max/min viscosity ratio: 1000.00
 DEAL:2d::
 DEAL:2d::Number of active cells: 1024 (on 6 levels)
 DEAL:2d::Number of degrees of freedom: 9539 (8450+1089)
-DEAL:2d::Solved-in 69 iterations, final residual: 5.58979e-08
+DEAL:2d::Solved-in 64 iterations, final residual: 5.31067e-08
 DEAL:2d::
 DEAL:2d::Number of active cells: 4096 (on 7 levels)
 DEAL:2d::Number of degrees of freedom: 37507 (33282+4225)
-DEAL:2d::Solved-in 69 iterations, final residual: 2.58019e-08
+DEAL:2d::Solved-in 64 iterations, final residual: 2.55967e-08
 DEAL:2d::
 DEAL:3d::Sinker problem in 3D.
 DEAL:3d::n_sinker: 4       max/min viscosity ratio: 1000.00
 DEAL:3d::
 DEAL:3d::Number of active cells: 512 (on 4 levels)
 DEAL:3d::Number of degrees of freedom: 15468 (14739+729)
-DEAL:3d::Solved-in 62 iterations, final residual: 1.99453e-08
+DEAL:3d::Solved-in 58 iterations, final residual: 2.08315e-08
 DEAL:3d::
index 002798281bd7041aa97b1a9aac1cc70982112b5e..e271e6ef0834a78ebae270ee050cd891f26c224a 100644 (file)
@@ -8,7 +8,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 9.41508
-DEAL:GMRES::Convergence step 10 value 7.45878e-13
+DEAL:GMRES::Convergence step 10 value 7.45877e-13
 DEAL::energy-error: 0.439211
 DEAL::L2-error:     0.0109342
 DEAL::Estimate 0.979555
@@ -22,7 +22,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 10.0041
-DEAL:GMRES::Convergence step 21 value 1.79039e-13
+DEAL:GMRES::Convergence step 21 value 1.79024e-13
 DEAL::energy-error: 0.332582
 DEAL::L2-error:     0.00548996
 DEAL::Estimate 0.838060
@@ -36,7 +36,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 10.8018
-DEAL:GMRES::Convergence step 28 value 2.71098e-13
+DEAL:GMRES::Convergence step 28 value 2.71118e-13
 DEAL::energy-error: 0.237705
 DEAL::L2-error:     0.00250869
 DEAL::Estimate 0.611449
@@ -50,7 +50,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 13.4185
-DEAL:GMRES::Convergence step 32 value 6.11150e-13
+DEAL:GMRES::Convergence step 31 value 9.79574e-13
 DEAL::energy-error: 0.170175
 DEAL::L2-error:     0.00120328
 DEAL::Estimate 0.448676
@@ -64,7 +64,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 16.6112
-DEAL:GMRES::Convergence step 32 value 7.81623e-13
+DEAL:GMRES::Convergence step 32 value 7.99534e-13
 DEAL::energy-error: 0.121878
 DEAL::L2-error:     0.000591777
 DEAL::Estimate 0.327296
@@ -78,7 +78,7 @@ DEAL::Assemble multilevel matrix
 DEAL::Assemble right hand side
 DEAL::Solve
 DEAL:GMRES::Starting value 19.3095
-DEAL:GMRES::Convergence step 31 value 8.40512e-13
+DEAL:GMRES::Convergence step 31 value 9.61737e-13
 DEAL::energy-error: 0.0863235
 DEAL::L2-error:     0.000292807
 DEAL::Estimate 0.232279
index 72980d7cfd3d9aea17b48de319be6b0ea34e8a52..23504a72e7187e30f3a1281af156519ea5161199 100644 (file)
@@ -6,29 +6,29 @@ DEAL:0::   Number of degrees of freedom: 9539 (8450+1089 -- 0+0)
 DEAL:0::Tracer particles: 1504
 DEAL:0::Solid particles: 9216
 DEAL:0:FGMRES::Starting value 30.9301
-DEAL:0:FGMRES::Convergence step 147 value 2.63430e-09
-DEAL:0::   Solved in 147 iterations.
+DEAL:0:FGMRES::Convergence step 148 value 2.35433e-09
+DEAL:0::   Solved in 148 iterations.
 DEAL:0::   Number of degrees of freedom: 9845 (8722+1123 -- 9216+1504)
 DEAL:0::Cycle 1:
 DEAL:0::Time : 0.00250000, time step: 0.00250000
 DEAL:0:FGMRES::Starting value 2.92442
-DEAL:0:FGMRES::Convergence step 141 value 5.01043e-09
-DEAL:0::   Solved in 141 iterations.
+DEAL:0:FGMRES::Convergence step 146 value 5.16889e-09
+DEAL:0::   Solved in 146 iterations.
 DEAL:0::Cycle 2:
 DEAL:0::Time : 0.00500000, time step: 0.00250000
 DEAL:0:FGMRES::Starting value 0.637099
-DEAL:0:FGMRES::Convergence step 121 value 4.79602e-09
-DEAL:0::   Solved in 121 iterations.
+DEAL:0:FGMRES::Convergence step 117 value 5.22134e-09
+DEAL:0::   Solved in 117 iterations.
 DEAL:0::Cycle 3:
 DEAL:0::Time : 0.00750000, time step: 0.00250000
 DEAL:0:FGMRES::Starting value 0.649739
-DEAL:0:FGMRES::Convergence step 126 value 5.21191e-09
-DEAL:0::   Solved in 126 iterations.
+DEAL:0:FGMRES::Convergence step 125 value 5.07375e-09
+DEAL:0::   Solved in 125 iterations.
 DEAL:0::Cycle 4:
 DEAL:0::Time : 0.0100000, time step: 0.00250000
 DEAL:0:FGMRES::Starting value 0.744112
-DEAL:0:FGMRES::Convergence step 131 value 4.88195e-09
-DEAL:0::   Solved in 131 iterations.
+DEAL:0:FGMRES::Convergence step 130 value 5.49751e-09
+DEAL:0::   Solved in 130 iterations.
 
 DEAL:1::Running StokesImmersedProblem<2> using Trilinos.
 DEAL:1::Cycle 0:
@@ -37,27 +37,27 @@ DEAL:1::   Number of degrees of freedom: 9539 (8450+1089 -- 0+0)
 DEAL:1::Tracer particles: 1504
 DEAL:1::Solid particles: 9216
 DEAL:1:FGMRES::Starting value 30.9301
-DEAL:1:FGMRES::Convergence step 147 value 2.63430e-09
-DEAL:1::   Solved in 147 iterations.
+DEAL:1:FGMRES::Convergence step 148 value 2.35433e-09
+DEAL:1::   Solved in 148 iterations.
 DEAL:1::   Number of degrees of freedom: 9845 (8722+1123 -- 9216+1504)
 DEAL:1::Cycle 1:
 DEAL:1::Time : 0.00250000, time step: 0.00250000
 DEAL:1:FGMRES::Starting value 2.92442
-DEAL:1:FGMRES::Convergence step 141 value 5.01043e-09
-DEAL:1::   Solved in 141 iterations.
+DEAL:1:FGMRES::Convergence step 146 value 5.16889e-09
+DEAL:1::   Solved in 146 iterations.
 DEAL:1::Cycle 2:
 DEAL:1::Time : 0.00500000, time step: 0.00250000
 DEAL:1:FGMRES::Starting value 0.637099
-DEAL:1:FGMRES::Convergence step 121 value 4.79602e-09
-DEAL:1::   Solved in 121 iterations.
+DEAL:1:FGMRES::Convergence step 117 value 5.22134e-09
+DEAL:1::   Solved in 117 iterations.
 DEAL:1::Cycle 3:
 DEAL:1::Time : 0.00750000, time step: 0.00250000
 DEAL:1:FGMRES::Starting value 0.649739
-DEAL:1:FGMRES::Convergence step 126 value 5.21191e-09
-DEAL:1::   Solved in 126 iterations.
+DEAL:1:FGMRES::Convergence step 125 value 5.07375e-09
+DEAL:1::   Solved in 125 iterations.
 DEAL:1::Cycle 4:
 DEAL:1::Time : 0.0100000, time step: 0.00250000
 DEAL:1:FGMRES::Starting value 0.744112
-DEAL:1:FGMRES::Convergence step 131 value 4.88195e-09
-DEAL:1::   Solved in 131 iterations.
+DEAL:1:FGMRES::Convergence step 130 value 5.49751e-09
+DEAL:1::   Solved in 130 iterations.
 
index 8570662278b3e621b82ed8a382bf17413376f709..ffa1065eb8793450c581922cca5511c93e63bf30 100644 (file)
@@ -82,7 +82,7 @@ main(int argc, char **argv)
 
     check_solver_within_range(solver.solve(A, u, f, preconditioner),
                               control.last_step(),
-                              74,
+                              71,
                               76);
   }
 }
index 5b6f34323b5067354d98a2ccde4355d5f9ce520d..3e4091158b3e03f47589f7267e0909f4c3f406fb 100644 (file)
@@ -1,4 +1,4 @@
 
 DEAL::Size 32 Unknowns 961
 DEAL::Solver type: N6dealii11SolverGMRESINS_16TrilinosWrappers3MPI6VectorEEE
-DEAL::Solver stopped within 74 - 76 iterations
+DEAL::Solver stopped within 71 - 76 iterations
index 56085099d12c96ae57e39b90a9e64acd2190b68d..1b926dae5c42d7561797ec752f181c3323ab340e 100644 (file)
@@ -75,8 +75,11 @@ main(int argc, char **argv)
     u.compress(VectorOperation::insert);
 
     GrowingVectorMemory<LinearAlgebra::TpetraWrappers::Vector<double>> mem;
+    SolverGMRES<LinearAlgebra::TpetraWrappers::Vector<double>>::AdditionalData
+                                                               sol_data(28);
     SolverGMRES<LinearAlgebra::TpetraWrappers::Vector<double>> solver(control,
-                                                                      mem);
+                                                                      mem,
+                                                                      sol_data);
     PreconditionIdentity                                       preconditioner;
 
     deallog << "Solver type: " << typeid(solver).name() << std::endl;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.