]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
preparations for ILU; comments for backward/forward
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Apr 1999 09:07:50 +0000 (09:07 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 15 Apr 1999 09:07:50 +0000 (09:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@1142 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h

index 7e3fc0ea23545596312672e235e1ad5bb60f78eb..ff677f8e0045f630a3b7bf3a0c4f83e4b40eb333 100644 (file)
@@ -455,13 +455,22 @@ class FullMatrix
     double residual (Vector<number2>& w, const Vector<number2>& v, const Vector<number3>& 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<typename number2>
     void forward (Vector<number2>& dst, const Vector<number2>& src) const;
 
                                     /**
-                                     *  Inversion of upper triangle .
+                                     * Backward elimination of upper triangle.
+                                     * @see forward
                                      */
     template<typename number2>
     void backward (Vector<number2>& dst, const Vector<number2>& src) const;
index df5b23e574c63d6ec6c663af76238a6fbd1c2f28..1289e6743ef293df62aed0a365f0275d3cda0e18 100644 (file)
@@ -422,14 +422,14 @@ double FullMatrix<number>::residual (Vector<number2>& dst, const Vector<number2>
 
 template <typename number>
 template <typename number2>
-void FullMatrix<number>::forward (Vector<number2>& dst, const Vector<number2>& src) const
+void FullMatrix<number>::forward (Vector<number2>& dst,
+                                 const Vector<number2>& 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()<n() ? m() : n());
+  unsigned int nu = ( (m()<n()) ? m() : n());
   double s;
   for (i=0; i<nu; ++i)
     {
@@ -443,7 +443,8 @@ void FullMatrix<number>::forward (Vector<number2>& dst, const Vector<number2>& s
 
 template <typename number>
 template <typename number2>
-void FullMatrix<number>::backward (Vector<number2>& dst, const Vector<number2>& src) const
+void FullMatrix<number>::backward (Vector<number2>& dst,
+                                  const Vector<number2>& src) const
 {
   unsigned int j;
   unsigned int nu = (m()<n() ? m() : n());
index a8d02d5a5e49e1a6a9a243d9000858a2ab8e665c..4c48052140222f26481fe0832cd9f6dbaebe0ff9 100644 (file)
@@ -69,8 +69,8 @@ class SparseMatrixStruct
 
                                     /**
                                      * Initialize a rectangular matrix with
-                                     * #m# rows and #n# columns,
-                                     * with at most #max_per_row#
+                                     * #m# rows and #n# columns.
+                                     * The matrix may contain at most #max_per_row#
                                      * nonzero entries per row.
                                      */
     SparseMatrixStruct (const unsigned int m,
@@ -85,6 +85,18 @@ class SparseMatrixStruct
     SparseMatrixStruct (const unsigned int n,
                        const unsigned int max_per_row);
 
+                                    /**
+                                     * Make a copy with extra off-diagonals.
+                                     *
+                                     * This constructs objects intended for
+                                     * the application of the ILU(n)-method.
+                                     * Therefore, additional to the original
+                                     * entry structure, space for #extra_cols#
+                                     * side-diagonals is provided.
+                                     */
+    SparseMatrixStruct(const SparseMatrixStruct& original,
+                      unsigned int extra_cols);
+    
                                     /**
                                      * Destructor.
                                      */
@@ -555,7 +567,36 @@ class SparseMatrix
                                      * #this#.
                                      */
     template <typename somenumber>
-    SparseMatrix<number> & copy_from (const SparseMatrix<somenumber> &);
+    SparseMatrix<number> & copy_from (const SparseMatrix<somenumber> &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<float> 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<float> ilu(ilustruct);
+                                     * ilu.ILU(source);
+                                     * \end{verbatim}
+                                     *
+                                     * After generating the ILU-decomposition,
+                                     * it can be applied to a vector by
+                                     * #backward_forward#.
+                                     */
+    template <typename somenumber>
+    SparseMatrix<number> & ILU (const SparseMatrix<somenumber> &source);
 
                                     /**
                                      * Add #matrix# scaled by #factor# to this
@@ -571,7 +612,8 @@ class SparseMatrix
                                      * data type of this matrix.
                                      */
     template <typename somenumber>
-    void add_scaled (const number factor, const SparseMatrix<somenumber> &matrix);
+    void add_scaled (const number factor,
+                    const SparseMatrix<somenumber> &matrix);
     
                                     /**
                                      * Return the value of the entry (i,j).
@@ -640,7 +682,13 @@ class SparseMatrix
     template <typename somenumber>
     void Tvmult (Vector<somenumber>& dst, const Vector<somenumber>& src) const;
   
-
+                                    /**
+                                     * Do backward-forward solution of a
+                                     * previously generated ILU.
+                                     */
+    template <typename somenumber>
+    void backward_forward (Vector<somenumber>& 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
index de6d70658c3b46106bfed9bd1f944111b1fb0864..4688022c34b4e396d4f4656ac74e064e5bb418ae 100644 (file)
@@ -21,11 +21,11 @@ template <typename number>
 SparseMatrix<number>::SparseMatrix () :
                cols(0),
                val(0),
-               max_len(0)
+               max_len(0),
+               is_ilu(false)
 {};
 
 
-
 template <typename number>
 SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
                cols(0),
@@ -41,7 +41,10 @@ SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
 
 template <typename number>
 SparseMatrix<number>::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<number>::reinit ()
 
   if (val)
     fill_n (&val[0], cols->vec_len, 0);
+
+  is_ilu = false;
 }
 
 
 
 template <typename number>
 void
-SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity) {
+SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity)
+{
   cols = &sparsity;
   reinit ();
 };
@@ -95,11 +101,13 @@ SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity) {
 
 template <typename number>
 void
-SparseMatrix<number>::clear () {
+SparseMatrix<number>::clear ()
+{
   cols = 0;
   if (val) delete[] val;
   val = 0;
   max_len = 0;
+  is_ilu = false;
 };
 
 
@@ -118,7 +126,8 @@ SparseMatrix<number>::empty () const
 
 template <typename number>
 unsigned int
-SparseMatrix<number>::n_nonzero_elements () const {
+SparseMatrix<number>::n_nonzero_elements () const
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   return cols->n_nonzero_elements ();
 };
@@ -128,7 +137,8 @@ SparseMatrix<number>::n_nonzero_elements () const {
 template <typename number>
 template <typename somenumber>
 SparseMatrix<number> &
-SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix) {
+SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix)
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
   Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
@@ -140,6 +150,7 @@ SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix) {
   while (val_ptr != end_ptr)
     *val_ptr++ = *matrix_ptr++;
   
+  is_ilu = false;
   return *this;
 };
 
@@ -149,7 +160,8 @@ template <typename number>
 template <typename somenumber>
 void
 SparseMatrix<number>::add_scaled (const number factor,
-                                 const SparseMatrix<somenumber> &matrix) {
+                                 const SparseMatrix<somenumber> &matrix)
+{
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
   Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
@@ -160,6 +172,7 @@ SparseMatrix<number>::add_scaled (const number factor,
 
   while (val_ptr != end_ptr)
     *val_ptr++ += factor * *matrix_ptr++;
+  is_ilu = false;
 };
 
 
@@ -314,7 +327,8 @@ SparseMatrix<number>::residual (Vector<somenumber>& dst,
 template <typename number>
 template <typename somenumber>
 void
-SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
+SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst,
+                                          const Vector<somenumber>& src,
                                           const number om) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
@@ -339,7 +353,8 @@ SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst, const Vector
 template <typename number>
 template <typename somenumber>
 void
-SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
+                                        const Vector<somenumber>& src,
                                         const number om) const
 {
                                   // to understand how this function works

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.