From e89c58bf2b0f3714104a76a248bd394f19d05e3a Mon Sep 17 00:00:00 2001 From: kronbichler Date: Sun, 11 May 2014 09:44:02 +0000 Subject: [PATCH] Compute eigenvalues with GMRES git-svn-id: https://svn.dealii.org/trunk@32897 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 ++ deal.II/include/deal.II/lac/solver_gmres.h | 48 +++++++++- tests/lac/gmres_eigenvalues.cc | 100 +++++++++++++++++++++ tests/lac/gmres_eigenvalues.output | 16 ++++ 4 files changed, 168 insertions(+), 3 deletions(-) create mode 100644 tests/lac/gmres_eigenvalues.cc create mode 100644 tests/lac/gmres_eigenvalues.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 92d4fcfcd6..8e6275ff0b 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,13 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/05/11) +
    2. +
  2. 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. diff --git a/deal.II/include/deal.II/lac/solver_gmres.h b/deal.II/include/deal.II/lac/solver_gmres.h index 448244f6eb..84133ceace 100644 --- a/deal.II/include/deal.II/lac/solver_gmres.h +++ b/deal.II/include/deal.II/lac/solver_gmres.h @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -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::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::solve (const MATRIX &A, // restart unsigned int accumulated_iterations = 0; + // for eigenvalue computation, need to collect the Hessenberg matrix (before + // applying Givens rotations) + FullMatrix 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::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::solve (const MATRIX &A, for (unsigned int j=0; j mat(dim,dim); + for (unsigned int i=0; i +#include +#include +#include + + + +template +void test (unsigned int variant) +{ + const unsigned int n = 64; + Vector rhs(n), sol(n); + rhs = 1.; + + FullMatrix 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::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 >::AdditionalData data; + data.max_n_tmp_vectors = 80; + data.compute_eigenvalues = true; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + if (variant == 0) + { + typename SolverCG >::AdditionalData cg_data; + cg_data.compute_eigenvalues = true; + SolverCG > 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(0); + test(1); + test(2); + test(3); + deallog.pop(); +} + diff --git a/tests/lac/gmres_eigenvalues.output b/tests/lac/gmres_eigenvalues.output new file mode 100644 index 0000000000..0c5aa5a46f --- /dev/null +++ b/tests/lac/gmres_eigenvalues.output @@ -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) -- 2.39.5