From 4265baf335bcf21477a0afbe300b492fd181c2f3 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
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
                                         * <tt>el(i,i)</tt> for accessing all
                                         * the elements in the
-                                        * diagonal. Otherwise, use the other
-                                        * <tt>initialize</tt> 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 <i>inverse</i> 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<const MATRIX> matrix_ptr;
 
-				       /**
-					* Stores the inverse of the diagonal
-					* of the underlying matrix.
-					*/
-  VECTOR diagonal_inverse;
-
 				       /**
 					* Internal vector used for
 					* <tt>vmult</tt> operation.
@@ -1544,28 +1531,19 @@ inline
 void 
 PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
 						  const AdditionalData &additional_data)
-{
-  VECTOR diagonal_inv (matrix.m());
-  for (unsigned int i=0; i<matrix.m(); ++i)
-    diagonal_inv(i) = 1./matrix.el(i,i);
-
-  initialize (matrix, diagonal_inv, additional_data);
-}
-
-
-template <class MATRIX, class VECTOR>
-inline
-void 
-PreconditionChebyshev<MATRIX,VECTOR>::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<matrix.m(); ++i)
+	data.matrix_diagonal_inverse(i) = 1./matrix.el(i,i);
+    }
+  update1.reinit (data.matrix_diagonal_inverse, true);
+  update2.reinit (data.matrix_diagonal_inverse, true);
+
 
 				 // calculate largest eigenvalue using a
 				 // hand-tuned CG iteration on the matrix
@@ -1578,16 +1556,16 @@ PreconditionChebyshev<MATRIX,VECTOR>::initialize (const MATRIX &matrix,
   double max_eigenvalue, min_eigenvalue;
   {
     double eigen_beta_alpha = 0;
-  
+
     std::vector<double> diagonal;
     std::vector<double> 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<MATRIX,VECTOR>::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<MATRIX,VECTOR>::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<MATRIX,VECTOR>::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<MATRIX,VECTOR>::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<MATRIX,VECTOR>::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<MATRIX,VECTOR>::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<VECTOR>::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<VECTOR>::solve (const MATRIX         &A,
 	conv = this->control().check(it,res);
 	if (conv)
 	  break;
-	
+
 	precondition.vmult(h,g);
 	
 	beta = gh;
-- 
2.39.5