]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement a scheme to free unused memory for dVector and dSMatrix. The memory is...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 3 Jun 1998 00:36:10 +0000 (00:36 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 3 Jun 1998 00:36:10 +0000 (00:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@373 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/dsmatrix.h
deal.II/lac/include/lac/dvector.h
deal.II/lac/source/dsmatrix.cc
deal.II/lac/source/dvector.cc

index 064cf26a48a504ff953d340d9dc1fcac079f7ac9..9291ff083399b7ffb9906629d4168a7aa59884e5 100644 (file)
@@ -28,13 +28,21 @@ template <int dim> class DoFHandler;
 CLASS
    dSMatrixStruct
 
-   @author Original version by Roland Becker, Guido Kanschat,
-      Franz-Theo Suttmeier; lots of enhancements, reorganisation
-      and documentation by Wolfgang Bangerth
+   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
    */
 class dSMatrixStruct
 {
   public:
+                                    /**
+                                     * Initialize the matrix empty, i.e. with
+                                     * no memory allocated. This is useful if
+                                     * you want such objects as member
+                                     * variables in other classes. You can make
+                                     * the structure usable by calling the
+                                     * #reinit# function.
+                                     */
+    dSMatrixStruct ();
+    
                                     /**
                                      * Initialize a rectangular matrix with
                                      * #m# rows and #n# columns,
@@ -50,7 +58,8 @@ class dSMatrixStruct
                                      * #n# with at most #max_per_row#
                                      * nonzero entries per row.
                                      */
-    dSMatrixStruct (const unsigned int n, const unsigned int max_per_row);
+    dSMatrixStruct (const unsigned int n,
+                   const unsigned int max_per_row);
 
                                     /**
                                      * Destructor.
@@ -63,6 +72,14 @@ class dSMatrixStruct
                                      * #m# rows and #n# columns,
                                      * with at most #max_per_row#
                                      * nonzero entries per row.
+                                     *
+                                     * If #m*n==0# all memory is freed,
+                                     * resulting in a total reinitialization
+                                     * of the object. If it is nonzero, new
+                                     * memory is only allocated if the new
+                                     * size extends the old one. This is done
+                                     * to save time and to avoid fragmentation
+                                     * of the heap.
                                      */
     void reinit (const unsigned int m,
                 const unsigned int n,
@@ -86,6 +103,15 @@ class dSMatrixStruct
                                      */
     void compress ();
 
+                                    /**
+                                     * Return whether the object is empty. It
+                                     * is empty if no memory is allocated,
+                                     * which is the same as that both
+                                     * dimensions are zero.
+                                     */
+    bool empty () const;
+    
+
                                     /**
                                      * Return the index of the matrix
                                      * element with row number #i# and
@@ -289,6 +315,10 @@ class dSMatrixStruct
                                     /**
                                      * Exception
                                      */
+    DeclException0 (ExcEmptyObject);
+                                    /**
+                                     * Exception
+                                     */
     DeclException0 (ExcInternalError);
 
   private:
@@ -351,7 +381,7 @@ class dSMatrix
                                      * long as #reinit# is not called with a
                                      * new sparsity structure.
                                      */
-    dSMatrix (dSMatrixStruct& c);
+    dSMatrix (const dSMatrixStruct &sparsity);
     
                                     /**
                                      * Destructor. Free all memory, but do not
@@ -361,15 +391,31 @@ class dSMatrix
     ~dSMatrix ();
     
 
-                                    //
+                                    /**
+                                     * Reinitialize the object but keep to
+                                     * the sparsity pattern previously used.
+                                     * This may be necessary if you #reinit#'d
+                                     * the sparsity structure and want to
+                                     * update the size of the matrix.
+                                     *
+                                     * Note that memory is only reallocated if
+                                     * the new size exceeds the old size. If
+                                     * that is not the case, the allocated
+                                     * memory is not reduced. However, if the
+                                     * sparsity structure is empty (i.e. the
+                                     * dimensions are zero), then all memory
+                                     * is freed.
+                                     */
     void reinit ();
-                                    //
-    void reinit (dSMatrixStruct &);
 
                                     /**
-                                     * Release all memory and return to a state
-                                     * just like after having called the
-                                     * default constructor.
+                                     * Reinitialize the sparse matrix with the
+                                     * given sparsity pattern. The latter tells
+                                     * the matrix how many nonzero elements
+                                     * there need to be reserved.
+                                     *
+                                     * Regarding memory allocation, the same
+                                     * applies as said above.
                                      *
                                      * You have to make sure that the lifetime
                                      * of the sparsity structure is at least
@@ -377,6 +423,15 @@ class dSMatrix
                                      * long as #reinit# is not called with a
                                      * new sparsity structure.
                                      */
+    void reinit (const dSMatrixStruct &sparsity);
+
+                                    /**
+                                     * Release all memory and return to a state
+                                     * just like after having called the
+                                     * default constructor. It also forgets the
+                                     * sparsity pattern it was previously tied
+                                     * to.
+                                     */
     void clear ();
     
                                     /**
index e306f96df45baa2c8b5206cae6ce889432d6e87b..a1286f3201cd5bca2c112b46112bb070e59e73ab 100644 (file)
@@ -108,6 +108,9 @@ class dVector : public VectorBase
                                     /**
                                      * Set all entries to zero. Equivalent to
                                      * #v = 0#, but more obvious and faster.
+                                     * Note that this function does not change
+                                     * the size of the vector, unlike the
+                                     * STL's #vector<>::clear# function.
                                      */
     void clear ();
     
@@ -160,10 +163,16 @@ class dVector : public VectorBase
                                     /**
                                      * Change the dimension of the vector to
                                      * #N#. The reserved memory for this vector
-                                     * remains unchanged, however, to make
+                                     * remains unchanged if possible, to make
                                      * things faster, but this may waste some
                                      * memory, so take this in the back of your
                                      * head.
+                                     * However, if #N==0# all memory is freed,
+                                     * i.e. if you want to resize the vector
+                                     * and release the memory not needed, you
+                                     * have to first call #reinit(0)# and then
+                                     * #reinit(N)#. This cited behaviour is
+                                     * analogous to that of the STL containers.
                                      *
                                      * On #fast==false#, the vector is filled by
                                      * zeros.
@@ -175,7 +184,9 @@ class dVector : public VectorBase
                                      * vector #V#. The same applies as for
                                      * the other #reinit# function.
                                      *
-                                     * The elements of #V# are not copied.
+                                     * The elements of #V# are not copied, i.e.
+                                     * this function is the same as calling
+                                     * #reinit (V.size(), fast)#.
                                      */
     void reinit (const dVector& V, const bool fast=false);
     
@@ -457,6 +468,10 @@ class dVector : public VectorBase
                                      * Exception
                                      */
     DeclException0 (ExcOutOfMemory);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcEmptyVector);
 };
 
 
index 95c0525077e2a56c5b079731e9179cd9ae82e044..f1109d6bf160b2dda3a15b397bd1bc446e54e5cf 100644 (file)
 
 
 
+dSMatrixStruct::dSMatrixStruct () :
+               max_dim(0),
+               max_vec_len(0),
+               rowstart(0),
+               colnums(0)
+{
+  reinit (0,0,0);
+};
+
+
+
 dSMatrixStruct::dSMatrixStruct (const unsigned int m, const unsigned int n,
                                const unsigned int max_per_row) 
                : max_dim(0),
@@ -37,10 +48,8 @@ dSMatrixStruct::dSMatrixStruct (const unsigned int n,
 
 dSMatrixStruct::~dSMatrixStruct ()
 {
-  if (rowstart != 0)
-    delete[] rowstart;
-  if (colnums != 0)
-    delete[] colnums;
+  if (rowstart != 0)  delete[] rowstart;
+  if (colnums != 0)   delete[] colnums;
 }
 
 
@@ -50,28 +59,37 @@ void
 dSMatrixStruct::reinit (const unsigned int m, const unsigned int n,
                        const unsigned int max_per_row)
 {
-  Assert (m>0, ExcInvalidNumber(m));
-  Assert (n>0, ExcInvalidNumber(n));
-  Assert (max_per_row>0, ExcInvalidNumber(max_per_row));
+  Assert ((max_per_row>0) || ((m==0) && (n==0)), ExcInvalidNumber(max_per_row));
   rows = m;
   cols = n;
   vec_len = m * max_per_row;
   max_row_len = max_per_row;
 
+  if (m*n == 0)
+    {
+      if (rowstart)  delete[] rowstart;
+      if (colnums)   delete[] colnums;
+      rowstart = 0;
+      colnums = 0;
+      max_vec_len = vec_len = max_dim = rows = cols = 0;
+      compressed = false;
+      return;
+    };
+  
   if (rows > max_dim)
     {
       if (rowstart) delete[] rowstart;
       max_dim = rows;
       rowstart = new unsigned int[max_dim+1];
-    }
-
+    };
+  
   if (vec_len > max_vec_len)
     {
       if (colnums) delete[] colnums;
       max_vec_len = vec_len;
       colnums = new int[max_vec_len];
-    }
-
+    };
+  
   for (unsigned int i=0; i<=rows; i++)
     rowstart[i] = i * max_per_row;
   fill_n (&colnums[0], vec_len, -1);
@@ -87,6 +105,8 @@ dSMatrixStruct::reinit (const unsigned int m, const unsigned int n,
 void
 dSMatrixStruct::compress ()
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
+  
   if (compressed) return;
   unsigned int next_free_entry = 0,
                next_row_start = 0,
@@ -160,9 +180,38 @@ dSMatrixStruct::compress ()
 
 
 
+bool
+dSMatrixStruct::empty () const {
+                                  // let's try to be on the safe side of
+                                  // life by using multiple possibilities in
+                                  // the check for emptiness... (sorry for
+                                  // this kludge -- emptying matrices and
+                                  // freeing memory was not present in the
+                                  // original implementation and I don't
+                                  // know at how many places I missed
+                                  // something in adding it, so I try to
+                                  // be cautious. wb)
+  if ((rowstart==0) || (rows==0) || (cols==0))
+    {
+      Assert (rowstart==0, ExcInternalError());
+      Assert (rows==0, ExcInternalError());
+      Assert (cols==0, ExcInternalError());
+      Assert (colnums==0, ExcInternalError());
+      Assert (vec_len==0, ExcInternalError());
+      Assert (max_vec_len==0, ExcInternalError());
+      Assert (vec_len==0, ExcInternalError());
+
+      return true;
+    };
+  return false;
+};
+
+
+
 int
 dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (i<rows, ExcInvalidIndex(i,rows));
   Assert (j<cols, ExcInvalidIndex(j,cols));
   Assert (compressed, ExcNotCompressed());
@@ -196,6 +245,7 @@ dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const
 void
 dSMatrixStruct::add (const unsigned int i, const unsigned int j)
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (i<rows, ExcInvalidIndex(i,rows));
   Assert (j<cols, ExcInvalidIndex(j,cols));
   Assert (compressed==false, ExcMatrixIsCompressed());
@@ -224,6 +274,7 @@ dSMatrixStruct::add (const unsigned int i, const unsigned int j)
 void
 dSMatrixStruct::add_matrix (const unsigned int n, const int* rowcols)
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   for (unsigned int i=0; i<n; ++i)
     for (unsigned int j=0; j<n; ++j)
       add(rowcols[i], rowcols[j]);
@@ -235,6 +286,7 @@ void
 dSMatrixStruct::add_matrix (const unsigned int m, const unsigned int n,
                            const int* rows, const int* cols)
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   for (unsigned i=0; i<m; ++i)
     for (unsigned j=0; j<n; ++j)
       add(rows[i], cols[j]);
@@ -245,6 +297,7 @@ dSMatrixStruct::add_matrix (const unsigned int m, const unsigned int n,
 void
 dSMatrixStruct::add_matrix (const iVector& rowcols)
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   unsigned int i,j;
   for (i=0;i<rowcols.n();i++)
     for (j=0;j<rowcols.n();j++)
@@ -256,6 +309,7 @@ dSMatrixStruct::add_matrix (const iVector& rowcols)
 void
 dSMatrixStruct::add_matrix (const iVector& rows, const iVector& cols)
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   unsigned int i,j;
   for (i=0;i<rows.n();i++)
     for (j=0;j<cols.n();j++)
@@ -267,6 +321,7 @@ dSMatrixStruct::add_matrix (const iVector& rows, const iVector& cols)
 void
 dSMatrixStruct::print_gnuplot (ostream &out) const
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
       if (colnums[j]>=0)
@@ -278,6 +333,7 @@ dSMatrixStruct::print_gnuplot (ostream &out) const
 unsigned int
 dSMatrixStruct::bandwidth () const
 {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   unsigned int b=0;
   for (unsigned int i=0; i<rows; ++i)
     for (unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j)
@@ -297,6 +353,7 @@ dSMatrixStruct::bandwidth () const
 
 unsigned int
 dSMatrixStruct::n_nonzero_elements () const {
+  Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (compressed, ExcNotCompressed());
   return colnums[rows]-colnums[0];
 };
@@ -315,7 +372,7 @@ dSMatrix::dSMatrix () :
 
 
 
-dSMatrix::dSMatrix (dSMatrixStruct &c)
+dSMatrix::dSMatrix (const dSMatrixStruct &c)
                : cols(&c), val(0), max_len(0)
 {
   reinit();
@@ -334,22 +391,31 @@ void
 dSMatrix::reinit ()
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
-  Assert (cols->compressed, ExcNotCompressed());
-  
+  Assert (cols->compressed || cols->empty(), ExcNotCompressed());
+
+  if (cols->empty()) 
+    {
+      if (val) delete[] val;
+      val = 0;
+      max_len = 0;
+      return;
+    };
+        
   if (max_len<cols->vec_len)
-  {
-    if (val) delete[] val;
-    val = new double[cols->vec_len];
-    max_len = cols->vec_len;
-  }
-//  memset(val, 0, sizeof(*val) * cols->vec_len);
-  for (int i = cols->vec_len-1 ; i>=0 ; i--) val[i] = 0;
+    {
+      if (val) delete[] val;
+      val = new double[cols->vec_len];
+      max_len = cols->vec_len;
+    };
+
+  if (val)
+    fill_n (&val[0], cols->vec_len, 0);
 }
 
 
 
 void
-dSMatrix::reinit (dSMatrixStruct &sparsity) {
+dSMatrix::reinit (const dSMatrixStruct &sparsity) {
   cols = &sparsity;
   reinit ();
 };
@@ -377,6 +443,7 @@ dSMatrix::n_nonzero_elements () const {
 dSMatrix &
 dSMatrix::copy_from (const dSMatrix &matrix) {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
 
   for (unsigned int i = 0 ; i<cols->vec_len; ++i)
@@ -390,6 +457,7 @@ dSMatrix::copy_from (const dSMatrix &matrix) {
 void
 dSMatrix::add_scaled (const double factor, const dSMatrix &matrix) {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
 
   for (unsigned int i = 0 ; i<cols->vec_len; ++i)
@@ -402,6 +470,7 @@ void
 dSMatrix::vmult (dVector& dst, const dVector& src) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size()));
 
@@ -422,6 +491,7 @@ dSMatrix::vmult (dVector& dst, const dVector& src) const
 void
 dSMatrix::Tvmult (dVector& dst, const dVector& src) const
 {
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert(n() == dst.size(), ExcDimensionsDontMatch(n(),dst.size()));
   Assert(m() == src.size(), ExcDimensionsDontMatch(m(),src.size()));
@@ -444,6 +514,7 @@ double
 dSMatrix::matrix_norm (const dVector& v) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert(m() == v.size(), ExcDimensionsDontMatch(m(),v.size()));
   Assert(n() == v.size(), ExcDimensionsDontMatch(n(),v.size()));
 
@@ -470,6 +541,7 @@ double
 dSMatrix::residual (dVector& dst, const dVector& u, const dVector& b)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
   Assert(m() == b.size(), ExcDimensionsDontMatch(m(),b.size()));
   Assert(n() == u.size(), ExcDimensionsDontMatch(n(),u.size()));
@@ -494,6 +566,7 @@ void
 dSMatrix::Jacobi_precond (dVector& dst, const dVector& src, const double om)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
 
   const unsigned int n = src.size();
 
@@ -507,6 +580,7 @@ void
 dSMatrix::SSOR_precond (dVector& dst, const dVector& src, const double om)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
 
   int p;
   const unsigned int  n = src.size();
@@ -541,6 +615,7 @@ void
 dSMatrix::SOR_precond (dVector& dst, const dVector& src, const double om)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   dst = src;
   SOR(dst,om);
 }
@@ -549,6 +624,7 @@ void
 dSMatrix::SOR (dVector& dst, const double om)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   Assert(n() == m(), ExcDimensionsDontMatch(n(),m()));
   Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
 
@@ -569,6 +645,7 @@ void
 dSMatrix::SSOR (dVector& dst, const double om)
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
 
   int p;
   const unsigned int  n = dst.size();
@@ -616,6 +693,7 @@ const dSMatrixStruct & dSMatrix::get_sparsity_pattern () const {
 
 void dSMatrix::print (ostream &out) const {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
 
   for (unsigned int i=0; i<cols->rows; ++i)
     for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
@@ -626,6 +704,7 @@ void dSMatrix::print (ostream &out) const {
 
 void dSMatrix::print_formatted (ostream &out, const unsigned int precision) const {
   Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
   out.precision (precision);
   out.setf (ios::scientific, ios::floatfield);   // set output format
   
@@ -642,3 +721,4 @@ void dSMatrix::print_formatted (ostream &out, const unsigned int precision) cons
 
   out.setf (0, ios::floatfield);                 // reset output format
 };
+
index 4b67a950ed0581a3f476dd156a726fc389a467d3..8543c19d4f00d1252b70090d2e1c630d20b7e658 100644 (file)
@@ -9,7 +9,7 @@
 #include <algorithm>
 
 
-inline double sqr (const double x) {
+static inline double sqr (const double x) {
   return x*x;
 };
 
@@ -24,17 +24,10 @@ dVector::dVector () :
 
 dVector::dVector (const unsigned int n) :
                dim(n),
-               maxdim(n),
+               maxdim(0),
                val(0)
 {
-  Assert (n>0, ExcInvalidNumber(n));
-
-  if (n)
-    {
-      val = new double[maxdim];
-      Assert (val != 0, ExcOutOfMemory());
-      clear ();
-    }
+  reinit (n, false);
 }
 
 
@@ -54,38 +47,33 @@ dVector::dVector (const dVector& v) :
 
 
 
-void dVector::reinit (const unsigned int n, const bool fast)
-{
-  Assert (n>0, ExcInvalidNumber(n));
-
+void dVector::reinit (const unsigned int n, const bool fast) {
+  if (n==0) 
+    {
+      if (val) delete[] val;
+      val = 0;
+      maxdim = dim = 0;
+      return;
+    };
+  
   if (n>maxdim)
-  {
-    if (val) delete[] val;
-    val = new double[n];
-    Assert (val != 0, ExcOutOfMemory());
-    maxdim = n;
-  }
+    {
+      if (val) delete[] val;
+      val = new double[n];
+      Assert (val != 0, ExcOutOfMemory());
+      maxdim = n;
+    };
   dim = n;
-  if (!fast)
+  if (fast == false)
     clear ();
 }
 
 
 
-void dVector::reinit (const dVector& v, const bool fast)
-{
-  const unsigned int n = v.size();
-  if (n>maxdim)
-  {
-    if (val) delete[] val;
-    val = new double[n];
-    Assert (val != 0, ExcOutOfMemory());
-    maxdim = n;
-  }
-  dim = n;
-  if (!fast)
-    clear ();
-}
+void dVector::reinit (const dVector& v, const bool fast) {
+  reinit (v.size(), fast);
+};
+
 
 
 
@@ -97,12 +85,15 @@ dVector::~dVector ()
 
 
 void dVector::clear () {
-  fill (begin(), end(), 0.);
+  if (dim>0)
+    fill (begin(), end(), 0.);
 }
 
 
 
 bool dVector::all_zero () const {
+  Assert (dim!=0, ExcEmptyVector());
+  
   const_iterator p = begin(),
                 e = end();
   while (p!=e)
@@ -115,6 +106,8 @@ bool dVector::all_zero () const {
 
 double dVector::operator * (const dVector& v) const
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   if (&v == this)
     return norm_sqr();
   
@@ -149,6 +142,8 @@ double dVector::operator * (const dVector& v) const
 
 double dVector::norm_sqr () const
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   double sum0 = 0,
         sum1 = 0,
         sum2 = 0,
@@ -177,6 +172,8 @@ double dVector::norm_sqr () const
 
 double dVector::mean_value () const
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   double sum0 = 0,
         sum1 = 0,
         sum2 = 0,
@@ -205,6 +202,8 @@ double dVector::mean_value () const
 
 double dVector::l1_norm () const
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   double sum0 = 0,
         sum1 = 0,
         sum2 = 0,
@@ -239,6 +238,8 @@ double dVector::l2_norm () const
 
 
 double dVector::linfty_norm () const {
+  Assert (dim!=0, ExcEmptyVector());
+
   double max0=0.,
         max1=0.,
         max2=0.,
@@ -265,6 +266,8 @@ double dVector::linfty_norm () const {
 
 dVector& dVector::operator += (const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   add (v);
   return *this;
 }
@@ -273,7 +276,9 @@ dVector& dVector::operator += (const dVector& v)
 
 dVector& dVector::operator -= (const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
+
   iterator i_ptr = begin(),
           i_end = end();
   const_iterator v_ptr = v.begin();
@@ -287,6 +292,8 @@ dVector& dVector::operator -= (const dVector& v)
 
 void dVector::add (const double v)
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   iterator i_ptr = begin(),
           i_end = end();
   while (i_ptr!=i_end)
@@ -297,7 +304,9 @@ void dVector::add (const double v)
 
 void dVector::add (const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
+
   iterator i_ptr = begin(),
           i_end = end();
   const_iterator v_ptr = v.begin();
@@ -309,7 +318,9 @@ void dVector::add (const dVector& v)
 
 void dVector::add (const double a, const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
+
   iterator i_ptr = begin(),
           i_end = end();
   const_iterator v_ptr = v.begin();
@@ -322,6 +333,7 @@ void dVector::add (const double a, const dVector& v)
 void dVector::add (const double a, const dVector& v,
                   const double b, const dVector& w)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
   iterator i_ptr = begin(),
@@ -336,6 +348,7 @@ void dVector::add (const double a, const dVector& v,
 
 void dVector::sadd (const double x, const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   iterator i_ptr = begin(),
           i_end = end();
@@ -348,6 +361,7 @@ void dVector::sadd (const double x, const dVector& v)
 
 void dVector::sadd (const double x, const double a, const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   iterator i_ptr = begin(),
           i_end = end();
@@ -361,6 +375,7 @@ void dVector::sadd (const double x, const double a, const dVector& v)
 void dVector::sadd (const double x, const double a,
                    const dVector& v, const double b, const dVector& w)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
   iterator i_ptr = begin(),
@@ -377,6 +392,7 @@ void dVector::sadd (const double x, const double a,
                    const dVector& v, const double b,
                    const dVector& w, const double c, const dVector& y)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   Assert (dim == w.dim, ExcDimensionsDontMatch(dim, w.dim));
   Assert (dim == y.dim, ExcDimensionsDontMatch(dim, y.dim));
@@ -394,6 +410,8 @@ void dVector::sadd (const double x, const double a,
 
 void dVector::scale (const double factor)
 {
+  Assert (dim!=0, ExcEmptyVector());
+
   iterator ptr=begin(), eptr=end();
   while (ptr!=eptr)
     *ptr++ *= factor;
@@ -404,6 +422,7 @@ void dVector::scale (const double factor)
 void dVector::equ (const double a, const dVector& u,
                   const double b, const dVector& v)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim));
   Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim));
   iterator i_ptr = begin(),
@@ -418,6 +437,7 @@ void dVector::equ (const double a, const dVector& u,
 
 void dVector::equ (const double a, const dVector& u)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (dim == u.dim, ExcDimensionsDontMatch(dim, u.dim));
   iterator i_ptr = begin(),
           i_end = end();
@@ -429,6 +449,7 @@ void dVector::equ (const double a, const dVector& u)
 
 
 void dVector::ratio (const dVector &a, const dVector &b) {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (a.dim == b.dim, ExcDimensionsDontMatch (a.dim, b.dim));
 
                                   // no need to reinit with zeros, since
@@ -446,6 +467,7 @@ void dVector::ratio (const dVector &a, const dVector &b) {
 
 dVector& dVector::operator = (const double s)
 {
+  Assert (dim!=0, ExcEmptyVector());
   fill (begin(), end(), s);
   return *this;
 }
@@ -457,7 +479,8 @@ dVector& dVector::operator = (const dVector& v)
   if (v.dim != dim)
     reinit (v.dim, true);
 
-  copy (v.begin(), v.end(), begin());
+  if (dim!=0)
+    copy (v.begin(), v.end(), begin());
   return *this;
 }
 
@@ -466,6 +489,7 @@ dVector& dVector::operator = (const dVector& v)
 void dVector::cadd (const unsigned int i, const VectorBase& V,
                    const double s, const unsigned int j)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -480,6 +504,7 @@ void dVector::cadd(unsigned int i, const VectorBase& V,
                   double s, unsigned int j,
                   double t, unsigned int k)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -497,6 +522,7 @@ void dVector::cadd (const unsigned int i, const VectorBase& V,
                    const double q, const unsigned int l,
                    const double r, const unsigned int m)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -512,6 +538,7 @@ void dVector::cadd (const unsigned int i, const VectorBase& V,
 
 void dVector::czero (const unsigned int i)
 {
+  Assert (dim!=0, ExcEmptyVector());
   Assert (i<dim, ExcInvalidIndex(i,dim));
   val[i] = 0.;
 }
@@ -521,6 +548,7 @@ void dVector::czero (const unsigned int i)
 void dVector::cequ (const unsigned int i, const VectorBase& V,
                    const double s, const unsigned int j)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -535,6 +563,7 @@ void dVector::cequ (const unsigned int i, const VectorBase& V,
                    const double s, const unsigned int j,
                    const double t, const unsigned int k)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -552,6 +581,7 @@ void dVector::cequ (const unsigned int i, const VectorBase& V,
                    const double q, const unsigned int l,
                    const double r, const unsigned int m)
 {
+  Assert (dim!=0, ExcEmptyVector());
   const dVector& v = (const dVector&) V;
 
   Assert (i<dim, ExcInvalidIndex(i,dim));
@@ -574,6 +604,7 @@ const char* dVector::name () const
 
 void dVector::print (FILE* f, const char* format) const
 {
+  Assert (dim!=0, ExcEmptyVector());
   if (!format) format = " %5.2f";
   for (unsigned int j=0;j<size();j++)
     fprintf(f, format, val[j]);
@@ -584,6 +615,7 @@ void dVector::print (FILE* f, const char* format) const
 
 void dVector::print (const char* format) const
 {
+  Assert (dim!=0, ExcEmptyVector());
   if (!format) format = " %5.2f";
   for (unsigned int j=0;j<size();j++)
     printf (format, val[j]);
@@ -593,6 +625,11 @@ void dVector::print (const char* format) const
 
 
 void dVector::print (ostream &out) const {
+  Assert (dim!=0, ExcEmptyVector());
   for (unsigned int i=0; i<size(); ++i)
     out << val[i] << endl;
 };
+
+
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.