From 4265baf335bcf21477a0afbe300b492fd181c2f3 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 31 Aug 2009 13:02:30 +0000 Subject: [PATCH] Make the initialize function in Chebyshev preconditioner unique. git-svn-id: https://svn.dealii.org/trunk@19351 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/precondition.h | 84 ++++++++++---------------- deal.II/lac/include/lac/solver_cg.h | 4 +- 2 files changed, 33 insertions(+), 55 deletions(-) diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h index 5103538aba..dd97d3b064 100644 --- a/deal.II/lac/include/lac/precondition.h +++ b/deal.II/lac/include/lac/precondition.h @@ -916,6 +916,12 @@ public: * eigenvalue. */ double eig_cg_residual; + + /** + * Stores the inverse of the diagonal + * of the underlying matrix. + */ + VECTOR matrix_diagonal_inverse; }; PreconditionChebyshev (); @@ -929,9 +935,9 @@ public: * matrix has an operator * el(i,i) for accessing all * the elements in the - * diagonal. Otherwise, use the other - * initialize function to - * manually provide the diagonal. + * diagonal. Alternatively, the + * diagonal can be supplied with the + * help of the AdditionalData field. * * This function calculates an * estimate of the eigenvalue range @@ -940,19 +946,6 @@ public: * iteration. */ void initialize (const MATRIX &matrix, - const AdditionalData &additional_data = AdditionalData()); - - /** - * Second initialize function. Takes - * the matrix which is used to form - * the preconditioner, a vector - * containing the inverse of - * the diagonal of the input matrix, - * and additional flags if there are - * any. - */ - void initialize (const MATRIX &matrix, - const VECTOR &diagonal_inverse, const AdditionalData &additional_data = AdditionalData()); /** @@ -986,12 +979,6 @@ private: */ SmartPointer matrix_ptr; - /** - * Stores the inverse of the diagonal - * of the underlying matrix. - */ - VECTOR diagonal_inverse; - /** * Internal vector used for * vmult operation. @@ -1544,28 +1531,19 @@ inline void PreconditionChebyshev::initialize (const MATRIX &matrix, const AdditionalData &additional_data) -{ - VECTOR diagonal_inv (matrix.m()); - for (unsigned int i=0; i -inline -void -PreconditionChebyshev::initialize (const MATRIX &matrix, - const VECTOR &diagonal_inv, - const AdditionalData &additional_data) { Assert (matrix.m() == matrix.n(), ExcMessage("Matrix not quadratic.")); matrix_ptr = &matrix; - update1.reinit (diagonal_inv, true); - update2.reinit (diagonal_inv, true); - diagonal_inverse = diagonal_inv; data = additional_data; + if (data.matrix_diagonal_inverse.size() != matrix.m()) + { + data.matrix_diagonal_inverse.reinit(matrix.m()); + for (unsigned int i=0; i::initialize (const MATRIX &matrix, double max_eigenvalue, min_eigenvalue; { double eigen_beta_alpha = 0; - + std::vector diagonal; std::vector offdiagonal; - + VECTOR rhs, g; - rhs.reinit(diagonal_inv, true); + rhs.reinit(data.matrix_diagonal_inverse, true); rhs = 1./sqrt(matrix.m()); - g.reinit(diagonal_inv, true); + g.reinit(data.matrix_diagonal_inverse, true); - int it=0; + unsigned int it=0; double res,gh,alpha,beta; g.equ(-1.,rhs); @@ -1599,7 +1577,7 @@ PreconditionChebyshev::initialize (const MATRIX &matrix, { it++; matrix.vmult (update1, update2); - update1.scale (diagonal_inverse); + update1.scale (data.matrix_diagonal_inverse); alpha = update2 * update1; alpha = gh/alpha; g.add (alpha, update1); @@ -1654,13 +1632,13 @@ PreconditionChebyshev::vmult (VECTOR &dst, matrix_ptr->vmult (update1, dst); update1 -= src; update1 /= theta; - update1.scale (diagonal_inverse); + update1.scale (data.matrix_diagonal_inverse); dst -= update1; } else { dst.equ (1./theta, src); - dst.scale (diagonal_inverse); + dst.scale (data.matrix_diagonal_inverse); update1.equ(-1.,dst); } @@ -1668,7 +1646,7 @@ PreconditionChebyshev::vmult (VECTOR &dst, { matrix_ptr->vmult (update2, dst); update2 -= src; - update2.scale (diagonal_inverse); + update2.scale (data.matrix_diagonal_inverse); const double rhokp = 1./(2.*sigma-rhok); const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta; rhok = rhokp; @@ -1692,13 +1670,13 @@ PreconditionChebyshev::Tvmult (VECTOR &dst, matrix_ptr->Tvmult (update1, dst); update1 -= src; update1 /= theta; - update1.scale (diagonal_inverse); + update1.scale (data.matrix_diagonal_inverse); dst -= update1; } else { dst.equ (1./theta, src); - dst.scale (diagonal_inverse); + dst.scale (data.matrix_diagonal_inverse); update1.equ(-1.,dst); } @@ -1706,7 +1684,7 @@ PreconditionChebyshev::Tvmult (VECTOR &dst, { matrix_ptr->Tvmult (update2, dst); update2 -= src; - update2.scale (diagonal_inverse); + update2.scale (data.matrix_diagonal_inverse); const double rhokp = 1./(2.*sigma-rhok); const double factor1 = rhokp * rhok, factor2 = 2.*rhokp/delta; rhok = rhokp; @@ -1723,7 +1701,7 @@ void PreconditionChebyshev::clear () { is_initialized = false; matrix_ptr = 0; - diagonal_inverse.reinit(0); + data.matrix_diagonal_inverse.reinit(0); update1.reinit(0); update2.reinit(0); } diff --git a/deal.II/lac/include/lac/solver_cg.h b/deal.II/lac/include/lac/solver_cg.h index 8833b36b39..81be4a62a7 100644 --- a/deal.II/lac/include/lac/solver_cg.h +++ b/deal.II/lac/include/lac/solver_cg.h @@ -347,7 +347,7 @@ SolverCG::solve (const MATRIX &A, else g.equ(-1.,b); res = g.l2_norm(); - + conv = this->control().check(0,res); if (conv) { @@ -378,7 +378,7 @@ SolverCG::solve (const MATRIX &A, conv = this->control().check(it,res); if (conv) break; - + precondition.vmult(h,g); beta = gh; -- 2.39.5