From 5a8420cbd67f570d5bfc9411aaa70ce5cf8ef202 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 22 Mar 2015 06:59:16 -0500 Subject: [PATCH] Improve the error messages we get from using Jacobi/SOR/SSOR with a matrix that has zeros on the diagonal. --- include/deal.II/lac/sparse_matrix.templates.h | 54 +++++++++++++++---- 1 file changed, 43 insertions(+), 11 deletions(-) diff --git a/include/deal.II/lac/sparse_matrix.templates.h b/include/deal.II/lac/sparse_matrix.templates.h index 2e80506e7e..94a0bebc29 100644 --- a/include/deal.II/lac/sparse_matrix.templates.h +++ b/include/deal.II/lac/sparse_matrix.templates.h @@ -1314,6 +1314,30 @@ SparseMatrix::residual (Vector &dst, } +namespace +{ + template + void check_nonzero_diagonal (const SparseMatrix &matrix) + { + for (typename SparseMatrix::size_type row=0; row template @@ -1328,6 +1352,8 @@ SparseMatrix::precondition_Jacobi (Vector &dst, AssertDimension (dst.size(), n()); AssertDimension (src.size(), n()); + check_nonzero_diagonal(*this); + const size_type n = src.size(); somenumber *dst_ptr = dst.begin(); const somenumber *src_ptr = src.begin(); @@ -1372,6 +1398,8 @@ SparseMatrix::precondition_SSOR (Vector &dst, AssertDimension (dst.size(), n()); AssertDimension (src.size(), n()); + check_nonzero_diagonal(*this); + const size_type n = src.size(); const std::size_t *rowstart_ptr = &cols->rowstart[0]; somenumber *dst_ptr = &dst(0); @@ -1398,9 +1426,8 @@ SparseMatrix::precondition_SSOR (Vector &dst, // divide by diagonal element *dst_ptr -= s * om; - Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero()); *dst_ptr /= val[*rowstart_ptr]; - }; + } rowstart_ptr = &cols->rowstart[0]; dst_ptr = &dst(0); @@ -1420,7 +1447,6 @@ SparseMatrix::precondition_SSOR (Vector &dst, s += val[j] * dst(cols->colnums[j]); *dst_ptr -= s * om; - Assert(val[*rowstart_ptr]!= 0., ExcDivideByZero()); *dst_ptr /= val[*rowstart_ptr]; }; return; @@ -1524,6 +1550,8 @@ SparseMatrix::SOR (Vector &dst, AssertDimension (m(), n()); AssertDimension (dst.size(), n()); + check_nonzero_diagonal(*this); + for (size_type row=0; row::SOR (Vector &dst, s -= val[j] * dst(col); } - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); dst(row) = s * om / val[cols->rowstart[row]]; } } @@ -1551,6 +1578,8 @@ SparseMatrix::TSOR (Vector &dst, AssertDimension (m(), n()); AssertDimension (dst.size(), n()); + check_nonzero_diagonal(*this); + size_type row=m()-1; while (true) { @@ -1559,7 +1588,6 @@ SparseMatrix::TSOR (Vector &dst, if (cols->colnums[j] > row) s -= val[j] * dst(cols->colnums[j]); - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); dst(row) = s * om / val[cols->rowstart[row]]; if (row == 0) @@ -1588,6 +1616,8 @@ SparseMatrix::PSOR (Vector &dst, Assert (m() == inverse_permutation.size(), ExcDimensionMismatch(m(), inverse_permutation.size())); + check_nonzero_diagonal(*this); + for (size_type urow=0; urow::PSOR (Vector &dst, } } - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); dst(row) = s * om / val[cols->rowstart[row]]; } } @@ -1626,6 +1655,8 @@ SparseMatrix::TPSOR (Vector &dst, Assert (m() == inverse_permutation.size(), ExcDimensionMismatch(m(), inverse_permutation.size())); + check_nonzero_diagonal(*this); + for (size_type urow=m(); urow != 0;) { --urow; @@ -1638,7 +1669,6 @@ SparseMatrix::TPSOR (Vector &dst, s -= val[j] * dst(col); } - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); dst(row) = s * om / val[cols->rowstart[row]]; } } @@ -1689,6 +1719,8 @@ SparseMatrix::SOR_step (Vector &v, Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); + check_nonzero_diagonal(*this); + for (size_type row=0; row::SOR_step (Vector &v, { s -= val[j] * v(cols->colnums[j]); } - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); v(row) += s * om / val[cols->rowstart[row]]; } } @@ -1716,6 +1747,8 @@ SparseMatrix::TSOR_step (Vector &v, Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); + check_nonzero_diagonal(*this); + for (int row=m()-1; row>=0; --row) { somenumber s = b(row); @@ -1723,7 +1756,6 @@ SparseMatrix::TSOR_step (Vector &v, { s -= val[j] * v(cols->colnums[j]); } - Assert(val[cols->rowstart[row]]!= 0., ExcDivideByZero()); v(row) += s * om / val[cols->rowstart[row]]; } } @@ -1757,6 +1789,8 @@ SparseMatrix::SSOR (Vector &dst, AssertDimension (m(), n()); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); + check_nonzero_diagonal(*this); + const size_type n = dst.size(); size_type j; somenumber s; @@ -1773,7 +1807,6 @@ SparseMatrix::SSOR (Vector &dst, } } dst(i) -= s * om; - Assert(val[cols->rowstart[i]]!= 0., ExcDivideByZero()); dst(i) /= val[cols->rowstart[i]]; } @@ -1788,7 +1821,6 @@ SparseMatrix::SSOR (Vector &dst, if (static_cast(i)rowstart[i]]!= 0., ExcDivideByZero()); dst(i) -= s * om / val[cols->rowstart[i]]; } } -- 2.39.5