From 97c2b4a348ac85dc72800460c6163f2f40189127 Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 29 Apr 1999 19:00:24 +0000 Subject: [PATCH] Add sparse ILU decompoition. git-svn-id: https://svn.dealii.org/trunk@1227 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_ilu.h | 202 +++++++++++ .../lac/include/lac/sparse_ilu.templates.h | 339 ++++++++++++++++++ deal.II/lac/source/sparse_ilu.cc | 27 ++ 3 files changed, 568 insertions(+) create mode 100644 deal.II/lac/include/lac/sparse_ilu.h create mode 100644 deal.II/lac/include/lac/sparse_ilu.templates.h create mode 100644 deal.II/lac/source/sparse_ilu.cc diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h new file mode 100644 index 0000000000..6a33d842e2 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -0,0 +1,202 @@ +/*---------------------------- matrix_decompositions.h ---------------------------*/ +/* $Id$ */ +#ifndef __sparse_ilu_H +#define __sparse_ilu_H +/*---------------------------- sparse_ilu.h ---------------------------*/ + + +#include + + + + +/** + * Incomplete LU decomposition of a sparse matrix into another sparse matrix. + * A given matrix is decomposed into a incomplete LU factorization, where + * by incomplete we mean that also a sparse decomposition is used and entries + * 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 #0# + * to #N-1#): + * \begin{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] + * \end{verbatim} + * Using this algorithm, we store the decomposition as a sparse matrix, for + * which the user has to give a sparsity pattern and which is why this + * class is derived from the #SparseMatrix#. Since it is not a matrix in + * the usual sense, the derivation is #protected# rather than #public#. + * + * Note that in the algorithm given, the lower left part of the matrix base + * class is used to store the #L# part of the decomposition, while + * the upper right part is used to store #U#. The diagonal is used to + * store the inverses of the diagonal elements of the decomposition; the + * latter makes the application of the decomposition faster, since inversion + * by the diagonal element has to be done only once, rather than at each + * application (multiplication is much faster than division). + * + * + * \subsection{Fill-in} + * + * The sparse ILU is frequently used with additional fill-in, i.e. the + * sparsity structure of the decomposition is denser than that of the matrix + * to be decomposed. The #decompose# function of this class allows this fill-in + * as long as all entries present in the original matrix are present in the + * decomposition also, i.e. the sparsity pattern of the decomposition is a + * superset of the sparsity pattern in the original matrix. + * + * Such fill-in can be accomplished by various ways, one of which is a + * copy-constructor of the #SparseMatrixStruct# class which allows the addition + * of side-diagonals to a given sparsity structure. + * + * + * \subsection{Use as a preconditioner} + * + * If you want to use an object of this class as a preconditioner for another + * matrix, you can do so by calling the solver function using the following + * sequence, for example (#ilu_sparsity# is some sparsity pattern to be used + * for the decomposition, which you have to create beforehand): + * \begin{verbatim} + * SparseILU ilu (ilu_sparsity); + * ilu.decompose (global_matrix); + * global_matrix.set_preconditioner (ilu); + * + * somesolver.solve (A, x, f, + * PreconditionUseMatrix,Vector > + * (ilu,&SparseILU::template apply_decomposition)); + * \end{verbatim} + * + * + * @author Wolfgang Bangerth, 1999, based on a similar implementation by Malte Braack + */ +template +class SparseILU : protected SparseMatrix +{ + public: + /** + * Constructor; initializes the decomposition + * to be empty, without any structure, i.e. + * it is not usable at all. This + * constructor is therefore only useful + * for objects which are members of a + * class. All other matrices should be + * created at a point in the data flow + * where all necessary information is + * available. + * + * You have to initialize + * the matrix before usage with + * #reinit(SparseMatrixStruct)#. + */ + SparseILU (); + + /** + * Constructor. Takes the given matrix + * sparsity structure to represent the + * sparsity pattern of this decomposition. + * You can change the sparsity pattern later + * on by calling the #reinit# function. + * + * You have to make sure that the lifetime + * of the sparsity structure is at least + * as long as that of this object or as + * long as #reinit# is not called with a + * new sparsity structure. + */ + SparseILU (const SparseMatrixStruct &sparsity); + + /** + * 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. + * + * This function does nothing more than + * passing down to the sparse matrix + * object the call for the same function, + * which is necessary however, since that + * function is not publically visible + * any more. + */ + void reinit (); + + /** + * Reinitialize the sparse matrix with the + * given sparsity pattern. The latter tells + * the matrix how many nonzero elements + * there need to be reserved. + * + * This function does nothing more than + * passing down to the sparse matrix + * object the call for the same function, + * which is necessary however, since that + * function is not publically visible + * any more. + */ + void reinit (const SparseMatrixStruct &sparsity); + + /** + * Perform the incomplete LU factorization + * of the given matrix. + * + * Note that the sparsity structures of + * the decomposition and the matrix passed + * to this function need not be equal, + * but that the pattern used by this + * matrix needs to contain all elements + * used by the matrix to be decomposed. + * Fill-in is thus allowed. + */ + template + void decompose (const SparseMatrix &matrix, + const double strengthen_diagonal=0.); + + /** + * Apply the incomplete decomposition, + * i.e. do one forward-backward step + * $dst=(LU)^{-1}src$. + */ + template + void apply_decomposition (Vector &dst, + const Vector &src) const; + + /** + * Exception + */ + DeclException0 (ExcMatrixNotSquare); + /** + * Exception + */ + DeclException2 (ExcSizeMismatch, + int, int, + << "The sizes " << arg1 << " and " << arg2 + << " of the given objects do not match."); + /** + * Exception + */ + DeclException1 (ExcInvalidStrengthening, + double, + << "The strengthening parameter " << arg1 + << " is not greater or equal than zero!"); +}; + + + + + + +/*---------------------------- sparse_ilu.h ---------------------------*/ +/* end of #ifndef __sparse_ilu_H */ +#endif +/*---------------------------- sparse_ilu.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h new file mode 100644 index 0000000000..3b1de729f5 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -0,0 +1,339 @@ +/*---------------------------- matrix_decompositions.templates.h ---------------------------*/ +/* $Id$ */ +#ifndef __matrix_decompositions_templates_H +#define __matrix_decompositions_templates_H +/*---------------------------- matrix_decompositions.templates.h ---------------------------*/ + + +#include +#include + +#include +#include + + + +template +SparseILU::SparseILU () +{}; + + + +template +SparseILU::SparseILU (const SparseMatrixStruct &sparsity) : + SparseMatrix (sparsity) +{}; + + + +template +void SparseILU::reinit () +{ + SparseMatrix::reinit (); +}; + + + +template +void SparseILU::reinit (const SparseMatrixStruct &sparsity) +{ + SparseMatrix::reinit (sparsity); +}; + + + + +template +template +void SparseILU::decompose (const SparseMatrix &matrix, + const double strengthen_diagonal) +{ + Assert (matrix.m()==matrix.n(), ExcMatrixNotSquare ()); + Assert (m()==n(), ExcMatrixNotSquare ()); + Assert (matrix.m()==m(), ExcSizeMismatch(matrix.m(), m())); + + Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal)); + + + // first thing: copy over all elements + // of #matrix# to the present object + // + // note that some elements in this + // matrix may not be in #matrix#, + // so we need to preset our matrix + // by zeroes. + if (true) + { + // preset the elements + fill_n (&global_entry(0), + n_nonzero_elements(), + 0); + + // note: pointers to the sparsity + // pattern of the old matrix! + const unsigned int * const rowstart_indices + = matrix.get_sparsity_pattern().get_rowstart_indices(); + const int * const column_numbers + = matrix.get_sparsity_pattern().get_column_numbers(); + + for (unsigned int row=0; row 0) + for (unsigned int row=0; row(row)); + + // k := *col_ptr + for (const int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr) + { + const unsigned int global_index_ik = col_ptr-column_numbers; + global_entry(global_index_ik) *= 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 + bool left_of_diagonal = true; + for (const int * j = col_ptr+1; + j<&column_numbers[rowstart_indices[row+1]]; + ++j) + { + // note: this inner loop could + // be made considerable faster + // if we consulted the row + // with number *col_ptr, + // instead of always asking + // sparsity(*col_ptr,*j), + // since we traverse this + // row linearly. I just didn't + // have the time to figure out + // the details. + + // check whether we have just + // traversed the diagonal + // + // note that j should never point + // to the diagonal itself! + if (left_of_diagonal && (*j > static_cast(row))) + { + Assert (*j != static_cast(row), ExcInternalError()); + + left_of_diagonal = false; + // a[i,i] -= a[i,k]*a[k,i] + const int global_index_ki = sparsity(*col_ptr,row); + + if (global_index_ki != -1) + diag_element(row) -= global_entry(global_index_ik) * + global_entry(global_index_ki); + + }; + + const int global_index_ij = j - &column_numbers[0], + global_index_kj = sparsity(*col_ptr,*j); + if ((global_index_ij != -1) && + (global_index_kj != -1)) + global_entry(global_index_ij) -= global_entry(global_index_ik) * + global_entry(global_index_kj); + }; + }; + }; + +/* + OLD CODE, rather crude first implementation with an algorithm taken + from 'W. Hackbusch, G. Wittum: Incomplete Decompositions (ILU)- + Algorithms, Theory, and Applications', page 6. + + for (unsigned int k=0; k +template +void SparseILU::apply_decomposition (Vector &dst, + const Vector &src) const +{ + Assert (dst.size() == src.size(), ExcSizeMismatch(dst.size(), src.size())); + Assert (dst.size() == m(), ExcSizeMismatch(dst.size(), m())); + + const unsigned int N=dst.size(); + const unsigned int * const rowstart_indices + = get_sparsity_pattern().get_rowstart_indices(); + const int * const column_numbers + = get_sparsity_pattern().get_column_numbers(); + // solve LUx=b in two steps: + // first Ly = b, then + // Ux = y + // + // first a forward solve. since + // the diagonal values of L are + // one, there holds + // y_i = b_i + // - sum_{j=0}^{i-1} L_{ij}y_j + // we split the y_i = b_i off and + // perform it at the outset of the + // loop + dst = src; + for (unsigned int row=0; row(row)); + + for (const int * col=rowstart; col!=first_after_diagonal; ++col) + dst(row) -= global_entry (col-column_numbers) * dst(*col); + }; + + // now the backward solve. same + // procedure, but we need not set + // dst before, since this is already + // done. + // + // note that we need to scale now, + // since the diagonal is not zero + // now + for (int row=N-1; row>=0; --row) + { + // get end of this row + const int * const rowend = &column_numbers[rowstart_indices[row+1]]; + // find the position where the part + // right of the diagonal starts + const int * const first_after_diagonal + = lower_bound (&column_numbers[rowstart_indices[row]+1], + &column_numbers[rowstart_indices[row+1]], + row); + + for (const int * col=first_after_diagonal; col!=rowend; ++col) + dst(row) -= global_entry (col-column_numbers) * dst(*col); + + // scale by the diagonal element. + // note that the diagonal element + // was stored inverted + dst(row) *= diag_element(row); + }; +}; + + + + +/*---------------------------- matrix_decompositions.templates.h ---------------------------*/ +/* end of #ifndef __matrix_decompositions_templates_H */ +#endif +/*---------------------------- matrix_decompositions.templates.h ---------------------------*/ diff --git a/deal.II/lac/source/sparse_ilu.cc b/deal.II/lac/source/sparse_ilu.cc new file mode 100644 index 0000000000..7d48252d6b --- /dev/null +++ b/deal.II/lac/source/sparse_ilu.cc @@ -0,0 +1,27 @@ +/* $Id$ */ + +#include + + +// explicit instantiations +template class SparseILU; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; + + +template class SparseILU; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; + -- 2.39.5