From 23a853a235abc40ba350715e8d07d6f348042aad Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 29 Feb 2008 23:03:09 +0000 Subject: [PATCH] The implementation of the ILU decomposition was rather ineffecient. It has been replaced by the one from the book by Saad. git-svn-id: https://svn.dealii.org/trunk@15821 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 9 + deal.II/lac/include/lac/sparse_ilu.h | 21 +- .../lac/include/lac/sparse_ilu.templates.h | 193 ++++++++---------- 3 files changed, 103 insertions(+), 120 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 0fcbdb44b3..210c49a359 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -214,6 +214,15 @@ an integer id of the current thread.

lac

    +
  1. Fixed: The implementation of SparseILU::decompose was rather +inefficient in that it accessed random elements of the matrix in its +inner loop. It has been replaced by the algorithm given in the book +by Yves Saad: "Iterative methods for sparse linear systems", second +edition, in section 10.3.2. +
    +(WB 2008/2/29) +
    2. +
  2. Fixed: The implementation of SparseILU::vmult very needlessly called SparseMatrix::vmult just to throw away the (nonsensical) result away immediately. This is now fixed. diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index 77b85ecdfa..28b8565290 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -31,21 +31,9 @@ DEAL_II_NAMESPACE_OPEN * in the decomposition that do not fit into the sparsity structure of this * object are discarded. * - * The algorithm used by this class is as follows (indices run from @p 0 - * to @p N-1): - * @verbatim - * copy original matrix into a[i,j] - * - * for i=1..N-1 - * a[i-1,i-1] = a[i-1,i-1]^{-1} - * - * for k=0..i-1 - * a[i,k] = a[i,k] * a[k,k] - * - * for j=k+1..N-1 - * if (a[i,j] exists & a[k,j] exists) - * a[i,j] -= a[i,k] * a[k,j] - * @endverbatim + * The algorithm used by this class is essentially a copy of the + * algorithm given in the book Y. Saad: "Iterative methods for sparse + * linear systems", second edition, in section 10.3.2. * * *

    Usage and state management

    @@ -58,8 +46,7 @@ DEAL_II_NAMESPACE_OPEN * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). * - * @author Wolfgang Bangerth, 1999, based on a similar implementation - * by Malte Braack; unified interface: Ralf Hartmann + * @author Wolfgang Bangerth, 2008; unified interface: Ralf Hartmann */ template class SparseILU : public SparseLUDecomposition diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index fb3d6722ab..c7cd70031b 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -20,6 +20,7 @@ #include #include + DEAL_II_NAMESPACE_OPEN template @@ -66,113 +67,99 @@ void SparseILU::decompose (const SparseMatrix &matrix, if (strengthen_diagonal>0) this->strengthen_diagonal_impl(); - const SparsityPattern &sparsity = this->get_sparsity_pattern(); - const std::size_t * const rowstart_indices = sparsity.get_rowstart_indices(); - const unsigned int * const column_numbers = sparsity.get_column_numbers(); - -/* - PSEUDO-ALGORITHM - (indices=0..N-1) - - for i=1..N-1 - a[i-1,i-1] = a[i-1,i-1]^{-1} - - for k=0..i-1 - a[i,k] = a[i,k] * a[k,k] - - for j=k+1..N-1 - if (a[i,j] exists & a[k,j] exists) - a[i,j] -= a[i,k] * a[k,j] -*/ + // in the following, we implement + // algorithm 10.4 in the book by + // Saad by translating in essence + // the algorithm given at the end + // of section 10.3.2, using the + // names of variables used there + const SparsityPattern &sparsity = this->get_sparsity_pattern(); + const std::size_t * const ia = sparsity.get_rowstart_indices(); + const unsigned int * const ja = sparsity.get_column_numbers(); + number * luval = &this->global_entry (0); - // i := row - for (unsigned int row=1; rowm(); ++row) + const unsigned int N = this->m(); + + std::vector iw (N, numbers::invalid_unsigned_int); + + for (unsigned int k=0; kglobal_entry(rowstart_indices[row-1]) != 0), - ExcDivideByZero()); + const unsigned int j1 = ia[k], + j2 = ia[k+1]-1; + + for (unsigned int j=j1; j<=j2; ++j) + iw[ja[j]] = j; + + // the algorithm in the book + // works on the elements of row + // k left of the + // diagonal. however, since we + // store the diagonal element + // at the first position, start + // at the element after the + // diagonal and run as long as + // we don't walk into the right + // half + unsigned int j = j1+1; + + label_150: - this->global_entry (rowstart_indices[row-1]) - = 1./this->global_entry (rowstart_indices[row-1]); - - // let k run over all lower-left - // elements of row i; skip - // diagonal element at start - const unsigned int * first_of_row - = &column_numbers[rowstart_indices[row]+1]; - const unsigned int * first_after_diagonal = this->prebuilt_lower_bound[row]; - - // k := *col_ptr - for (const unsigned int * col_ptr = first_of_row; - col_ptr!=first_after_diagonal; ++col_ptr) - { - const unsigned int global_index_ik = col_ptr-column_numbers; - this->global_entry(global_index_ik) *= this->diag_element(*col_ptr); - - // now do the inner loop over - // j. note that we need to do - // it in the right order, i.e. - // taking into account that the - // columns are sorted within each - // row correctly, but excluding - // the main diagonal entry - // - // the explicit use of operator() - // works around a bug in some gcc - // versions (see PR 18803) - const int global_index_ki = sparsity.operator()(*col_ptr,row); - - if (global_index_ki != -1) - this->diag_element(row) -= this->global_entry(global_index_ik) * - this->global_entry(global_index_ki); - - for (const unsigned int * j = col_ptr+1; - j<&column_numbers[rowstart_indices[row+1]]; - ++j) - { - // get the locations of - // entries ij and kj in - // the matrix. note that kglobal_entry(global_index_ij) -= this->global_entry(global_index_ik) * - this->global_entry(global_index_kj); - }; - }; - }; - - // Here the very last diagonal - // element still has to be inverted - // because the for-loop doesn't do - // it... - this->diag_element(this->m()-1) = 1./this->diag_element(this->m()-1); + unsigned int jrow = ja[j]; + if (jrow >= k) + goto label_200; + + // actual computations: + { + number t1 = luval[j] * luval[ia[jrow]]; + luval[j] = t1; + + // jj runs from just right of + // the diagonal to the end of + // the row + unsigned int jj = ia[jrow]+1; + while (ja[jj] < jrow) + ++jj; + for (; jj k) || (j==ia[k+1]), ExcInternalError()); + + // now we have to deal with the + // diagonal element. in the + // book it is located at + // position 'j', but here we + // use the convention of + // storing the diagonal element + // first, so instead of j we + // use uptr[k]=ia[k] + Assert (luval[ia[k]] != 0, ExcInternalError()); + + luval[ia[k]] = 1./luval[ia[k]]; + + for (unsigned int j=j1; j<=j2; ++j) + iw[ja[j]] = numbers::invalid_unsigned_int; + } } -- 2.39.5