From: guido Date: Thu, 15 Apr 1999 09:07:50 +0000 (+0000) Subject: preparations for ILU; comments for backward/forward X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2674ec8ee38844b861dae81b164e65d378a96923;p=dealii-svn.git preparations for ILU; comments for backward/forward git-svn-id: https://svn.dealii.org/trunk@1142 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index 7e3fc0ea23..ff677f8e00 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -455,13 +455,22 @@ class FullMatrix double residual (Vector& w, const Vector& v, const Vector& b) const; /** - * Inversion of lower triangle . + * Forward elimination of lower triangle. + * Inverts the lower triangle of a + * quadratic matrix. + * + * If the matrix has more columns than rows, + * this function only operates on the left + * square submatrix. If there are more rows, + * the upper square part of the matrix + * is considered */ template void forward (Vector& dst, const Vector& src) const; /** - * Inversion of upper triangle . + * Backward elimination of upper triangle. + * @see forward */ template void backward (Vector& dst, const Vector& src) const; diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index df5b23e574..1289e6743e 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -422,14 +422,14 @@ double FullMatrix::residual (Vector& dst, const Vector template template -void FullMatrix::forward (Vector& dst, const Vector& src) const +void FullMatrix::forward (Vector& dst, + const Vector& src) const { - Assert(n() == m(), ExcNotQuadratic()); - Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); + Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); unsigned int i,j; - unsigned int nu = (m()::forward (Vector& dst, const Vector& s template template -void FullMatrix::backward (Vector& dst, const Vector& src) const +void FullMatrix::backward (Vector& dst, + const Vector& src) const { unsigned int j; unsigned int nu = (m() - SparseMatrix & copy_from (const SparseMatrix &); + SparseMatrix & copy_from (const SparseMatrix &source); + + /** + * Generate ILU. + * The matrix will entries contain the + * incomplete LU-factorization of + * the matrix #source#. Having a matrix + * #source# and a matrix structure object + * #struct#, the code for generating + * an ILU factorization reads + * \begin{verbatim} + * SparseMatrix ilu(struct); + * ilu.ILU(source); + * \end{verbatim} + * + * If additional side diagonals are + * needed (ILU(n)-algorithm), you have to + * construct a second matrix structure: + * \begin{verbatim} + * SparseMatrixStruct ilustruct(struct,n); + * SparseMatrix ilu(ilustruct); + * ilu.ILU(source); + * \end{verbatim} + * + * After generating the ILU-decomposition, + * it can be applied to a vector by + * #backward_forward#. + */ + template + SparseMatrix & ILU (const SparseMatrix &source); /** * Add #matrix# scaled by #factor# to this @@ -571,7 +612,8 @@ class SparseMatrix * data type of this matrix. */ template - void add_scaled (const number factor, const SparseMatrix &matrix); + void add_scaled (const number factor, + const SparseMatrix &matrix); /** * Return the value of the entry (i,j). @@ -640,7 +682,13 @@ class SparseMatrix template void Tvmult (Vector& dst, const Vector& src) const; - + /** + * Do backward-forward solution of a + * previously generated ILU. + */ + template + void backward_forward (Vector& v); + /** * Return the norm of the vector #v# with * respect to the norm induced by this @@ -817,12 +865,18 @@ class SparseMatrix /** * Exception */ + DeclException0 (ExcNoILU); + /** + * Exception + */ DeclException0 (ExcInvalidConstructorCall); private: const SparseMatrixStruct * cols; number* val; unsigned int max_len; + bool is_ilu; + // make all other sparse matrices // friends diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index de6d70658c..4688022c34 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -21,11 +21,11 @@ template SparseMatrix::SparseMatrix () : cols(0), val(0), - max_len(0) + max_len(0), + is_ilu(false) {}; - template SparseMatrix::SparseMatrix (const SparseMatrix &m) : cols(0), @@ -41,7 +41,10 @@ SparseMatrix::SparseMatrix (const SparseMatrix &m) : template SparseMatrix::SparseMatrix (const SparseMatrixStruct &c) - : cols(&c), val(0), max_len(0) + : cols(&c), + val(0), + max_len(0), + is_ilu(false) { reinit(); }; @@ -80,13 +83,16 @@ SparseMatrix::reinit () if (val) fill_n (&val[0], cols->vec_len, 0); + + is_ilu = false; } template void -SparseMatrix::reinit (const SparseMatrixStruct &sparsity) { +SparseMatrix::reinit (const SparseMatrixStruct &sparsity) +{ cols = &sparsity; reinit (); }; @@ -95,11 +101,13 @@ SparseMatrix::reinit (const SparseMatrixStruct &sparsity) { template void -SparseMatrix::clear () { +SparseMatrix::clear () +{ cols = 0; if (val) delete[] val; val = 0; max_len = 0; + is_ilu = false; }; @@ -118,7 +126,8 @@ SparseMatrix::empty () const template unsigned int -SparseMatrix::n_nonzero_elements () const { +SparseMatrix::n_nonzero_elements () const +{ Assert (cols != 0, ExcMatrixNotInitialized()); return cols->n_nonzero_elements (); }; @@ -128,7 +137,8 @@ SparseMatrix::n_nonzero_elements () const { template template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &matrix) { +SparseMatrix::copy_from (const SparseMatrix &matrix) +{ Assert (cols != 0, ExcMatrixNotInitialized()); Assert (val != 0, ExcMatrixNotInitialized()); Assert (cols == matrix.cols, ExcDifferentSparsityPatterns()); @@ -140,6 +150,7 @@ SparseMatrix::copy_from (const SparseMatrix &matrix) { while (val_ptr != end_ptr) *val_ptr++ = *matrix_ptr++; + is_ilu = false; return *this; }; @@ -149,7 +160,8 @@ template template void SparseMatrix::add_scaled (const number factor, - const SparseMatrix &matrix) { + const SparseMatrix &matrix) +{ Assert (cols != 0, ExcMatrixNotInitialized()); Assert (val != 0, ExcMatrixNotInitialized()); Assert (cols == matrix.cols, ExcDifferentSparsityPatterns()); @@ -160,6 +172,7 @@ SparseMatrix::add_scaled (const number factor, while (val_ptr != end_ptr) *val_ptr++ += factor * *matrix_ptr++; + is_ilu = false; }; @@ -314,7 +327,8 @@ SparseMatrix::residual (Vector& dst, template template void -SparseMatrix::precondition_Jacobi (Vector& dst, const Vector& src, +SparseMatrix::precondition_Jacobi (Vector& dst, + const Vector& src, const number om) const { Assert (cols != 0, ExcMatrixNotInitialized()); @@ -339,7 +353,8 @@ SparseMatrix::precondition_Jacobi (Vector& dst, const Vector template template void -SparseMatrix::precondition_SSOR (Vector& dst, const Vector& src, +SparseMatrix::precondition_SSOR (Vector& dst, + const Vector& src, const number om) const { // to understand how this function works