From: guido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 29 May 2002 10:58:56 +0000 (+0000)
Subject: vector initialization in 1.47 fixed
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=61216bb447ebb15a9c60daeb7fba247b14b61e4d;p=dealii-svn.git

vector initialization in 1.47 fixed


git-svn-id: https://svn.dealii.org/trunk@5924 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/lac/include/lac/solver_gmres.h b/deal.II/lac/include/lac/solver_gmres.h
index e20129f2bc..9509af21e2 100644
--- a/deal.II/lac/include/lac/solver_gmres.h
+++ b/deal.II/lac/include/lac/solver_gmres.h
@@ -158,25 +158,142 @@ class SolverGMRES : public Solver<VECTOR>
     void givens_rotation (Vector<double> &h,  Vector<double> &b,
 			  Vector<double> &ci, Vector<double> &si,
 			  int col) const;
-				   /**
-				    * Projected system matrix
-				    */
-  FullMatrix<double> H;
-				   /**
-				    * Auxiliary matrix for inverting @p{H}
-				    */
-  FullMatrix<double> H1;
-
+				     /**
+				      * Projected system matrix
+				      */
+    FullMatrix<double> H;
+				     /**
+				      * Auxiliary matrix for inverting @p{H}
+				      */
+    FullMatrix<double> H1;
+    
   private:
-				   /**
-				    * No copy constructor.
-				    */
-  SolverGMRES (const SolverGMRES<VECTOR>&);
+				     /**
+				      * 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.
+				      */
+    
+    class TmpVectors
+    {
+      public:
+					 /**
+					  * Constructor. Prepares an
+					  * array of @p{VECTOR} of
+					  * length @p{max_size}.
+					  */
+	TmpVectors(unsigned int max_size,
+		   VectorMemory<VECTOR>& vmem);
+					 /**
+					  * Delete all allocated vectors.
+					  */
+	~TmpVectors();
+
+					 /**
+					  * Get vector number
+					  * @p{i}. If this vector was
+					  * unused before, an error
+					  * occurs.
+					  */
+	VECTOR& operator[] (unsigned int i) const;
+	
+					 /**
+					  * Get vector number
+					  * @p{i}. Allocate it if
+					  * necessary.
+					  *
+					  * If a vector must be
+					  * allocated, @p{temp} is
+					  * used to reinit it to the
+					  * proper dimensions.
+					  */
+	VECTOR& operator() (unsigned int i,
+			    const VECTOR& temp);
+	
+      private:
+					 /**
+					  * Pool were vectors are obtained from.
+					  */
+	VectorMemory<VECTOR>& mem;
+					 /**
+					  * Field for storing the vectors.
+					  */
+	std::vector<VECTOR*> data;
+					 /**
+					  * Offset of the first
+					  * vector. This is for later
+					  * when vector rotation will
+					  * be implemented.
+					  */
+	unsigned int offset;
+    };
+				     /**
+				      * No copy constructor.
+				      */
+    SolverGMRES (const SolverGMRES<VECTOR>&);
 };
 
 
 /* --------------------- Inline and template functions ------------------- */
 
