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 e20129f2bc2267a8911e533902b2234e8d006070..9509af21e23ea239602940ab83328de65dfa8656 100644 (file)
--- 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)