From 3a81c2b7fa8a969d26ce4ab3bd7e2678f8eecef0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 15 Mar 2006 21:21:08 +0000 Subject: [PATCH] Check values written into matrices and some vectors to make sure they're finite git-svn-id: https://svn.dealii.org/trunk@12606 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.html | 10 +++ deal.II/lac/include/lac/block_matrix_base.h | 27 ++++--- .../lac/include/lac/block_sparse_matrix_ez.h | 32 ++++++--- deal.II/lac/include/lac/block_vector.h | 8 +++ deal.II/lac/include/lac/block_vector_base.h | 65 +++++++++++++++++ deal.II/lac/include/lac/full_matrix.h | 2 + .../lac/include/lac/full_matrix.templates.h | 13 +++- deal.II/lac/include/lac/petsc_matrix_base.h | 8 ++- deal.II/lac/include/lac/sparse_matrix.h | 22 ++++-- deal.II/lac/include/lac/sparse_matrix_ez.h | 14 ++++ deal.II/lac/include/lac/vector.h | 11 +++ deal.II/lac/include/lac/vector.templates.h | 45 ++++++++++++ deal.II/lac/source/petsc_matrix_base.cc | 10 ++- deal.II/lac/source/petsc_vector_base.cc | 71 ++++++++++++++++++- 14 files changed, 308 insertions(+), 30 deletions(-) diff --git a/deal.II/doc/news/changes.html b/deal.II/doc/news/changes.html index 5d6ed58050..20e0db42fc 100644 --- a/deal.II/doc/news/changes.html +++ b/deal.II/doc/news/changes.html @@ -409,6 +409,16 @@ inconvenience this causes.

lac

    +
  1. + 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. +
    + (Stephan Kramer, WB 2006/03/15) +

    +
  2. Changed: There are new FullMatrix::(i,j) * to value. Throws an * error if the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. + * exist or if value 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 value to the element - * (i,j). Throws an error if - * the entry does not - * exist. Still, it is allowed to - * store zero values in + * Add value to the + * element (i,j). + * Throws an error if the entry + * does not exist or if + * value 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::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 row_index = row_block_indices.global_to_local (i), col_index = column_block_indices.global_to_local (j); @@ -1586,6 +1593,10 @@ BlockMatrixBase::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 diff --git a/deal.II/lac/include/lac/block_sparse_matrix_ez.h b/deal.II/lac/include/lac/block_sparse_matrix_ez.h index 12517a8141..4c5ed096d9 100644 --- a/deal.II/lac/include/lac/block_sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/block_sparse_matrix_ez.h @@ -214,12 +214,13 @@ class BlockSparseMatrixEZ : public Subscriptor unsigned int n () const; /** - * Set the element (i,j) 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 (i,j) + * to @p value. Throws an error + * if the entry does not exist or + * if value 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 - * (i,j). Throws an error if - * the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. + * (i,j). Throws an + * error if the entry does not + * exist or if value 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::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 row_index = row_indices.global_to_local (i), col_index = column_indices.global_to_local (j); @@ -432,6 +438,10 @@ BlockSparseMatrixEZ::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 row_index = row_indices.global_to_local (i), col_index = column_indices.global_to_local (j); diff --git a/deal.II/lac/include/lac/block_vector.h b/deal.II/lac/include/lac/block_vector.h index 1417f16abd..495e127337 100644 --- a/deal.II/lac/include/lac/block_vector.h +++ b/deal.II/lac/include/lac/block_vector.h @@ -422,6 +422,10 @@ inline BlockVector & BlockVector::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::operator = (const BlockVector &v) template void BlockVector::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; in_blocks();++i) this->components[i].scale(factor); } diff --git a/deal.II/lac/include/lac/block_vector_base.h b/deal.II/lac/include/lac/block_vector_base.h index 07d4682b89..eb97e4bad4 100644 --- a/deal.II/lac/include/lac/block_vector_base.h +++ b/deal.II/lac/include/lac/block_vector_base.h @@ -1842,6 +1842,10 @@ BlockVectorBase::operator -= (const BlockVectorBase& v) template void BlockVectorBase::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 void BlockVectorBase::add (const value_type a, const BlockVectorBase& 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::add (const value_type a, const value_type b, const BlockVectorBase& 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 void BlockVectorBase::sadd (const value_type x, const BlockVectorBase& 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 void BlockVectorBase::sadd (const value_type x, const value_type a, const BlockVectorBase& 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::sadd (const value_type x, const value_type a, const value_type b, const BlockVectorBase& 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::sadd (const value_type x, const value_type a, const value_type c, const BlockVectorBase& 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::equ (const value_type a, const value_type b, const BlockVectorBase& 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 void BlockVectorBase::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 BlockVectorBase& BlockVectorBase::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 & BlockVectorBase::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 & BlockVectorBase::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::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); } diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index ec667beb4f..ae08654650 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -97,6 +97,10 @@ template FullMatrix & FullMatrix::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 FullMatrix & FullMatrix::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; diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index d6ddcd9b4f..8936ea256b 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -338,7 +338,10 @@ namespace PETScWrappers * contrast to the SparseMatrix * class which throws an error * if the entry does not exist. - */ + * If value 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 value is not a + * finite number an exception + * is thrown. */ void add (const unsigned int i, const unsigned int j, diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index eb75aa6c62..4c37510981 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -686,9 +686,10 @@ class SparseMatrix : public virtual Subscriptor * Set the element (i,j) * to value. Throws an * error if the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. + * exist or if value 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 value to the * element (i,j). Throws * an error if the entry does not - * exist. Still, it is allowed to - * store zero values in - * non-existent fields. + * exist or if value 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::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::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); diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index 2af0484815..b112250626 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -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 value 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 * (i,j). Allocates the entry * if it does not exist. Filters * out zeroes automatically. + * If value 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::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::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 inline Vector & Vector::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 inline Vector & Vector::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 inline Vector & Vector::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; } diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index db9fbf7c78..dc188728ae 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -43,6 +43,8 @@ namespace internal template 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::add (const Vector& v) template void Vector::add (const Number a, const Vector& 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 void Vector::add (const Number a, const Vector& v, const Number b, const Vector& 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::add (const Number a, const Vector& v, template void Vector::sadd (const Number x, const Vector& 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 void Vector::sadd (const Number x, const Number a, const Vector& 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::sadd (const Number x, const Number a, const Vector& v, const Number b, const Vector& 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::sadd (const Number x, const Number a, const Vector& w, const Number c, const Vector& 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::sadd (const Number x, const Number a, template void Vector::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 template void Vector::equ (const Number a, const Vector& 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 void Vector::equ (const Number a, const Vector& u, const Number b, const Vector& 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)); diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index e7774dc451..2a0af2d2a3 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -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; diff --git a/deal.II/lac/source/petsc_vector_base.cc b/deal.II/lac/source/petsc_vector_base.cc index 252a96cbf1..cc91844fb7 100644 --- a/deal.II/lac/source/petsc_vector_base.cc +++ b/deal.II/lac/source/petsc_vector_base.cc @@ -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())); -- 2.39.5