]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Check values written into matrices and some vectors to make sure they're finite
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 15 Mar 2006 21:21:08 +0000 (21:21 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 15 Mar 2006 21:21:08 +0000 (21:21 +0000)
git-svn-id: https://svn.dealii.org/trunk@12606 0785d39b-7218-0410-832d-ea1e28bc413d

14 files changed:
deal.II/doc/news/changes.html
deal.II/lac/include/lac/block_matrix_base.h
deal.II/lac/include/lac/block_sparse_matrix_ez.h
deal.II/lac/include/lac/block_vector.h
deal.II/lac/include/lac/block_vector_base.h
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix_ez.h
deal.II/lac/include/lac/vector.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/petsc_matrix_base.cc
deal.II/lac/source/petsc_vector_base.cc

index 5d6ed580506bd9c7013b454d194010059dc0f91a..20e0db42fc382d39808fb4122901e4f04fdd9454 100644 (file)
@@ -409,6 +409,16 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       Improved: All matrix (and some vector) classes now check whether
+       entries written into them represent finite floating point
+       values. This should catch some bugs earlier where one writes
+       infinite or NaN values into a matrix only to realize later that
+       the linear solver fails.
+       <br>
+       (Stephan Kramer, WB 2006/03/15)
+       </p>  
+
   <li> <p>
        Changed: There are new <code
        class="class">FullMatrix</code>::<code
index a99a311a8f4d4684f4e4b8d30d81befb8071d1b6..7a6348ea65a797c00778619a459d7a8de80bdd5e 100644 (file)
@@ -465,20 +465,23 @@ class BlockMatrixBase : public Subscriptor
                                      * Set the element <tt>(i,j)</tt>
                                      * to <tt>value</tt>. Throws an
                                      * error if the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
-                                     * non-existent fields.
+                                     * exist or if <tt>value</tt> is
+                                     * not a finite number. Still, it
+                                     * is allowed to store zero
+                                     * values in non-existent fields.
                                      */
     void set (const unsigned int i,
              const unsigned int j,
              const value_type value);
     
                                     /**
-                                     * Add <tt>value</tt> to the element
-                                     * <tt>(i,j)</tt>.  Throws an error if
-                                     * the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
+                                     * Add <tt>value</tt> to the
+                                     * element <tt>(i,j)</tt>.
+                                     * Throws an error if the entry
+                                     * does not exist or if
+                                     * <tt>value</tt> is not a finite
+                                     * number. Still, it is allowed
+                                     * to store zero values in
                                      * non-existent fields.
                                      */
     void add (const unsigned int i,
@@ -1569,6 +1572,10 @@ BlockMatrixBase<MatrixType>::set (const unsigned int i,
                                   const unsigned int j,
                                   const value_type value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   const std::pair<unsigned int,unsigned int>
     row_index = row_block_indices.global_to_local (i),
     col_index = column_block_indices.global_to_local (j);
@@ -1586,6 +1593,10 @@ BlockMatrixBase<MatrixType>::add (const unsigned int i,
                                   const unsigned int j,
                                   const value_type value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
                                    // save some cycles for zero additions, but
                                    // only if it is safe for the matrix we are
                                    // working with
index 12517a8141c537f9111de4c6806cd1652437f8d4..4c5ed096d92761346afc8faca0a325202e6b2c82 100644 (file)
@@ -214,12 +214,13 @@ class BlockSparseMatrixEZ : public Subscriptor
     unsigned int n () const;
 
                                     /**
-                                     * Set the element <tt>(i,j)</tt> to
-                                     * @p value.  Throws an error if
-                                     * the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
-                                     * non-existent fields.
+                                     * Set the element <tt>(i,j)</tt>
+                                     * to @p value.  Throws an error
+                                     * if the entry does not exist or
+                                     * if <tt>value</tt> is not a
+                                     * finite number. Still, it is
+                                     * allowed to store zero values
+                                     * in non-existent fields.
                                      */
     void set (const unsigned int i,
              const unsigned int j,
@@ -227,11 +228,12 @@ class BlockSparseMatrixEZ : public Subscriptor
     
                                     /**
                                      * Add @p value to the element
-                                     * <tt>(i,j)</tt>.  Throws an error if
-                                     * the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
-                                     * non-existent fields.
+                                     * <tt>(i,j)</tt>.  Throws an
+                                     * error if the entry does not
+                                     * exist or if <tt>value</tt> is
+                                     * not a finite number. Still, it
+                                     * is allowed to store zero
+                                     * values in non-existent fields.
                                      */
     void add (const unsigned int i, const unsigned int j,
              const Number value);
@@ -415,6 +417,10 @@ BlockSparseMatrixEZ<Number>::set (const unsigned int i,
                                  const unsigned int j,
                                  const Number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   const std::pair<unsigned int,unsigned int>
     row_index = row_indices.global_to_local (i),
     col_index = column_indices.global_to_local (j);
@@ -432,6 +438,10 @@ BlockSparseMatrixEZ<Number>::add (const unsigned int i,
                                  const unsigned int j,
                                  const Number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   const std::pair<unsigned int,unsigned int>
     row_index = row_indices.global_to_local (i),
     col_index = column_indices.global_to_local (j);
index 1417f16abd05e195e0f714672dbaf0a60d7f2a45..495e1273378fbae8811f9afd679709f7ae70fbc1 100644 (file)
@@ -422,6 +422,10 @@ inline
 BlockVector<Number> &
 BlockVector<Number>::operator = (const value_type s)
 {
+
+  Assert (deal_II_numbers::is_finite(s), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
   BaseClass::operator = (s);
   return *this;
 }
@@ -453,6 +457,10 @@ BlockVector<Number>::operator = (const BlockVector<Number2> &v)
 template <typename Number>
 void BlockVector<Number>::scale (const value_type factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   for (unsigned int i=0; i<this->n_blocks();++i)
     this->components[i].scale(factor);
 }
index 07d4682b8923a2a850c10a97d8cb60967af9069a..eb97e4bad41427d5a13b1caa5cd79c37f2cf616c 100644 (file)
@@ -1842,6 +1842,10 @@ BlockVectorBase<VectorType>::operator -= (const BlockVectorBase<VectorType>& v)
 template <class VectorType>
 void BlockVectorBase<VectorType>::add (const value_type a)
 {
+
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   for (unsigned int i=0;i<n_blocks();++i)
     {
       components[i].add(a);
@@ -1868,6 +1872,10 @@ template <class VectorType>
 void BlockVectorBase<VectorType>::add (const value_type a,
                                        const BlockVectorBase<VectorType>& v)
 {
+
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   
@@ -1885,6 +1893,12 @@ void BlockVectorBase<VectorType>::add (const value_type a,
                                        const value_type b,
                                        const BlockVectorBase<VectorType>& w)
 {
+
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   Assert (n_blocks() == w.n_blocks(),
@@ -1903,6 +1917,10 @@ template <class VectorType>
 void BlockVectorBase<VectorType>::sadd (const value_type x,
                                         const BlockVectorBase<VectorType>& v)
 {
+
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   
@@ -1918,6 +1936,12 @@ template <class VectorType>
 void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                         const BlockVectorBase<VectorType>& v)
 {
+
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   
@@ -1935,6 +1959,14 @@ void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                         const value_type b,
                                         const BlockVectorBase<VectorType>& w)
 {
+
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   Assert (n_blocks() == w.n_blocks(),
@@ -1956,6 +1988,16 @@ void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                         const value_type c,
                                         const BlockVectorBase<VectorType>& y)
 {
+
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(c), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));  
   Assert (n_blocks() == w.n_blocks(),
@@ -1990,6 +2032,12 @@ void BlockVectorBase<VectorType>::equ (const value_type a,
                                        const value_type b,
                                        const BlockVectorBase<VectorType>& w)
 {
+
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   Assert (n_blocks() == w.n_blocks(),
@@ -2007,6 +2055,10 @@ template <class BlockVector2>
 void BlockVectorBase<VectorType>::equ (const value_type    a,
                                        const BlockVector2 &v)
 {
+
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
   
@@ -2019,6 +2071,10 @@ template <class VectorType>
 BlockVectorBase<VectorType>&
 BlockVectorBase<VectorType>::operator = (const value_type s)
 {
+
+  Assert (deal_II_numbers::is_finite(s), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   for (unsigned int i=0;i<n_blocks();++i)
     components[i] = s;
     
@@ -2081,6 +2137,10 @@ inline
 BlockVectorBase<VectorType> &
 BlockVectorBase<VectorType>::operator *= (const value_type factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   for (unsigned int i=0;i<n_blocks();++i)
     components[i] *= factor;
 
@@ -2094,6 +2154,11 @@ inline
 BlockVectorBase<VectorType> &
 BlockVectorBase<VectorType>::operator /= (const value_type factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (factor > 0., ExcZero() );
+
   for (unsigned int i=0;i<n_blocks();++i)
     components[i] /= factor;
 
index 8e4b1b6b5d87fbf148586deb3d5714c67a62385a..9d14820f82de347aae0f3a4a9ec984bf6f0832ae 100644 (file)
@@ -1258,6 +1258,8 @@ inline
 number
 FullMatrix<number>::Accessor::value() const
 {
+  Assert (deal_II_numbers::is_finite( matrix->el(a_row, a_col) ), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
   return matrix->el(a_row, a_col);
 }
 
index ec667beb4fa2d65145117c9a2be2f412cd04f675..ae08654650121444e95a86058d9b28d1eb50cccd 100644 (file)
@@ -97,6 +97,10 @@ template <typename number>
 FullMatrix<number> &
 FullMatrix<number>::operator *= (const double factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   number       *p = &this->el(0,0);
   const number *e = &this->el(0,0) + n()*m();
   while (p != e)
@@ -111,11 +115,18 @@ template <typename number>
 FullMatrix<number> &
 FullMatrix<number>::operator /= (const double factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   number       *p = &this->el(0,0);
   const number *e = &this->el(0,0) + n()*m();
 
   const number factor_inv = 1./factor;
-  
+
+  Assert (deal_II_numbers::is_finite(factor_inv), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
   while (p != e)
     *p++ *= factor_inv;
 
index d6ddcd9b4fee574fceb8a0126c947acf8e0f1564..8936ea256bb198b0135ab1784803497f0d0ff31e 100644 (file)
@@ -338,7 +338,10 @@ namespace PETScWrappers
                                        * contrast to the SparseMatrix
                                        * class which throws an error
                                        * if the entry does not exist.
-                                        */
+                                       * If <tt>value</tt> is not a
+                                       * finite number an exception
+                                       * is thrown.
+                                       */
       void set (const unsigned int i,
                 const unsigned int j,
                 const PetscScalar value);
@@ -357,6 +360,9 @@ namespace PETScWrappers
                                        * contrast to the SparseMatrix
                                        * class which throws an error
                                        * if the entry does not exist.
+                                       * If <tt>value</tt> is not a
+                                       * finite number an exception
+                                       * is thrown.
                                         */
       void add (const unsigned int i,
                 const unsigned int j,
index eb75aa6c62bae337a0215f11e1eeb692a32eedea..4c375109817c80e1574284f67a799e87bbf050d3 100644 (file)
@@ -686,9 +686,10 @@ class SparseMatrix : public virtual Subscriptor
                                      * Set the element (<i>i,j</i>)
                                      * to <tt>value</tt>. Throws an
                                      * error if the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
-                                     * non-existent fields.
+                                     * exist or if <tt>value</tt> is
+                                     * not a finite number. Still, it
+                                     * is allowed to store zero
+                                     * values in non-existent fields.
                                      */
     void set (const unsigned int i,
               const unsigned int j,
@@ -710,9 +711,10 @@ class SparseMatrix : public virtual Subscriptor
                                      * Add <tt>value</tt> to the
                                      * element (<i>i,j</i>).  Throws
                                      * an error if the entry does not
-                                     * exist. Still, it is allowed to
-                                     * store zero values in
-                                     * non-existent fields.
+                                     * exist or if <tt>value</tt> is
+                                     * not a finite number. Still, it
+                                     * is allowed to store zero
+                                     * values in non-existent fields.
                                      */
     void add (const unsigned int i,
               const unsigned int j,
@@ -1746,6 +1748,10 @@ void SparseMatrix<number>::set (const unsigned int i,
                                const unsigned int j,
                                const number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (cols != 0, ExcNotInitialized());
                                   // it is allowed to set elements of
                                   // the matrix that are not part of
@@ -1768,6 +1774,10 @@ void SparseMatrix<number>::add (const unsigned int i,
                                const unsigned int j,
                                const number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (cols != 0, ExcNotInitialized());
 
   const unsigned int index = cols->operator()(i,j);
index 2af048481502d5c1251200f844100b0d7aab129f..b112250626ad2d3c0af80939c0e09b65d5b12dc4 100644 (file)
@@ -461,6 +461,9 @@ class SparseMatrixEZ : public Subscriptor
                                      * @p value. Allocates the entry,
                                      * if it does not exist and
                                      * @p value is non-zero.
+                                     * If <tt>value</tt> is not a
+                                     * finite number an exception
+                                     * is thrown.
                                      */
     void set (const unsigned int i, const unsigned int j,
              const number value);
@@ -470,6 +473,9 @@ class SparseMatrixEZ : public Subscriptor
                                      * <tt>(i,j)</tt>. Allocates the entry
                                      * if it does not exist. Filters
                                      * out zeroes automatically.
+                                     * If <tt>value</tt> is not a
+                                     * finite number an exception
+                                     * is thrown.
                                      */
     void add (const unsigned int i, const unsigned int j,
              const number value);
@@ -1474,6 +1480,10 @@ void SparseMatrixEZ<number>::set (const unsigned int i,
                                  const unsigned int j,
                                  const number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (i<m(), ExcIndexRange(i,0,m()));
   Assert (j<n(), ExcIndexRange(j,0,n()));
 
@@ -1500,6 +1510,10 @@ void SparseMatrixEZ<number>::add (const unsigned int i,
                                  const unsigned int j,
                                  const number value)
 {
+
+  Assert (deal_II_numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (i<m(), ExcIndexRange(i,0,m()));
   Assert (j<n(), ExcIndexRange(j,0,n()));
 
index df78d9f203a6713b6a3e3f158914b8f0cfd16730..4f2fcf79127b06add3b3115110811e5e57aebf0c 100644 (file)
@@ -842,6 +842,9 @@ template <typename Number>
 inline
 Vector<Number> & Vector<Number>::operator = (const Number s)
 {
+  Assert (deal_II_numbers::is_finite(s), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  
   if (s != 0.)
     Assert (vec_size!=0, ExcEmptyObject());
   if (vec_size!=0)
@@ -924,6 +927,10 @@ template <typename Number>
 inline
 Vector<Number> & Vector<Number>::operator *= (const Number factor)
 {
+
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   scale (factor);
   return *this;
 }
@@ -934,6 +941,10 @@ template <typename Number>
 inline
 Vector<Number> & Vector<Number>::operator /= (const Number factor)
 {
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (factor > 0., ExcZero() );
+
   *this *= (1./factor);
   return *this;
 }
index db9fbf7c78671a67ad9dc40823042d23ad8f0229..dc188728aee2ebf109c2842bd6f70ec9e315b1c7 100644 (file)
@@ -43,6 +43,8 @@ namespace internal
     template <typename Number>
     inline Number sqr (const Number x)
     {
+      Assert (deal_II_numbers::is_finite(x), 
+             ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
       return x*x;
     }
   }
@@ -437,6 +439,9 @@ void Vector<Number>::add (const Vector<Number>& v)
 template <typename Number>
 void Vector<Number>::add (const Number a, const Vector<Number>& v)
 {
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
 
@@ -452,6 +457,11 @@ template <typename Number>
 void Vector<Number>::add (const Number a, const Vector<Number>& v,
                          const Number b, const Vector<Number>& w)
 {
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
@@ -467,6 +477,9 @@ void Vector<Number>::add (const Number a, const Vector<Number>& v,
 template <typename Number>
 void Vector<Number>::sadd (const Number x, const Vector<Number>& v)
 {
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   iterator i_ptr = begin(),
@@ -481,6 +494,11 @@ template <typename Number>
 void Vector<Number>::sadd (const Number x, const Number a,
                            const Vector<Number>& v)
 {
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   iterator i_ptr = begin(),
@@ -496,6 +514,13 @@ void Vector<Number>::sadd (const Number x, const Number a,
                           const Vector<Number>& v, const Number b,
                            const Vector<Number>& w)
 {
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
@@ -514,6 +539,15 @@ void Vector<Number>::sadd (const Number x, const Number a,
                           const Vector<Number>& w, const Number c,
                            const Vector<Number>& y)
 {
+  Assert (deal_II_numbers::is_finite(x), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(c), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
   Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size));
@@ -533,6 +567,9 @@ void Vector<Number>::sadd (const Number x, const Number a,
 template <typename Number>
 void Vector<Number>::scale (const Number factor)
 {
+  Assert (deal_II_numbers::is_finite(factor), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
 
   iterator             ptr  = begin();
@@ -562,6 +599,9 @@ template <typename Number>
 template <typename Number2>
 void Vector<Number>::equ (const Number a, const Vector<Number2>& u)
 {
+  Assert (deal_II_numbers::is_finite(a), 
+         ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
   iterator i_ptr = begin(),
@@ -576,6 +616,11 @@ template <typename Number>
 void Vector<Number>::equ (const Number a, const Vector<Number>& u,
                          const Number b, const Vector<Number>& v)
 {
+  Assert (deal_II_numbers::is_finite(a), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+  Assert (deal_II_numbers::is_finite(b), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
   Assert (vec_size!=0, ExcEmptyObject());
   Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size));
   Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
index e7774dc451b2bfc25c009590209b8c1713ce68c5..2a0af2d2a34952fd1ec9a91a2bcf557f4cf76083 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -136,6 +136,10 @@ namespace PETScWrappers
                    const unsigned int j,
                    const PetscScalar value)
   {
+
+    Assert (!std::isnan(value) && !std::isinf(value), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
     if (last_action != LastAction::insert)
       {
         int ierr;
@@ -164,6 +168,10 @@ namespace PETScWrappers
                    const unsigned int j,
                    const PetscScalar value)
   {
+
+    Assert (deal_II_numbers::is_finite(value), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    
     if (last_action != LastAction::add)
       {
         int ierr;
index 252a96cbf1d1ebeb60226fa7668ff72da7cf33c3..cc91844fb7aed1d355fd3fe73ad76ddbb2fafcf5 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004, 2005 by the deal.II authors
+//    Copyright (C) 2004, 2005, 2006 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -121,6 +121,10 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator = (const PetscScalar s)
   {
+
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
                                      // flush previously cached elements. this
                                      // seems to be necessary since petsc
                                      // 2.2.1, at least for parallel vectors
@@ -490,6 +494,10 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator *= (const PetscScalar a)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecScale (&a, vector);
 #else
@@ -506,8 +514,15 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator /= (const PetscScalar a)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
     const PetscScalar factor = 1./a;
     
+    Assert (deal_II_numbers::is_finite(factor), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecScale (&factor, vector);
 #else
@@ -558,7 +573,11 @@ namespace PETScWrappers
   void
   VectorBase::add (const PetscScalar s)
   {
- #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
+
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
+#if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecShift (&s, vector);
 #else
     const int ierr = VecShift (vector, s);
@@ -581,6 +600,10 @@ namespace PETScWrappers
   VectorBase::add (const PetscScalar a,
                    const VectorBase     &v)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecAXPY (&a, v, vector);
 #else
@@ -598,6 +621,12 @@ namespace PETScWrappers
                    const PetscScalar b,
                    const VectorBase &w)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(b), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
     const PetscScalar weights[2] = {a,b};
     Vec               addends[2] = {v.vector, w.vector};
     
@@ -616,6 +645,10 @@ namespace PETScWrappers
   VectorBase::sadd (const PetscScalar s,
                     const VectorBase &v)
   {
+
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
 #if (PETSC_VERSION_MAJOR <= 2) && (PETSC_VERSION_MINOR < 3)
     const int ierr = VecAYPX (&s, v, vector);
 #else
@@ -632,6 +665,12 @@ namespace PETScWrappers
                     const PetscScalar a,
                     const VectorBase     &v)
   {
+
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
                                      // there is nothing like a AXPAY
                                      // operation in Petsc, so do it in two
                                      // steps
@@ -648,6 +687,14 @@ namespace PETScWrappers
                     const PetscScalar b,
                     const VectorBase     &w)
   {
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(b), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
                                      // there is no operation like MAXPAY, so
                                      // do it in two steps
     *this *= s;
@@ -675,6 +722,16 @@ namespace PETScWrappers
                     const PetscScalar c,
                     const VectorBase     &x)
   {
+    Assert (deal_II_numbers::is_finite(s), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(b), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(c), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
                                      // there is no operation like MAXPAY, so
                                      // do it in two steps
     *this *= s;
@@ -707,6 +764,10 @@ namespace PETScWrappers
   VectorBase::equ (const PetscScalar a,
                    const VectorBase &v)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
     Assert (size() == v.size(),
             ExcDimensionMismatch (size(), v.size()));
 
@@ -727,6 +788,12 @@ namespace PETScWrappers
                    const PetscScalar b,
                    const VectorBase &w)
   {
+
+    Assert (deal_II_numbers::is_finite(a), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+    Assert (deal_II_numbers::is_finite(b), 
+           ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
     Assert (size() == v.size(),
             ExcDimensionMismatch (size(), v.size()));
 

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.