<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.
#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>
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
* 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;
};
/**
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)
{}
// 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);
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);
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)
--- /dev/null
+// ---------------------------------------------------------------------
+// $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();
+}
+
--- /dev/null
+
+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)
projects / dealii.git / commitdiff