]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Compute eigenvalues with GMRES
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 11 May 2014 09:44:02 +0000 (09:44 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 11 May 2014 09:44:02 +0000 (09:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@32897 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/solver_gmres.h
tests/lac/gmres_eigenvalues.cc [new file with mode: 0644]
tests/lac/gmres_eigenvalues.output [new file with mode: 0644]

index 92d4fcfcd68719ac3c7214323150c9210ec45979..8e6275ff0b9f450c755948cd302e17d4ce40fe8f 100644 (file)
@@ -148,6 +148,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: The GMRES solver of deal.II can now write an estimate of
+  eigenvalues to the log file, in analogy to the CG solver. This is enabled
+  by the flag SolverGMRES<>::AdditionalData::compute_eigenvalues.
+  <br>
+  (Martin Kronbichler, 2014/05/11)
+  </li>  
+
   <li> New: The GridIn::read_vtk() function places fewer restrictions
   on the VTK files it wants to read and should, consequently, be able
   to read more correctly formatted VTK files than before.
index 448244f6eb93813f2383a33f058ac3b3554cfd18..84133ceacefafe1d5944c720d25d1d7307d98c0e 100644 (file)
@@ -26,6 +26,7 @@
 #include <deal.II/lac/solver.h>
 #include <deal.II/lac/solver_control.h>
 #include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/lapack_full_matrix.h>
 #include <deal.II/lac/vector.h>
 
 #include <vector>
@@ -170,7 +171,8 @@ public:
     AdditionalData (const unsigned int max_n_tmp_vectors = 30,
                     const bool right_preconditioning = false,
                     const bool use_default_residual = true,
-                    const bool force_re_orthogonalization = false);
+                    const bool force_re_orthogonalization = false,
+                    const bool compute_eigenvalues = false);
 
     /**
      * Maximum number of temporary vectors. This parameter controls the size
@@ -200,6 +202,17 @@ public:
      * if necessary.
      */
     bool force_re_orthogonalization;
+
+    /**
+     * Compute all eigenvalues of the Hessenberg matrix generated while
+     * solving, i.e., the projected system matrix. This gives an approximation
+     * of the eigenvalues of the (preconditioned) system matrix. Since the
+     * Hessenberg matrix is thrown away at restart, the eigenvalues are
+     * printed for every 30 iterations.
+     *
+     * @note Requires LAPACK support.
+     */
+    bool compute_eigenvalues;
   };
 
   /**
@@ -439,12 +452,14 @@ SolverGMRES<VECTOR>::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         force_re_orthogonalization,
+                const bool         compute_eigenvalues)
   :
   max_n_tmp_vectors(max_n_tmp_vectors),
   right_preconditioning(right_preconditioning),
   use_default_residual(use_default_residual),
-  force_re_orthogonalization(force_re_orthogonalization)
+  force_re_orthogonalization(force_re_orthogonalization),
+  compute_eigenvalues (compute_eigenvalues)
 {}
 
 
@@ -581,6 +596,12 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
   // restart
   unsigned int accumulated_iterations = 0;
 
+  // for eigenvalue computation, need to collect the Hessenberg matrix (before
+  // applying Givens rotations)
+  FullMatrix<double> H_orig;
+  if (additional_data.compute_eigenvalues)
+    H_orig.reinit(n_tmp_vectors, n_tmp_vectors-1);
+
   // matrix used for the orthogonalization process later
   H.reinit(n_tmp_vectors, n_tmp_vectors-1);
 
@@ -727,6 +748,12 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
           if (numbers::is_finite(1./s))
             vv *= 1./s;
 
+          // for eigenvalues, get the resulting coefficients from the
+          // orthogonalization process
+          if (additional_data.compute_eigenvalues)
+            for (unsigned int i=0; i<dim+1; ++i)
+              H_orig(i,inner_iteration) = h(i);
+
           //  Transformation into triagonal structure
           givens_rotation(h,gamma,ci,si,inner_iteration);
 
@@ -795,6 +822,21 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
         for (unsigned int j=0; j<dim; ++j)
           H1(i,j) = H(i,j);
 
+      // compute eigenvalues from the Hessenberg matrix generated during the
+      // inner iterations
+      if (additional_data.compute_eigenvalues)
+        {
+          LAPACKFullMatrix<double> mat(dim,dim);
+          for (unsigned int i=0; i<dim; ++i)
+            for (unsigned int j=0; j<dim; ++j)
+              mat(i,j) = H_orig(i,j);
+          mat.compute_eigenvalues();
+          deallog << "Eigenvalue estimate: ";
+          for (unsigned int i=0; i<mat.n(); ++i)
+            deallog << ' ' << mat.eigenvalue(i);
+          deallog << std::endl;
+        }
+
       H1.backward(h,gamma);
 
       if (left_precondition)
diff --git a/tests/lac/gmres_eigenvalues.cc b/tests/lac/gmres_eigenvalues.cc
new file mode 100644 (file)
index 0000000..688190f
--- /dev/null
@@ -0,0 +1,100 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// test eigenvalue approximation by GMRES algorithm
+
+#include "../tests.h"
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/solver_gmres.h>
+#include <deal.II/lac/precondition.h>
+
+
+
+template <typename number>
+void test (unsigned int variant)
+{
+  const unsigned int n = 64;
+  Vector<number> rhs(n), sol(n);
+  rhs = 1.;
+
+  FullMatrix<number> matrix(n, n);
+
+  // put diagonal entries of different strengths. these are very challenging
+  // for GMRES and will usually take a lot of iterations until the Krylov
+  // subspace is complete enough
+  if (variant == 0)
+    for (unsigned int i=0; i<n; ++i)
+      matrix(i,i) = (i+1);
+  else if (variant == 1)
+    for (unsigned int i=0; i<n; ++i)
+      matrix(i,i) = (i+1) * (i+1) * (i+1) * (i+1);
+  else if (variant == 2)
+    for (unsigned int i=0; i<n; ++i)
+      matrix(i,i) = (i%2?1.:-1.)*(i+1);
+  else if (variant == 3)
+    for (unsigned int i=0; i<n; ++i)
+      {
+        matrix(i,i) = (i+1);
+        if (i<n-1)
+          matrix(i,i+1) = 1.5+i;
+        if (i<n-2)
+          matrix(i,i+2) = -0.25;
+      }
+  else
+    Assert(false, ExcMessage("Invalid variant"));
+  if (types_are_equal<number,float>::value == true)
+    Assert(variant < 4, ExcMessage("Invalid_variant"));
+
+  deallog.push(Utilities::int_to_string(variant,1));
+
+  SolverControl control(1000, variant==1?1e-4:1e-13);
+  typename SolverGMRES<Vector<number> >::AdditionalData data;
+  data.max_n_tmp_vectors = 80;
+  data.compute_eigenvalues = true;
+
+  SolverGMRES<Vector<number> > solver(control, data);
+  solver.solve(matrix, sol, rhs, PreconditionIdentity());
+
+  if (variant == 0)
+    {
+      typename SolverCG<Vector<number> >::AdditionalData cg_data;
+      cg_data.compute_eigenvalues = true;
+      SolverCG<Vector<number> > solver_cg(control, cg_data);
+      sol = 0;
+      solver_cg.solve(matrix, sol, rhs, PreconditionIdentity());
+    }
+
+  deallog.pop();
+}
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  deallog.push("double");
+  test<double>(0);
+  test<double>(1);
+  test<double>(2);
+  test<double>(3);
+  deallog.pop();
+}
+
diff --git a/tests/lac/gmres_eigenvalues.output b/tests/lac/gmres_eigenvalues.output
new file mode 100644 (file)
index 0000000..0c5aa5a
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:double:0:GMRES::Starting value 8.00
+DEAL:double:0:GMRES::Convergence step 56 value 0
+DEAL:double:0:GMRES::Eigenvalue estimate:  (64.0,0.00) (63.0,0.00) (62.0,0.00) (61.0,0.00) (60.0,0.00) (59.0,0.00) (58.0,0.00) (57.0,0.00) (56.0,0.00) (55.0,0.00) (1.00,0.00) (2.00,0.00) (54.0,0.00) (3.00,0.00) (53.0,0.00) (52.0,0.00) (51.0,0.00) (50.0,0.00) (49.0,0.00) (47.9,0.00) (46.7,0.00) (45.5,0.00) (4.00,0.00) (5.00,0.00) (44.3,0.00) (43.0,0.00) (41.7,0.00) (40.4,0.00) (6.00,0.00) (7.00,0.00) (39.0,0.00) (37.6,0.00) (8.00,0.00) (9.00,0.00) (36.1,0.00) (34.7,0.00) (10.0,0.00) (33.2,0.00) (11.0,0.00) (12.0,0.00) (31.8,0.00) (30.3,0.00) (28.9,0.00) (13.0,0.00) (14.0,0.00) (15.0,0.00) (16.0,0.00) (17.1,0.00) (18.2,0.00) (27.5,0.00) (19.5,0.00) (20.7,0.00) (22.0,0.00) (26.0,0.00) (23.3,0.00) (24.7,0.00)
+DEAL:double:0:cg::Starting value 8.00
+DEAL:double:0:cg::Convergence step 56 value 0
+DEAL:double:0:cg:: 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00 9.00 10.0 11.0 12.0 13.0 14.0 15.0 16.1 17.3 18.4 19.7 21.0 22.3 23.7 25.1 26.6 28.0 29.5 31.0 32.5 34.0 35.5 37.0 38.4 39.9 41.3 42.7 44.0 45.3 46.6 47.7 48.9 50.0 51.0 52.0 53.0 54.0 55.0 56.0 57.0 58.0 59.0 60.0 61.0 62.0 63.0 64.0
+DEAL:double:1:GMRES::Starting value 8.00
+DEAL:double:1:GMRES::Convergence step 64 value 4.41e-10
+DEAL:double:1:GMRES::Eigenvalue estimate:  (1.68e+07,0.00) (1.58e+07,0.00) (1.48e+07,0.00) (1.38e+07,0.00) (1.30e+07,0.00) (1.21e+07,0.00) (1.13e+07,0.00) (1.06e+07,0.00) (9.83e+06,0.00) (9.15e+06,0.00) (8.50e+06,0.00) (7.89e+06,0.00) (7.31e+06,0.00) (6.77e+06,0.00) (6.25e+06,0.00) (5.76e+06,0.00) (5.31e+06,0.00) (4.88e+06,0.00) (4.48e+06,0.00) (4.10e+06,0.00) (3.75e+06,0.00) (3.42e+06,0.00) (3.11e+06,0.00) (2.83e+06,0.00) (2.56e+06,0.00) (2.31e+06,0.00) (2.09e+06,0.00) (1.87e+06,0.00) (1.68e+06,0.00) (1.50e+06,0.00) (1.34e+06,0.00) (1.19e+06,0.00) (1.05e+06,0.00) (9.24e+05,0.00) (8.10e+05,0.00) (7.07e+05,0.00) (6.15e+05,0.00) (5.31e+05,0.00) (4.57e+05,0.00) (3.91e+05,0.00) (3.32e+05,0.00) (2.80e+05,0.00) (2.34e+05,0.00) (1.94e+05,0.00) (1.60e+05,0.00) (1.30e+05,0.00) (1.05e+05,0.00) (8.35e+04,0.00) (6.55e+04,0.00) (5.06e+04,0.00) (3.84e+04,0.00) (2.86e+04,0.00) (2.07e+04,0.00) (1.46e+04,0.00) (1.00e+04,0.00) (6.56e+03,0.00) (4.10e+03,0.00) (2.40e+03,0.00) (1.30e+03,0.00) (625.,0.00) (256.,0.00) (81.0,0.00) (1.00,0.00) (16.0,0.00)
+DEAL:double:2:GMRES::Starting value 8.00
+DEAL:double:2:GMRES::Convergence step 64 value 0
+DEAL:double:2:GMRES::Eigenvalue estimate:  (64.0,0.00) (62.0,0.00) (60.0,0.00) (58.0,0.00) (-63.0,0.00) (56.0,0.00) (54.0,0.00) (-61.0,0.00) (52.0,0.00) (-59.0,0.00) (50.0,0.00) (48.0,0.00) (-57.0,0.00) (46.0,0.00) (-55.0,0.00) (44.0,0.00) (42.0,0.00) (-53.0,0.00) (40.0,0.00) (-51.0,0.00) (-49.0,0.00) (38.0,0.00) (-47.0,0.00) (36.0,0.00) (-45.0,0.00) (-43.0,0.00) (34.0,0.00) (32.0,0.00) (30.0,0.00) (-41.0,0.00) (28.0,0.00) (-39.0,0.00) (26.0,0.00) (-37.0,0.00) (-35.0,0.00) (24.0,0.00) (22.0,0.00) (20.0,0.00) (18.0,0.00) (16.0,0.00) (-33.0,0.00) (14.0,0.00) (-31.0,0.00) (-29.0,0.00) (-27.0,0.00) (-25.0,0.00) (12.0,0.00) (-23.0,0.00) (-21.0,0.00) (10.0,0.00) (8.00,0.00) (6.00,0.00) (4.00,0.00) (2.00,0.00) (-19.0,0.00) (-17.0,0.00) (-15.0,0.00) (-13.0,0.00) (-1.00,0.00) (-3.00,0.00) (-11.0,0.00) (-5.00,0.00) (-9.00,0.00) (-7.00,0.00)
+DEAL:double:3:GMRES::Starting value 8.00
+DEAL:double:3:GMRES::Convergence step 64 value 0
+DEAL:double:3:GMRES::Eigenvalue estimate:  (67.4,1.30) (67.4,-1.30) (66.6,3.85) (66.6,-3.85) (65.0,6.22) (65.0,-6.22) (62.8,8.25) (62.8,-8.25) (60.3,9.81) (60.3,-9.81) (57.9,11.1) (57.9,-11.1) (55.2,12.4) (55.2,-12.4) (52.2,13.5) (52.2,-13.5) (48.8,14.1) (48.8,-14.1) (45.4,14.3) (45.4,-14.3) (41.9,13.9) (41.9,-13.9) (40.4,12.1) (40.4,-12.1) (37.7,13.7) (37.7,-13.7) (38.7,2.89) (38.7,-2.89) (34.3,13.4) (34.3,-13.4) (31.2,12.7) (31.2,-12.7) (28.4,11.8) (28.4,-11.8) (25.8,10.8) (25.8,-10.8) (23.4,9.65) (23.4,-9.65) (25.2,0.00) (21.1,8.46) (21.1,-8.46) (19.0,7.19) (19.0,-7.19) (17.1,5.86) (17.1,-5.86) (15.3,4.42) (15.3,-4.42) (14.3,2.78) (14.3,-2.78) (13.2,2.58) (13.2,-2.58) (11.5,1.65) (11.5,-1.65) (10.1,0.640) (10.1,-0.640) (8.92,0.00) (8.01,0.00) (7.00,0.00) (6.00,0.00) (5.00,0.00) (4.00,0.00) (3.00,0.00) (2.00,0.00) (1.00,0.00)

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.