+template <class VECTOR>
+inline
+SolverGMRES<VECTOR>::TmpVectors::TmpVectors (unsigned int max_size,
+				 VectorMemory<VECTOR>& vmem)
+		:
+		mem(vmem),
+		data (max_size, 0),
+		offset(0)
+{}
+
+
+template <class VECTOR>
+inline
+SolverGMRES<VECTOR>::TmpVectors::~TmpVectors ()
+{
+  for (typename std::vector<VECTOR*>::iterator v = data.begin();
+       v != data.end(); ++v)
+    if (*v != 0)
+      mem.free(*v);
+}
+
+
+template <class VECTOR>
+inline VECTOR&
+SolverGMRES<VECTOR>::TmpVectors::operator[] (unsigned int i) const
+{
+  Assert (i+offset<data.size(),
+	  ExcIndexRange(i,-offset, data.size()-offset));
+  
+  i -= offset;
+
+  Assert (data[i] != 0, ExcNotInitialized());
+  return *data[i];
+}
+
+
+template <class VECTOR>
+inline VECTOR&
+SolverGMRES<VECTOR>::TmpVectors::operator() (unsigned int i,
+					     const VECTOR& temp)
+{
+  Assert (i+offset<data.size(),
+	  ExcIndexRange(i,-offset, data.size()-offset));
+  i -= offset;
+  if (data[i] == 0)
+    {
+      data[i] = mem.alloc();
+      data[i]->reinit(temp);
+    }
+  return *data[i];
+}
+
 
 template <class VECTOR>
 SolverGMRES<VECTOR>::SolverGMRES (SolverControl        &cn,
@@ -238,22 +355,10 @@ SolverGMRES<VECTOR>::solve (const MATRIX         &A,
   deallog.push("GMRES");
   const unsigned int n_tmp_vectors = additional_data.max_n_tmp_vectors;
 
-				   // allocate an array of n_tmp_vectors
-				   // temporary vectors from the VectorMemory
-				   // object; resize them but do not set their
-				   // values since they will be overwritten soon
-				   // anyway.
-
-				   // This is really bad since vectors
-				   // should only be allocated if
-				   // really needed. (GK)
-  std::vector<VECTOR*> tmp_vectors (n_tmp_vectors, 0);
-  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
-    {
-      tmp_vectors[tmp] = memory.alloc();
-      tmp_vectors[tmp]->reinit (x, true);
-    };
-
+				   // Generate an object where basis
+				   // vectors are stored.
+  TmpVectors tmp_vectors (n_tmp_vectors, memory);
+  
 				   // number of the present iteration; this
 				   // number is not reset to zero upon a
 				   // restart
@@ -282,8 +387,8 @@ unsigned int dim = 0;
   const bool left_precondition = ! additional_data.right_preconditioning;
 
 				   // define two aliases
-  VECTOR &v = *tmp_vectors[0];
-  VECTOR &p = *tmp_vectors[n_tmp_vectors-1];
+  VECTOR &v = tmp_vectors(0, x);
+  VECTOR &p = tmp_vectors(n_tmp_vectors-1, x);
 
 
 ///////////////////////////////////
@@ -339,14 +444,14 @@ unsigned int dim = 0;
 	{
 	  ++accumulated_iterations;
 	  // yet another alias
-	  VECTOR& vv = *tmp_vectors[inner_iteration+1];
+	  VECTOR& vv = tmp_vectors(inner_iteration+1, x);
 	  
 	  if (left_precondition)
 	    {
-	      A.vmult(p, *tmp_vectors[inner_iteration]);
+	      A.vmult(p, tmp_vectors[inner_iteration]);
 	      precondition.vmult(vv,p);
 	    } else {
-	      precondition.vmult(p,*tmp_vectors[inner_iteration]);
+	      precondition.vmult(p, tmp_vectors[inner_iteration]);
 	      A.vmult(vv,p);
 	    };
 	  
@@ -355,8 +460,8 @@ unsigned int dim = 0;
 					   /* Orthogonalization */
 	  for (unsigned int i=0 ; i<dim ; ++i)
 	    {
-	      h(i) = vv * *tmp_vectors[i];
-	      vv.add(-h(i),*tmp_vectors[i]);
+	      h(i) = vv * tmp_vectors[i];
+	      vv.add(-h(i), tmp_vectors[i]);
 	    };
       
 	  double s = vv.l2_norm();
@@ -365,9 +470,9 @@ unsigned int dim = 0;
 					   /* Re-orthogonalization */
 	  for (unsigned i=0; i<dim; ++i)
 	    {
-	      double htmp = vv * *tmp_vectors[i];
+	      double htmp = vv * tmp_vectors[i];
 	      h(i) += htmp;
-	      vv.add(-htmp,*tmp_vectors[i]);
+	      vv.add(-htmp, tmp_vectors[i]);
 	    };
 	  
 	  s = vv.l2_norm();
@@ -403,12 +508,12 @@ unsigned int dim = 0;
 
       if (left_precondition)
 	for (unsigned int i=0 ; i<dim; ++i)
-	  x.add(h(i), *tmp_vectors[i]);
+	  x.add(h(i), tmp_vectors[i]);
       else
 	{
 	  p = 0.;
 	  for (unsigned int i=0; i<dim; ++i)
-	    p.add(h(i), *tmp_vectors[i]);
+	    p.add(h(i), tmp_vectors[i]);
 	  precondition.vmult(v,p);
 	  x.add(1.,v);
 	};
@@ -418,13 +523,8 @@ unsigned int dim = 0;
     }
   while (iteration_state == SolverControl::iterate);
 
-				   // free the allocated memory before
-				   // leaving
-  for (unsigned int tmp=0; tmp<n_tmp_vectors; ++tmp)
-    memory.free (tmp_vectors[tmp]);
-
+  
   deallog.pop();
-
 				   // in case of failure: throw
 				   // exception
   if (control().last_check() != SolverControl::success)