From 9947df29cf366889faf946279e5ec59d3f7e24e2 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 6 May 2013 19:51:19 +0000 Subject: [PATCH] Choose to adaptively re-orthogonalize in GMRES, which cuts the number of inner products into one half for most applications. Add proper tests. git-svn-id: https://svn.dealii.org/trunk@29464 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/solver_gmres.h | 356 +++++++++--------- tests/lac/gmres_reorthogonalize_01.cc | 98 +++++ .../lac/gmres_reorthogonalize_01/cmp/generic | 26 ++ tests/lac/gmres_reorthogonalize_02.cc | 65 ++++ .../lac/gmres_reorthogonalize_02/cmp/generic | 5 + tests/lac/gmres_reorthogonalize_03.cc | 98 +++++ .../lac/gmres_reorthogonalize_03/cmp/generic | 21 ++ tests/lac/gmres_reorthogonalize_04.cc | 87 +++++ .../lac/gmres_reorthogonalize_04/cmp/generic | 5 + tests/lac/gmres_reorthogonalize_05.cc | 88 +++++ .../lac/gmres_reorthogonalize_05/cmp/generic | 5 + 11 files changed, 675 insertions(+), 179 deletions(-) create mode 100644 tests/lac/gmres_reorthogonalize_01.cc create mode 100644 tests/lac/gmres_reorthogonalize_01/cmp/generic create mode 100644 tests/lac/gmres_reorthogonalize_02.cc create mode 100644 tests/lac/gmres_reorthogonalize_02/cmp/generic create mode 100644 tests/lac/gmres_reorthogonalize_03.cc create mode 100644 tests/lac/gmres_reorthogonalize_03/cmp/generic create mode 100644 tests/lac/gmres_reorthogonalize_04.cc create mode 100644 tests/lac/gmres_reorthogonalize_04/cmp/generic create mode 100644 tests/lac/gmres_reorthogonalize_05.cc create mode 100644 tests/lac/gmres_reorthogonalize_05/cmp/generic diff --git a/deal.II/include/deal.II/lac/solver_gmres.h b/deal.II/include/deal.II/lac/solver_gmres.h index f60a5b0ff5..b52355ae10 100644 --- a/deal.II/include/deal.II/lac/solver_gmres.h +++ b/deal.II/include/deal.II/lac/solver_gmres.h @@ -34,21 +34,16 @@ DEAL_II_NAMESPACE_OPEN namespace internal { /** - * A namespace for a helper class - * to the GMRES solver. + * A namespace for a helper class to the GMRES solver. */ namespace SolverGMRES { /** - * Class to hold temporary - * vectors. This class - * automatically allocates a new - * vector, once it is needed. + * Class to hold temporary vectors. This class automatically allocates a + * new vector, once it is needed. * - * A future version should also - * be able to shift through - * vectors automatically, - * avoiding restart. + * A future version should also be able to shift through vectors + * automatically, avoiding restart. */ template @@ -56,9 +51,7 @@ namespace internal { public: /** - * Constructor. Prepares an - * array of @p VECTOR of - * length @p max_size. + * Constructor. Prepares an array of @p VECTOR of length @p max_size. */ TmpVectors(const unsigned int max_size, VectorMemory &vmem); @@ -69,21 +62,15 @@ namespace internal ~TmpVectors(); /** - * Get vector number - * @p i. If this vector was - * unused before, an error + * Get vector number @p i. If this vector was unused before, an error * occurs. */ VECTOR &operator[] (const unsigned int i) const; /** - * Get vector number - * @p i. Allocate it if - * necessary. + * Get vector number @p i. Allocate it if necessary. * - * If a vector must be - * allocated, @p temp is - * used to reinit it to the + * If a vector must be allocated, @p temp is used to reinit it to the * proper dimensions. */ VECTOR &operator() (const unsigned int i, @@ -91,22 +78,18 @@ namespace internal private: /** - * Pool were vectors are - * obtained from. + * Pool were vectors are obtained from. */ VectorMemory &mem; /** - * Field for storing the - * vectors. + * Field for storing the vectors. */ std::vector data; /** - * Offset of the first - * vector. This is for later - * when vector rotation will - * be implemented. + * Offset of the first vector. This is for later when vector rotation + * will be implemented. */ unsigned int offset; }; @@ -169,56 +152,49 @@ class SolverGMRES : public Solver { public: /** - * Standardized data struct to - * pipe additional data to the - * solver. + * Standardized data struct to pipe additional data to the solver. */ 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 and - * the residual of the stopping - * criterion to the default residual. + * 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. */ AdditionalData (const unsigned int max_n_tmp_vectors = 30, const bool right_preconditioning = false, - const bool use_default_residual = true); + const bool use_default_residual = true, + const bool force_re_orthogonalization = false); /** - * Maximum number of - * temporary vectors. This - * parameter controls the - * size of the Arnoldi basis, - * which for historical - * reasons is + * Maximum number of temporary vectors. This parameter controls the size + * of the Arnoldi basis, which for historical reasons is * #max_n_tmp_vectors-2. */ unsigned int max_n_tmp_vectors; /** - * Flag for right - * preconditioning. + * Flag for right preconditioning. * - * @note Change between left - * and right preconditioning - * will also change the way - * residuals are - * evaluated. See the - * corresponding section in - * the SolverGMRES. + * @note Change between left and right preconditioning will also change + * the way residuals are evaluated. See the corresponding section in the + * SolverGMRES. */ bool right_preconditioning; /** - * Flag for the default - * residual that is used to - * measure convergence. + * Flag for the default residual that is used to measure convergence. */ bool use_default_residual; + + /** + * Flag to force re-orthogonalization of orthonormal basis in every + * step. If set to false, the solver automatically checks for loss of + * orthogonality every 10 iterations and enables re-orthogonalization only + * if necessary. + */ + bool force_re_orthogonalization; }; /** @@ -229,16 +205,14 @@ public: const AdditionalData &data=AdditionalData()); /** - * Constructor. Use an object of - * type GrowingVectorMemory as - * a default to allocate memory. + * Constructor. Use an object of type GrowingVectorMemory as a default to + * allocate memory. */ SolverGMRES (SolverControl &cn, const AdditionalData &data=AdditionalData()); /** - * Solve the linear system $Ax=b$ - * for x. + * Solve the linear system $Ax=b$ for x. */ template void @@ -255,35 +229,51 @@ public: protected: /** - * Includes the maximum number of - * tmp vectors. + * Includes the maximum number of tmp vectors. */ AdditionalData additional_data; /** - * Implementation of the computation of - * the norm of the residual. + * Implementation of the computation of the norm of the residual. */ virtual double criterion(); /** - * Transformation of an upper - * Hessenberg matrix into - * tridiagonal structure by givens - * rotation of the last column + * Transformation of an upper Hessenberg matrix into tridiagonal structure + * by givens rotation of the last column */ void givens_rotation (Vector &h, Vector &b, Vector &ci, Vector &si, int col) const; + + /** + * Orthogonalize the vector @p vv against the @p dim (orthogonal) vectors + * given by the first argument using the modified Gram-Schmidt + * algorithm. The factors used for orthogonalization are stored in @p h. The + * boolean @p re_orthogonalize specifies whether the modified Gram-Schmidt + * algorithm should be applied twice. The algorithm checks loss of + * orthogonality in the procedure every fifth step and sets the flag to true + * in that case. All subsequent iterations use re-orthogonalization. + */ + static double + modified_gram_schmidt (const internal::SolverGMRES::TmpVectors &orthogonal_vectors, + const unsigned int dim, + const unsigned int accumulated_iterations, + VECTOR &vv, + Vector &h, + bool &re_orthogonalize); + /** * Projected system matrix */ FullMatrix H; + /** * Auxiliary matrix for inverting @p H */ FullMatrix H1; + private: /** * No copy constructor. @@ -315,21 +305,14 @@ class SolverFGMRES : public Solver { public: /** - * Standardized data struct to - * pipe additional data to the - * solver. + * Standardized data struct to pipe additional data to the solver. */ struct AdditionalData { /** - * Constructor. By default, - * set the number of - * temporary vectors to 30, - * preconditioning from left - * and the residual of the - * stopping criterion to the - * default residual - * (cf. class documentation). + * Constructor. By default, set the number of temporary vectors to 30, + * preconditioning from left and the residual of the stopping criterion to + * the default residual (cf. class documentation). */ AdditionalData(const unsigned int max_basis_size = 30, const bool /*use_default_residual*/ = true) @@ -338,8 +321,7 @@ public: {} /** - * Maximum number of - * tmp vectors. + * Maximum number of tmp vectors. */ unsigned int max_basis_size; }; @@ -352,16 +334,14 @@ public: const AdditionalData &data=AdditionalData()); /** - * Constructor. Use an object of - * type GrowingVectorMemory as - * a default to allocate memory. + * Constructor. Use an object of type GrowingVectorMemory as a default to + * allocate memory. */ SolverFGMRES (SolverControl &cn, const AdditionalData &data=AdditionalData()); /** - * Solve the linear system $Ax=b$ - * for x. + * Solve the linear system $Ax=b$ for x. */ template void @@ -453,11 +433,13 @@ inline SolverGMRES::AdditionalData:: AdditionalData (const unsigned int max_n_tmp_vectors, const bool right_preconditioning, - const bool use_default_residual) + const bool use_default_residual, + const bool force_re_orthogonalization) : max_n_tmp_vectors(max_n_tmp_vectors), right_preconditioning(right_preconditioning), - use_default_residual(use_default_residual) + use_default_residual(use_default_residual), + force_re_orthogonalization(force_re_orthogonalization) {} @@ -509,6 +491,65 @@ SolverGMRES::givens_rotation (Vector &h, +template +inline +double +SolverGMRES::modified_gram_schmidt (const internal::SolverGMRES::TmpVectors &orthogonal_vectors, + const unsigned int dim, + const unsigned int accumulated_iterations, + VECTOR &vv, + Vector &h, + bool &re_orthogonalize) +{ + const unsigned int inner_iteration = dim - 1; + + // need initial norm for detection of re-orthogonalization, see below + double norm_vv_start = 0; + if (re_orthogonalize == false && inner_iteration % 5 == 4) + norm_vv_start = vv.l2_norm(); + + // Orthogonalization + for (unsigned int i=0 ; i 10. * norm_vv_start * + std::sqrt(std::numeric_limits::epsilon())) + return norm_vv; + + else + { + re_orthogonalize = true; + deallog << "Re-orthogonalization enabled at step " + << accumulated_iterations << std::endl; + } + } + + if (re_orthogonalize == true) + for (unsigned int i=0 ; i template void @@ -517,12 +558,9 @@ SolverGMRES::solve (const MATRIX &A, const VECTOR &b, const PRECONDITIONER &precondition) { - // this code was written a very - // long time ago by people not - // associated with deal.II. we - // don't make any guarantees to its - // optimality or that it even works - // as expected... + // this code was written a very long time ago by people not associated with + // deal.II. we don't make any guarantees to its optimality or that it even + // works as expected... //TODO:[?] Check, why there are two different start residuals. //TODO:[GK] Make sure the parameter in the constructor means maximum basis size @@ -530,8 +568,7 @@ SolverGMRES::solve (const MATRIX &A, deallog.push("GMRES"); const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors; - // Generate an object where basis - // vectors are stored. + // Generate an object where basis vectors are stored. internal::SolverGMRES::TmpVectors tmp_vectors (n_tmp_vectors, this->memory); // number of the present iteration; this @@ -539,12 +576,10 @@ SolverGMRES::solve (const MATRIX &A, // restart unsigned int accumulated_iterations = 0; - // matrix used for the orthogonalization - // process later + // matrix used for the orthogonalization process later H.reinit(n_tmp_vectors, n_tmp_vectors-1); - // some additional vectors, also used - // in the orthogonalization + // some additional vectors, also used in the orthogonalization dealii::Vector gamma(n_tmp_vectors), ci (n_tmp_vectors-1), @@ -556,17 +591,13 @@ SolverGMRES::solve (const MATRIX &A, SolverControl::State iteration_state = SolverControl::iterate; - // switch to determine whether we want a - // left or a right preconditioner. at - // present, left is default, but both - // ways are implemented + // switch to determine whether we want a left or a right preconditioner. at + // present, left is default, but both ways are implemented const bool left_precondition = !additional_data.right_preconditioning; - // Per default the left - // preconditioned GMRes uses the - // preconditioned residual and the - // right preconditioned GMRes uses - // the unpreconditioned residual as - // stopping criterion. + + // Per default the left preconditioned GMRes uses the preconditioned + // residual and the right preconditioned GMRes uses the unpreconditioned + // residual as stopping criterion. const bool use_default_residual = additional_data.use_default_residual; // define two aliases @@ -589,16 +620,15 @@ SolverGMRES::solve (const MATRIX &A, gamma_ = new dealii::Vector (gamma.size()); } - /////////////////////////////////// - // outer iteration: loop until we - // either reach convergence or the - // maximum number of iterations is - // exceeded. each cycle of this - // loop amounts to one restart + bool re_orthogonalize = additional_data.force_re_orthogonalization; + + /////////////////////////////////////////////////////////////////////////// + // outer iteration: loop until we either reach convergence or the maximum + // number of iterations is exceeded. each cycle of this loop amounts to one + // restart do { - // reset this vector to the - // right size + // reset this vector to the right size h.reinit (n_tmp_vectors-1); if (left_precondition) @@ -615,14 +645,9 @@ SolverGMRES::solve (const MATRIX &A, double rho = v.l2_norm(); - // 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 + // 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) { iteration_state = this->control().check ( @@ -661,13 +686,9 @@ SolverGMRES::solve (const MATRIX &A, 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 + // 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) && @@ -684,44 +705,31 @@ SolverGMRES::solve (const MATRIX &A, precondition.vmult(vv,p); } else - { + { precondition.vmult(p, tmp_vectors[inner_iteration]); A.vmult(vv,p); }; dim = inner_iteration+1; - /* Orthogonalization */ - for (unsigned int i=0 ; i::solve (const MATRIX &A, }; A.vmult(*r,*x_); r->sadd(-1.,1.,b); - // Now *r contains the - // unpreconditioned - // residual!! + // Now *r contains the unpreconditioned residual!! if (left_precondition) { const double res=r->l2_norm(); @@ -775,9 +781,8 @@ SolverGMRES::solve (const MATRIX &A, } } }; - // end of inner iteration. now - // calculate the solution from - // the temporary vectors + // end of inner iteration. now calculate the solution from the temporary + // vectors h.reinit(dim); H1.reinit(dim+1,dim); @@ -798,8 +803,7 @@ SolverGMRES::solve (const MATRIX &A, precondition.vmult(v,p); x.add(1.,v); }; - // end of outer iteration. restart if - // no convergence and the number of + // end of outer iteration. restart if no convergence and the number of // iterations is not exceeded } while (iteration_state == SolverControl::iterate); @@ -813,8 +817,7 @@ SolverGMRES::solve (const MATRIX &A, } deallog.pop(); - // in case of failure: throw - // exception + // in case of failure: throw exception if (this->control().last_check() != SolverControl::success) throw SolverControl::NoConvergence (this->control().last_step(), this->control().last_value()); @@ -827,9 +830,8 @@ template double SolverGMRES::criterion () { - // dummy implementation. this function is - // not needed for the present implementation - // of gmres + // dummy implementation. this function is not needed for the present + // implementation of gmres Assert (false, ExcInternalError()); return 0; } @@ -873,18 +875,15 @@ SolverFGMRES::solve ( const unsigned int basis_size = additional_data.max_basis_size; - // Generate an object where basis - // vectors are stored. + // Generate an object where basis vectors are stored. typename internal::SolverGMRES::TmpVectors v (basis_size, this->memory); typename internal::SolverGMRES::TmpVectors z (basis_size, this->memory); - // number of the present iteration; this - // number is not reset to zero upon a + // number of the present iteration; this number is not reset to zero upon a // restart unsigned int accumulated_iterations = 0; - // matrix used for the orthogonalization - // process later + // matrix used for the orthogonalization process later H.reinit(basis_size+1, basis_size); // Vectors for projected system @@ -953,8 +952,7 @@ SolverFGMRES::solve ( this->memory.free(aux); deallog.pop(); - // in case of failure: throw - // exception + // in case of failure: throw exception if (this->control().last_check() != SolverControl::success) throw SolverControl::NoConvergence (this->control().last_step(), this->control().last_value()); diff --git a/tests/lac/gmres_reorthogonalize_01.cc b/tests/lac/gmres_reorthogonalize_01.cc new file mode 100644 index 0000000000..68fbeab708 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_01.cc @@ -0,0 +1,98 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// tests that GMRES builds an orthonormal basis properly for a few difficult +// test matrices + +#include "../tests.h" +#include +#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); + 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, 1e2*std::numeric_limits::epsilon()); + typename SolverGMRES >::AdditionalData data; + data.max_n_tmp_vectors = 80; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + deallog.pop(); +} + +int main() +{ + std::ofstream logfile("gmres_reorthogonalize_01/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); + test(4); + test(5); + deallog.pop(); + deallog.threshold_double(1.e-4); + deallog.push("float"); + test(0); + test(1); + test(2); + test(3); + deallog.pop(); +} + diff --git a/tests/lac/gmres_reorthogonalize_01/cmp/generic b/tests/lac/gmres_reorthogonalize_01/cmp/generic new file mode 100644 index 0000000000..f1e7ab9648 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_01/cmp/generic @@ -0,0 +1,26 @@ + +DEAL:double:0:GMRES::Starting value 8.00 +DEAL:double:0:GMRES::Convergence step 57 value 0 +DEAL:double:1:GMRES::Starting value 8.00 +DEAL:double:1:GMRES::Re-orthogonalization enabled at step 65 +DEAL:double:1:GMRES::Convergence step 65 value 0 +DEAL:double:2:GMRES::Starting value 8.00 +DEAL:double:2:GMRES::Convergence step 57 value 0 +DEAL:double:3:GMRES::Starting value 8.00 +DEAL:double:3:GMRES::Re-orthogonalization enabled at step 65 +DEAL:double:3:GMRES::Convergence step 65 value 0 +DEAL:double:4:GMRES::Starting value 8.00 +DEAL:double:4:GMRES::Convergence step 57 value 0 +DEAL:double:5:GMRES::Starting value 8.00 +DEAL:double:5:GMRES::Re-orthogonalization enabled at step 65 +DEAL:double:5:GMRES::Convergence step 65 value 0 +DEAL:float:0:GMRES::Starting value 8.00 +DEAL:float:0:GMRES::Convergence step 37 value 0 +DEAL:float:1:GMRES::Starting value 8.00 +DEAL:float:1:GMRES::Re-orthogonalization enabled at step 65 +DEAL:float:1:GMRES::Convergence step 65 value 0 +DEAL:float:2:GMRES::Starting value 8.00 +DEAL:float:2:GMRES::Convergence step 37 value 0 +DEAL:float:3:GMRES::Starting value 8.00 +DEAL:float:3:GMRES::Re-orthogonalization enabled at step 65 +DEAL:float:3:GMRES::Convergence step 66 value 0 diff --git a/tests/lac/gmres_reorthogonalize_02.cc b/tests/lac/gmres_reorthogonalize_02.cc new file mode 100644 index 0000000000..98a1a508cf --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_02.cc @@ -0,0 +1,65 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// tests that GMRES builds an orthonormal basis properly for a larger matrix +// and basis size than gmres_orthogonalize_01 + +#include "../tests.h" +#include +#include +#include +#include + + + +template +void test () +{ + const unsigned int n = 200; + Vector rhs(n), sol(n); + rhs = 1.; + + // only add diagonal entries + SparsityPattern sp(n, n); + sp.compress(); + SparseMatrix matrix(sp); + + for (unsigned int i=0; i::epsilon()); + typename SolverGMRES >::AdditionalData data; + data.max_n_tmp_vectors = 202; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + deallog.pop(); +} + +int main() +{ + std::ofstream logfile("gmres_reorthogonalize_02/output"); + deallog << std::setprecision(4); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + deallog.push("double"); + test(); + deallog.pop(); + deallog.threshold_double(1.e-4); + deallog.push("float"); + test(); + deallog.pop(); +} + diff --git a/tests/lac/gmres_reorthogonalize_02/cmp/generic b/tests/lac/gmres_reorthogonalize_02/cmp/generic new file mode 100644 index 0000000000..561b5e490e --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_02/cmp/generic @@ -0,0 +1,5 @@ + +DEAL:double:GMRES::Starting value 14.14 +DEAL:double:GMRES::Convergence step 109 value 0 +DEAL:float:GMRES::Starting value 14.14 +DEAL:float:GMRES::Convergence step 68 value 0 diff --git a/tests/lac/gmres_reorthogonalize_03.cc b/tests/lac/gmres_reorthogonalize_03.cc new file mode 100644 index 0000000000..1ccf7b25f5 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_03.cc @@ -0,0 +1,98 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// same as gmres_reorthogonalize_01 but forces re-orthogonalization. + +#include "../tests.h" +#include +#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); + 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, 1e2*std::numeric_limits::epsilon()); + typename SolverGMRES >::AdditionalData data; + data.max_n_tmp_vectors = 80; + data.force_re_orthogonalization = true; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); + + deallog.pop(); +} + +int main() +{ + std::ofstream logfile("gmres_reorthogonalize_03/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); + test(4); + test(5); + deallog.pop(); + deallog.threshold_double(1.e-4); + deallog.push("float"); + test(0); + test(1); + test(2); + test(3); + deallog.pop(); +} + diff --git a/tests/lac/gmres_reorthogonalize_03/cmp/generic b/tests/lac/gmres_reorthogonalize_03/cmp/generic new file mode 100644 index 0000000000..54cff0cfe8 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_03/cmp/generic @@ -0,0 +1,21 @@ + +DEAL:double:0:GMRES::Starting value 8.00 +DEAL:double:0:GMRES::Convergence step 57 value 0 +DEAL:double:1:GMRES::Starting value 8.00 +DEAL:double:1:GMRES::Convergence step 64 value 0 +DEAL:double:2:GMRES::Starting value 8.00 +DEAL:double:2:GMRES::Convergence step 57 value 0 +DEAL:double:3:GMRES::Starting value 8.00 +DEAL:double:3:GMRES::Convergence step 64 value 0 +DEAL:double:4:GMRES::Starting value 8.00 +DEAL:double:4:GMRES::Convergence step 57 value 0 +DEAL:double:5:GMRES::Starting value 8.00 +DEAL:double:5:GMRES::Convergence step 64 value 0 +DEAL:float:0:GMRES::Starting value 8.00 +DEAL:float:0:GMRES::Convergence step 37 value 0 +DEAL:float:1:GMRES::Starting value 8.00 +DEAL:float:1:GMRES::Convergence step 64 value 0 +DEAL:float:2:GMRES::Starting value 8.00 +DEAL:float:2:GMRES::Convergence step 37 value 0 +DEAL:float:3:GMRES::Starting value 8.00 +DEAL:float:3:GMRES::Convergence step 64 value 0 diff --git a/tests/lac/gmres_reorthogonalize_04.cc b/tests/lac/gmres_reorthogonalize_04.cc new file mode 100644 index 0000000000..fc06e0d1e4 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_04.cc @@ -0,0 +1,87 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// same as gmres_reorthogonalize_02 but with worse inner product + +#include "../tests.h" +#include +#include +#include +#include + + + +// reimplements things from vector.templates.h in another way which is more +// prone to roundoff, the linker will select this version instead of the one +// in the deal.II library +namespace dealii{ +template +template +Number Vector::operator* (const Vector &v) const +{ + Number sum = 0; + for (unsigned int i=0; i +typename Vector::real_type Vector::l2_norm () const +{ + real_type sum = 0; + for (unsigned int i=0; i +void test () +{ + const unsigned int n = 200; + Vector rhs(n), sol(n); + rhs = 1.; + + // only add diagonal entries + SparsityPattern sp(n, n); + sp.compress(); + SparseMatrix matrix(sp); + + for (unsigned int i=0; i::epsilon()); + typename SolverGMRES >::AdditionalData data; + data.max_n_tmp_vectors = 202; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); +} + +int main() +{ + std::ofstream logfile("gmres_reorthogonalize_04/output"); + deallog << std::setprecision(3); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + deallog.push("double"); + test(); + deallog.pop(); + deallog.threshold_double(1.e-4); + deallog.push("float"); + test(); + deallog.pop(); +} + diff --git a/tests/lac/gmres_reorthogonalize_04/cmp/generic b/tests/lac/gmres_reorthogonalize_04/cmp/generic new file mode 100644 index 0000000000..eb5639ce73 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_04/cmp/generic @@ -0,0 +1,5 @@ + +DEAL:double:GMRES::Starting value 14.1 +DEAL:double:GMRES::Convergence step 202 value 0 +DEAL:float:GMRES::Starting value 14.1 +DEAL:float:GMRES::Convergence step 201 value 0 diff --git a/tests/lac/gmres_reorthogonalize_05.cc b/tests/lac/gmres_reorthogonalize_05.cc new file mode 100644 index 0000000000..9d8e2ad10a --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_05.cc @@ -0,0 +1,88 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// same as gmres_reorthogonalize_04 but with forced orthogonalization + +#include "../tests.h" +#include +#include +#include +#include + + + +// reimplements things from vector.templates.h in another way which is more +// prone to roundoff, the linker will select this version instead of the one +// in the deal.II library +namespace dealii{ +template +template +Number Vector::operator* (const Vector &v) const +{ + Number sum = 0; + for (unsigned int i=0; i +typename Vector::real_type Vector::l2_norm () const +{ + real_type sum = 0; + for (unsigned int i=0; i +void test () +{ + const unsigned int n = 200; + Vector rhs(n), sol(n); + rhs = 1.; + + // only add diagonal entries + SparsityPattern sp(n, n); + sp.compress(); + SparseMatrix matrix(sp); + + for (unsigned int i=0; i::epsilon()); + typename SolverGMRES >::AdditionalData data; + data.max_n_tmp_vectors = 202; + data.force_re_orthogonalization = true; + + SolverGMRES > solver(control, data); + solver.solve(matrix, sol, rhs, PreconditionIdentity()); +} + +int main() +{ + std::ofstream logfile("gmres_reorthogonalize_05/output"); + deallog << std::setprecision(3); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + deallog.push("double"); + test(); + deallog.pop(); + deallog.threshold_double(1.e-4); + deallog.push("float"); + test(); + deallog.pop(); +} + diff --git a/tests/lac/gmres_reorthogonalize_05/cmp/generic b/tests/lac/gmres_reorthogonalize_05/cmp/generic new file mode 100644 index 0000000000..2e79ab3e16 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_05/cmp/generic @@ -0,0 +1,5 @@ + +DEAL:double:GMRES::Starting value 14.1 +DEAL:double:GMRES::Convergence step 109 value 0 +DEAL:float:GMRES::Starting value 14.1 +DEAL:float:GMRES::Convergence step 66 value 0 -- 2.39.5