]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add sparse ILU decompoition.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Apr 1999 19:00:24 +0000 (19:00 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Apr 1999 19:00:24 +0000 (19:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@1227 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_ilu.h [new file with mode: 0644]
deal.II/lac/include/lac/sparse_ilu.templates.h [new file with mode: 0644]
deal.II/lac/source/sparse_ilu.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h
new file mode 100644 (file)
index 0000000..6a33d84
--- /dev/null
@@ -0,0 +1,202 @@
+/*----------------------------   matrix_decompositions.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __sparse_ilu_H
+#define __sparse_ilu_H
+/*----------------------------   sparse_ilu.h     ---------------------------*/
+
+
+#include <lac/sparsematrix.h>
+
+
+
+
+/**
+ * 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<double> ilu (ilu_sparsity);
+ *   ilu.decompose (global_matrix);
+ *   global_matrix.set_preconditioner (ilu);
+ *
+ *   somesolver.solve (A, x, f,
+ *                     PreconditionUseMatrix<SparseILU<double>,Vector<double> >
+ *                           (ilu,&SparseILU<double>::template apply_decomposition<double>));
+ * \end{verbatim}
+ *
+ *
+ * @author Wolfgang Bangerth, 1999, based on a similar implementation by Malte Braack
+ */
+template <typename number>
+class SparseILU : protected SparseMatrix<number>
+{
+  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 <typename somenumber>
+    void decompose (const SparseMatrix<somenumber> &matrix,
+                   const double                    strengthen_diagonal=0.);
+
+                                    /**
+                                     * Apply the incomplete decomposition,
+                                     * i.e. do one forward-backward step
+                                     * $dst=(LU)^{-1}src$.
+                                     */
+    template <typename somenumber>
+    void apply_decomposition (Vector<somenumber>       &dst,
+                             const Vector<somenumber> &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 (file)
index 0000000..3b1de72
--- /dev/null
@@ -0,0 +1,339 @@
+/*----------------------------   matrix_decompositions.templates.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __matrix_decompositions_templates_H
+#define __matrix_decompositions_templates_H
+/*----------------------------   matrix_decompositions.templates.h     ---------------------------*/
+
+
+#include <lac/vector.h>
+#include <lac/sparse_ilu.h>
+
+#include <algorithm>
+#include <cmath>
+
+
+
+template <typename number>
+SparseILU<number>::SparseILU () 
+{};
+
+
+
+template <typename number>
+SparseILU<number>::SparseILU (const SparseMatrixStruct &sparsity) :
+               SparseMatrix<number> (sparsity)
+{};
+
+
+
+template <typename number>
+void SparseILU<number>::reinit ()
+{
+  SparseMatrix<number>::reinit ();
+};
+
+
+
+template <typename number>
+void SparseILU<number>::reinit (const SparseMatrixStruct &sparsity)
+{
+  SparseMatrix<number>::reinit (sparsity);
+};
+
+
+
+
+template <typename number>
+template <typename somenumber>
+void SparseILU<number>::decompose (const SparseMatrix<somenumber> &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<m(); ++row)
+       for (const int * col = &column_numbers[rowstart_indices[row]];
+            col != &column_numbers[rowstart_indices[row+1]]; ++col)
+         set (row, *col, matrix.global_entry(col-column_numbers));
+    };
+
+  if (strengthen_diagonal > 0)
+    for (unsigned int row=0; row<m(); ++row)
+      {
+                                        // get the length of the row
+                                        // (without the diagonal element)
+       const unsigned int rowlength = get_sparsity_pattern().get_rowstart_indices()[row+1]
+                                      -get_sparsity_pattern().get_rowstart_indices()[row]
+                                      -1;
+       
+                                        // get the global index of the first
+                                        // non-diagonal element in this row
+       const unsigned int rowstart
+         = get_sparsity_pattern().get_rowstart_indices()[row] + 1;
+       number * const diagonal_element = &global_entry(rowstart-1);
+
+       number rowsum = 0;
+       for (unsigned int global_index=rowstart;
+            global_index<rowstart+rowlength; ++global_index)
+         rowsum += fabs(global_entry(global_index));
+
+       *diagonal_element += strengthen_diagonal * rowsum;
+      };
+  
+         
+                                  // now work only on this
+                                  // matrix
+  const SparseMatrixStruct &sparsity = get_sparsity_pattern();
+  const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices();
+  const 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]
+*/
+
+  
+                                  // i := row
+  for (unsigned int row=1; row<m(); ++row)
+    {
+                                      // invert diagonal element of the
+                                      // previous row. this is a hack,
+                                      // which is possible since this
+                                      // element is not needed any more
+                                      // in the process of decomposition
+                                      // and since it makes the backward
+                                      // step when applying the decomposition
+                                      // significantly faster
+      global_entry (rowstart_indices[row-1])
+       = 1./global_entry (rowstart_indices[row-1]);
+
+                                      // let k run over all lower-left
+                                      // elements of row i; skip
+                                      // diagonal element at start
+      const int * first_of_row
+       = &column_numbers[rowstart_indices[row]+1];
+      const int * first_after_diagonal
+       = lower_bound (&column_numbers[rowstart_indices[row]+1],
+                      &column_numbers[rowstart_indices[row+1]],
+                      static_cast<int>(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<int>(row)))
+               {
+                 Assert (*j != static_cast<int>(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<m()-1; ++k)
+    for (unsigned int i=k+1; i<m(); ++i)
+      {
+                                        // get the global index
+                                        // of the element (i,k)
+       const int global_index_ik = get_sparsity_pattern()(i,k);
+
+                                        // if this element is zero,
+                                        // then we continue with the
+                                        // next i, since e would be
+                                        // zero and nothing would happen
+                                        // in this loop
+       if (global_index_ik == -1)
+         continue;
+       
+       const number e = global_entry(global_index_ik) / diag_element(k);
+       global_entry(global_index_ik) = e;
+
+       for (unsigned int j=k+1; j<m(); ++j)
+         {
+                                            // find out about a_kj
+                                            // if this does not exist,
+                                            // then the updates within
+                                            // this innermost loop would
+                                            // be zero, invariable of the
+                                            // fact of whether a_ij is a
+                                            // nonzero or a zero element
+           const int global_index_kj = get_sparsity_pattern()(k,j);
+           if (global_index_kj == -1)
+             continue;
+
+           const int global_index_ij = get_sparsity_pattern()(i,j);
+           if (global_index_ij != -1)
+             global_entry(global_index_ij) -= e*global_entry(global_index_kj);
+           else
+             diag_element(i) -= e*global_entry(global_index_kj);
+         };
+      };
+*/      
+};
+
+
+
+
+template <typename number>
+template <typename somenumber>
+void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
+                                            const Vector<somenumber> &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<N; ++row)
+    {
+                                      // get start of this row. skip the
+                                      // diagonal element
+      const int * const rowstart = &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 (rowstart,
+                      &column_numbers[rowstart_indices[row+1]],
+                      static_cast<int>(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 (file)
index 0000000..7d48252
--- /dev/null
@@ -0,0 +1,27 @@
+/* $Id$ */
+
+#include <lac/sparse_ilu.templates.h>
+
+
+// explicit instantiations
+template class SparseILU<double>;
+template void SparseILU<double>::decompose (const SparseMatrix<double> &,
+                                           const double);
+template void SparseILU<double>::apply_decomposition (Vector<double> &,
+                                                     const Vector<double> &) const;
+template void SparseILU<double>::decompose (const SparseMatrix<float> &,
+                                           const double);
+template void SparseILU<double>::apply_decomposition (Vector<float> &,
+                                                     const Vector<float> &) const;
+
+
+template class SparseILU<float>;
+template void SparseILU<float>::decompose (const SparseMatrix<double> &,
+                                          const double);
+template void SparseILU<float>::apply_decomposition (Vector<double> &,
+                                                    const Vector<double> &) const;
+template void SparseILU<float>::decompose (const SparseMatrix<float> &,
+                                          const double);
+template void SparseILU<float>::apply_decomposition (Vector<float> &,
+                                                    const Vector<float> &) const;
+

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.