]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add Stephens implementation of sparse decompositions.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Dec 2002 17:28:05 +0000 (17:28 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 Dec 2002 17:28:05 +0000 (17:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@6792 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/lac/include/lac/sparse_decomposition.h [new file with mode: 0644]
deal.II/lac/include/lac/sparse_decomposition.templates.h [new file with mode: 0644]
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.h.x [new file with mode: 0644]
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/include/lac/sparse_ilu.templates.h.x [new file with mode: 0644]
deal.II/lac/include/lac/sparse_mic.h [new file with mode: 0644]
deal.II/lac/include/lac/sparse_mic.templates.h [new file with mode: 0644]
deal.II/lac/source/sparse_decomposition.cc [new file with mode: 0644]
deal.II/lac/source/sparse_ilu.cc
deal.II/lac/source/sparse_ilu.cc.x [new file with mode: 0644]
deal.II/lac/source/sparse_mic.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/sparse_decomposition.h b/deal.II/lac/include/lac/sparse_decomposition.h
new file mode 100644 (file)
index 0000000..24b1364
--- /dev/null
@@ -0,0 +1,327 @@
+//----------------------  sparse_decomposition.h  ---------------------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------  sparse_decomposition.h  ---------------------------
+#ifndef __deal2__sparse_decomposition_h
+#define __deal2__sparse_decomposition_h
+
+#include <base/config.h>
+#include <lac/sparse_matrix.h>
+
+#include <cmath>
+
+/**
+ * Abstract base class for sparse LU decompositions of a sparse matrix
+ * into another sparse matrix.
+ *
+ * The decomposition is stored as a sparse matrix, for
+ * which the user has to give a sparsity pattern and which is why this
+ * class is derived from the @p{SparseMatrix}. Since it is not a matrix in
+ * the usual sense, the derivation is @p{protected} rather than @p{public}.
+ *
+ * @sect3{Fill-in}
+ *
+ * The sparse LU decompositions are 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 @p{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 @p{SparsityPattern} class which allows the addition
+ * of side-diagonals to a given sparsity structure.
+ *
+ *
+ * @sect3{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 (@p{lu_sparsity} is some sparsity pattern to be used
+ * for the decomposition, which you have to create beforehand):
+ * @begin{verbatim}
+ *   SparseLUImplementation<double> lu (lu_sparsity);
+ *   lu.decompose (global_matrix);
+ *
+ *   somesolver.solve (A, x, f,  lu);
+ * @end{verbatim}
+ *
+ * @sect2{State management}
+ *
+ * In order to prevent users from applying decompositions before the
+ * decomposition itself has been built, and to introduce some
+ * optimization of common "sparse idioms", this class introduces a
+ * simple state management.  A SparseLUdecomposition instance is
+ * considered @p{not decomposed} if the decompose method has not yet
+ * been invoked since the last time the underlying @ref{SparseMatrix}
+ * had changed. The underlying sparse matrix is considered changed
+ * when one of this class reinit methods, constructors or destructors
+ * are invoked.  The @p{not decomposed} state is indicated by a false
+ * value returned by @p{is_decomposed} method.  It is illegal to apply
+ * this decomposition (@p{vmult} method) in not decomposed state; in
+ * this case, the @p{vmult} method throws an @p{ExcInvalidState}
+ * exception. This object turns into decomposed state immediately
+ * after its @p{decompose} method is invoked. The @p{decomposed}
+ * state is indicated by true value returned by @p{is_decomposed}
+ * method. It is legal to apply this decomposition (@p{vmult} method) in
+ * decompoed state.
+ *
+ *
+ * @sect2{Particular implementations}
+ *
+ * It is enough to override the @p{decompose} and @p{vmult} methods to
+ * implement particular LU decompositions, like the true LU, or the
+ * Cholesky decomposition. Additionally, if that decomposition needs
+ * fine tuned diagonal strengthening on a per row basis, it may override the
+ * @p{get_strengthen_diagonal} method. You should invoke the non-abstract
+ * base class method to employ the state management. Implementations
+ * may choose more restrictive definition of what is legal or illegal
+ * state; but they must conform to the @p{is_decomposed} method
+ * specification above.
+ *
+ * If an exception is thrown by method other than @p{vmult}, this
+ * object may be left in an inconsistent state.
+ *
+ * @author Stephen "Cheffo" Kolaroff, 2002, based on SparseILU implementation by Wolfgang Bangerth
+ */
+template <typename number>
+class SparseLUDecomposition : 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
+                                     * @p{reinit(SparsityPattern)}.
+                                     */
+
+    SparseLUDecomposition ();
+
+                                    /**
+                                     * 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 @p{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 @p{reinit} is not
+                                     * called with a new sparsity
+                                     * structure.
+                                     */
+    SparseLUDecomposition (const SparsityPattern& sparsity);
+
+                                     /**
+                                      * Destruction.
+                                      */
+    virtual ~SparseLUDecomposition ();
+
+                                    /**
+                                     * Reinitialize the object but
+                                     * keep to the sparsity pattern
+                                     * previously used.  This may be
+                                     * necessary if you @p{reinit}'d
+                                     * the sparsity structure and
+                                     * want to update the size of the
+                                     * matrix.
+                                     *
+                                     * After this method is invoked,
+                                     * this object is out of synch
+                                     * (not decomposed state).
+                                     *
+                                     * This function only releases
+                                     * some memory and calls the
+                                     * respective function of the
+                                     * base class.
+                                     */
+    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 only releases
+                                     * some memory and calls the
+                                     * respective function of the
+                                     * base class.
+                                     */
+    void reinit (const SparsityPattern &sparsity);
+
+                                    /**
+                                     * Perform the sparse LU
+                                     * factorization of the given
+                                     * matrix. After this method
+                                     * invokation, and before
+                                     * consequtive reinit invokation
+                                     * this object is in decomposed
+                                     * state.
+                                     *
+                                     * 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.);
+
+                                     /**
+                                      * Determines if this object is
+                                      * in synch with the underlying
+                                      * @ref{SparsityPattern}.
+                                      */ 
+    virtual bool is_decomposed () const;       
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     */
+    virtual unsigned int memory_consumption () const;
+
+                                     /**
+                                      * Exception
+                                      */
+    DeclException1 (ExcInvalidStrengthening,
+                   double,
+                   << "The strengthening parameter " << arg1
+                   << " is not greater or equal than zero!");
+
+                                     /**
+                                      * Exception. Indicates violation
+                                      * of a @p{state rule}.
+                                      */
+    DeclException0 (ExcInvalidState);
+
+  protected:
+                                     /**
+                                      * Copies the passed SparseMatrix
+                                      * onto this object. This
+                                      * object's sparsity pattern
+                                      * remains unchanged.
+                                      */
+    template<typename somenumber>
+    void copy_from (const SparseMatrix<somenumber>& matrix);
+
+                                     /**
+                                      * Performs the strengthening
+                                      * loop. For each row calculates
+                                      * the sum of absolute values of
+                                      * its elements, determines the
+                                      * strengthening factor (through
+                                      * @p{get_strengthen_diagonal})
+                                      * sf and multiplies the diagonal
+                                      * entry with @p{sf+1}.
+                                      */
+    virtual void strengthen_diagonal_impl ();
+
+                                     /**
+                                      * In the decomposition phase,
+                                      * computes a strengthening
+                                      * factor for the diagonal entry
+                                      * in row @p{row} with sum of
+                                      * absolute values of its
+                                      * elements @p{rowsum}.<br> Note:
+                                      * The default implementation in
+                                      * @ref{SparseLUDecomposition}
+                                      * returns
+                                      * @p{strengthen_diagonal}'s
+                                      * value.
+                                      */
+    virtual number get_strengthen_diagonal(const number rowsum, const unsigned int row) const;
+    
+                                     /**
+                                      * State flag. If not in
+                                      * @em{decomposed} state, it is
+                                      * unlegal to apply the
+                                      * decomposition.  This flag is
+                                      * cleared when the underlaying
+                                      * @ref{SparseMatrix}
+                                      * @ref{SparsityPattern} is
+                                      * changed, and set by
+                                      * @p{decompose}.
+                                      */
+    bool decomposed;
+
+                                     /**
+                                      * The default strenghtening
+                                      * value, returned by
+                                      * @p{get_strengthen_diagonal}.
+                                      */
+    double  strengthen_diagonal;
+
+                                     /**
+                                      * For every row in the
+                                      * underlying
+                                      * @ref{SparsityPattern}, this
+                                      * array contains a pointer
+                                      * to the row's first
+                                      * afterdiagonal entry. Becomes
+                                      * available after invokation of
+                                      * @p{decompose}.
+                                      */
+    std::vector<const unsigned int*> prebuilt_lower_bound;
+    
+  private:
+                                     /**
+                                      * Fills the
+                                      * @ref{prebuilt_lower_bound}
+                                      * array.
+                                      */
+    void prebuild_lower_bound ();
+    
+};
+
+
+
+template <typename number>
+inline number
+SparseLUDecomposition<number>::
+get_strengthen_diagonal(const number rowsum, const unsigned int row) const
+{
+  return strengthen_diagonal;
+};
+
+
+
+template <typename number>
+inline bool
+SparseLUDecomposition<number>::is_decomposed () const
+{
+  return decomposed;
+}
+
+#endif // __deal2__sparse_decomposition_h
diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h
new file mode 100644 (file)
index 0000000..52911d6
--- /dev/null
@@ -0,0 +1,170 @@
+//---------------------  sparse_decomposition.templates.h  ----------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------  sparse_decomposition.templates.h  ----------------
+
+#include <base/memory_consumption.h>
+#include <lac/sparse_decomposition.h>
+#include <algorithm>
+
+
+template<typename number>
+SparseLUDecomposition<number>::SparseLUDecomposition()
+                :
+                SparseMatrix<number>(),
+                decomposed(false)
+{}
+
+
+
+template<typename number>
+SparseLUDecomposition<number>::
+SparseLUDecomposition (const SparsityPattern& sparsity) :
+                SparseMatrix<number>(sparsity),
+                decomposed(false)
+{}
+
+
+
+template<typename number>
+SparseLUDecomposition<number>::~SparseLUDecomposition()
+{}
+
+
+
+template<typename number>
+template<typename somenumber>
+void
+SparseLUDecomposition<number>::
+decompose (const SparseMatrix<somenumber> &matrix,
+           const double                    strengthen_diagonal)
+{
+  decomposed = false;
+  
+  this->strengthen_diagonal = strengthen_diagonal;
+  prebuild_lower_bound ();
+  copy_from (matrix);
+  decomposed = true;
+}
+
+
+
+template <typename number>
+void
+SparseLUDecomposition<number>::reinit ()
+{
+  decomposed = false;
+  if (true)
+    {
+      std::vector<const unsigned int*> tmp;
+      tmp.swap (prebuilt_lower_bound);
+    };
+  SparseMatrix<number>::reinit ();
+}
+
+
+
+template <typename number>
+void SparseLUDecomposition<number>::reinit (const SparsityPattern& sparsity)
+{
+  decomposed = false;
+  if (true)
+    {
+      std::vector<const unsigned int*> tmp;
+      tmp.swap (prebuilt_lower_bound);
+    };
+  SparseMatrix<number>::reinit (sparsity);
+}
+
+
+
+template<typename number>
+void
+SparseLUDecomposition<number>::prebuild_lower_bound()
+{
+  const unsigned int* const column_numbers = get_sparsity_pattern().get_column_numbers();
+  const unsigned int* const rowstart_indices = get_sparsity_pattern().get_rowstart_indices();
+  const unsigned int N = m();
+
+  prebuilt_lower_bound.resize (N);
+
+  for(unsigned int row=0; row<N; row++) {
+    prebuilt_lower_bound[row]
+      = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
+                          &column_numbers[rowstart_indices[row+1]],
+                          row);
+  }
+}
+
+template <typename number>
+template <typename somenumber>
+void
+SparseLUDecomposition<number>::copy_from (const SparseMatrix<somenumber>& matrix)
+{
+                                   // preset the elements
+  std::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 unsigned int * const column_numbers
+    = matrix.get_sparsity_pattern().get_column_numbers();
+
+  for (unsigned int row=0; row<m(); ++row)
+    for (const unsigned int * col = &column_numbers[rowstart_indices[row]];
+         col != &column_numbers[rowstart_indices[row+1]]; ++col)
+      set (row, *col, matrix.global_entry(col-column_numbers));
+}
+
+
+
+template <typename number>
+void
+SparseLUDecomposition<number>::strengthen_diagonal_impl ()
+{
+  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 += std::fabs(global_entry(global_index));
+
+      *diagonal_element += get_strengthen_diagonal (rowsum, row)  * rowsum;
+    }
+}
+
+
+
+template <typename number>
+inline unsigned int
+SparseLUDecomposition<number>::memory_consumption () const
+{
+  unsigned int
+    res = (SparseMatrix<number>::memory_consumption () +
+           MemoryConsumption::memory_consumption(prebuilt_lower_bound));
+  return res;
+}
+
+
index 00359dcdc94632dbe7ea84af0a2a986dac54040d..e4a37461141c065c3d83bf39f0e54b8875f4f640 100644 (file)
@@ -1,8 +1,6 @@
 //----------------------------  sparse_ilu.h  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 //    further information on this license.
 //
 //----------------------------  sparse_ilu.h  ---------------------------
+
 #ifndef __deal2__sparse_ilu_h
 #define __deal2__sparse_ilu_h
 
 
 #include <base/config.h>
 #include <lac/sparse_matrix.h>
+#include <lac/sparse_decomposition.h>
 
 
 /**
  *      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 @p{SparseMatrix}. Since it is not a matrix in
- * the usual sense, the derivation is @p{protected} rather than @p{public}.
- *
- * Note that in the algorithm given, the lower left part of the matrix base
- * class is used to store the @p{L} part of the decomposition, while 
- * the upper right part is used to store @p{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).
- *
- *
- * @sect3{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 @p{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 @p{SparsityPattern} class which allows the addition
- * of side-diagonals to a given sparsity structure.
- *
- *
- * @sect3{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 (@p{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);
  *
- *   somesolver.solve (A, x, f,
- *                     PreconditionUseMatrix<SparseILU<double>,Vector<double> >
- *                           (ilu,&SparseILU<double>::template apply_decomposition<double>));
- * @end{verbatim}
+ * 
+ * @sect2{Usage and state management}
  *
+ * Refer to @ref{SparseLUDecomposition} documentation for suggested
+ * usage and state management.
  *
+ * 
  * @sect2{On template instantiations}
  *
  * Member functions of this class are either implemented in this file
  * @author Wolfgang Bangerth, 1999, based on a similar implementation by Malte Braack
  */
 template <typename number>
-class SparseILU : protected SparseMatrix<number>
+class SparseILU : public SparseLUDecomposition<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
-                                     * @p{reinit(SparsityPattern)}.
-                                     */
+                                     /**
+                                      * Constructor. Does nothing, so
+                                      * you have to call @p{reinit}
+                                      * sometimes afterwards.
+                                      */
     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 @p{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 @p{reinit} is not called with a
-                                     * new sparsity structure.
-                                     */
+                                     /**
+                                      * Constructor. Initialize the
+                                      * sparsity pattern of this
+                                      * object with the given
+                                      * argument.
+                                      */
     SparseILU (const SparsityPattern &sparsity);
 
-                                    /**
-                                     * Reinitialize the object but keep to
-                                     * the sparsity pattern previously used.
-                                     * This may be necessary if you @p{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 SparsityPattern &sparsity);
-
                                     /**
                                      * Perform the incomplete LU
                                      * factorization of the given
@@ -179,24 +95,42 @@ class SparseILU : protected SparseMatrix<number>
                                      * used by the matrix to be
                                      * decomposed.  Fill-in is thus
                                      * allowed.
+                                     *
+                                     * If @p{strengthen_diagonal}
+                                     * parameter is greater than
+                                     * zero, this method invokes
+                                     * @p{get_strengthen_diagonal_impl
+                                     * ()}.
+                                     *
+                                     * Refer to
+                                     * @ref{SparseLUDecomposition}
+                                     * documentation for state
+                                     * management.
                                      */
     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$.
+                                     * Same as @p{vmult}. This method
+                                     * is deprecated, and left for
+                                     * backward compability. It may
+                                     * be removed in later versions.
+                                     *
                                      */
     template <typename somenumber>
     void apply_decomposition (Vector<somenumber>       &dst,
                              const Vector<somenumber> &src) const;
 
                                     /**
-                                     * Same as
-                                     * @p{apply_decomposition},
-                                     * format for LAC.
+                                     * Apply the incomplete decomposition,
+                                     * i.e. do one forward-backward step
+                                     * $dst=(LU)^{-1}src$.
+                                     *
+                                     * Refer to
+                                     * @ref{SparseLUDecomposition}
+                                     * documentation for state
+                                     * management.
                                      */
     template <typename somenumber>
     void vmult (Vector<somenumber>       &dst,
@@ -229,14 +163,16 @@ class SparseILU : protected SparseMatrix<number>
                    << " is not greater or equal than zero!");
 };
 
+
+
 template <typename number>
 template <typename somenumber>
+inline
 void
-SparseILU<number>::vmult (Vector<somenumber>       &dst,
-                         const Vector<somenumber> &src) const
+SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
+                                        const Vector<somenumber> &src) const
 {
-  apply_decomposition(dst, src);
-}
-
+  vmult (dst, src);
+};
 
-#endif
+#endif // __deal2__sparse_ilu_h
diff --git a/deal.II/lac/include/lac/sparse_ilu.h.x b/deal.II/lac/include/lac/sparse_ilu.h.x
new file mode 100644 (file)
index 0000000..00359dc
--- /dev/null
@@ -0,0 +1,242 @@
+//----------------------------  sparse_ilu.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_ilu.h  ---------------------------
+#ifndef __deal2__sparse_ilu_h
+#define __deal2__sparse_ilu_h
+
+
+#include <base/config.h>
+#include <lac/sparse_matrix.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 @p{0}
+ * to @p{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 @p{SparseMatrix}. Since it is not a matrix in
+ * the usual sense, the derivation is @p{protected} rather than @p{public}.
+ *
+ * Note that in the algorithm given, the lower left part of the matrix base
+ * class is used to store the @p{L} part of the decomposition, while 
+ * the upper right part is used to store @p{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).
+ *
+ *
+ * @sect3{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 @p{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 @p{SparsityPattern} class which allows the addition
+ * of side-diagonals to a given sparsity structure.
+ *
+ *
+ * @sect3{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 (@p{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);
+ *
+ *   somesolver.solve (A, x, f,
+ *                     PreconditionUseMatrix<SparseILU<double>,Vector<double> >
+ *                           (ilu,&SparseILU<double>::template apply_decomposition<double>));
+ * @end{verbatim}
+ *
+ *
+ * @sect2{On template instantiations}
+ *
+ * Member functions of this class are either implemented in this file
+ * or in a file of the same name with suffix ``.templates.h''. For the
+ * most common combinations of the template parameters, instantiations
+ * of this class are provided in a file with suffix ``.cc'' in the
+ * ``source'' directory. If you need an instantiation that is not
+ * listed there, you have to include this file along with the
+ * corresponding ``.templates.h'' file and instantiate the respective
+ * class yourself.
+ *
+ * @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
+                                     * @p{reinit(SparsityPattern)}.
+                                     */
+    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 @p{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 @p{reinit} is not called with a
+                                     * new sparsity structure.
+                                     */
+    SparseILU (const SparsityPattern &sparsity);
+
+                                    /**
+                                     * Reinitialize the object but keep to
+                                     * the sparsity pattern previously used.
+                                     * This may be necessary if you @p{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 SparsityPattern &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;
+
+                                    /**
+                                     * Same as
+                                     * @p{apply_decomposition},
+                                     * format for LAC.
+                                     */
+    template <typename somenumber>
+    void vmult (Vector<somenumber>       &dst,
+               const Vector<somenumber> &src) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     */
+    unsigned int memory_consumption () 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!");
+};
+
+template <typename number>
+template <typename somenumber>
+void
+SparseILU<number>::vmult (Vector<somenumber>       &dst,
+                         const Vector<somenumber> &src) const
+{
+  apply_decomposition(dst, src);
+}
+
+
+#endif
index 608d889a0ccfadbce011650ea7def10f3e20703a..617da1df7f426c163937adbbd640a233cccecb45 100644 (file)
@@ -1,8 +1,6 @@
 //----------------------------  sparse_ilu.templates.h  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -10,8 +8,8 @@
 //    further information on this license.
 //
 //----------------------------  sparse_ilu.templates.h  ---------------------------
-#ifndef __deal2__sparse_ilu_templates_h
-#define __deal2__sparse_ilu_templates_h
+#ifndef _sparse_ilu_templates_h
+#define _sparse_ilu_templates_h
 
 
 
 
 
 template <typename number>
-SparseILU<number>::SparseILU () 
+SparseILU<number>::SparseILU ()
 {};
 
 
 
 template <typename number>
 SparseILU<number>::SparseILU (const SparsityPattern &sparsity) :
-               SparseMatrix<number> (sparsity)
+               SparseLUDecomposition<number> (sparsity)
 {};
 
 
 
-template <typename number>
-void SparseILU<number>::reinit ()
-{
-  SparseMatrix<number>::reinit ();
-};
-
-
-
-template <typename number>
-void SparseILU<number>::reinit (const SparsityPattern &sparsity)
-{
-  SparseMatrix<number>::reinit (sparsity);
-};
-
-
-
 template <typename number>
 template <typename somenumber>
 void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                   const double                    strengthen_diagonal)
 {
+  SparseLUDecomposition<number>::decompose (matrix, strengthen_diagonal);
   Assert (matrix.m()==matrix.n(), ExcMatrixNotSquare ());
   Assert (this->m()==this->n(),   ExcMatrixNotSquare ());
   Assert (matrix.m()==this->m(),  ExcSizeMismatch(matrix.m(), this->m()));
   
   Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal));
 
+  copy_from (matrix);
 
-                                  // first thing: copy over all elements
-                                  // of @p{matrix} to the present object
-                                  //
-                                  // note that some elements in this
-                                  // matrix may not be in @p{matrix},
-                                  // so we need to preset our matrix
-                                  // by zeroes.
-  if (true)
-    {
-                                      // preset the elements
-      std::fill_n (&this->global_entry(0),
-                  this->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 unsigned int * const column_numbers
-       = matrix.get_sparsity_pattern().get_column_numbers();
-      
-      for (unsigned int row=0; row<this->m(); ++row)
-       for (const unsigned 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)
+    strengthen_diagonal_impl();
 
-  if (strengthen_diagonal > 0)
-    for (unsigned int row=0; row<this->m(); ++row)
-      {
-                                        // get the length of the row
-                                        // (without the diagonal element)
-       const unsigned int
-         rowlength = (this->get_sparsity_pattern().get_rowstart_indices()[row+1]
-                      -
-                      this->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
-         = this->get_sparsity_pattern().get_rowstart_indices()[row] + 1;
-       number * const diagonal_element = &this->global_entry(rowstart-1);
-
-       number rowsum = 0;
-       for (unsigned int global_index=rowstart;
-            global_index<rowstart+rowlength; ++global_index)
-         rowsum += std::fabs(this->global_entry(global_index));
-
-       *diagonal_element += strengthen_diagonal * rowsum;
-      };
-
-
-                                  // now work only on this
-                                  // matrix
   const SparsityPattern             &sparsity = this->get_sparsity_pattern();
   const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices();
   const unsigned int * const column_numbers   = sparsity.get_column_numbers();
@@ -129,14 +60,14 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
   (indices=0..N-1)
   
   for i=1..N-1
-    a[i-1,i-1] = a[i-1,i-1]^{-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 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]
+  for j=k+1..N-1
+  if (a[i,j] exists & a[k,j] exists)
+  a[i,j] -= a[i,k] * a[k,j]
 */
 
 
@@ -162,13 +93,11 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                       // diagonal element at start
       const unsigned int * first_of_row
        = &column_numbers[rowstart_indices[row]+1];
-      const unsigned int * first_after_diagonal
-       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
-                           &column_numbers[rowstart_indices[row+1]],
-                           row);
+      const unsigned int * first_after_diagonal = prebuilt_lower_bound[row];
 
                                       // k := *col_ptr
-      for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++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);
@@ -192,7 +121,7 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
            {
 //TODO:[WB] make code faster by using the following comment          
                                               // note: this inner loop could
-                                              // be made considerable faster
+                                              // be made considerably faster
                                               // if we consulted the row
                                               // with number *col_ptr,
                                               // instead of always asking
@@ -215,61 +144,18 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                   // 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);
-
-/*
-  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);
-         };
-      };
-*/      
+  this->diag_element(this->m()-1) = 1./this->diag_element(this->m()-1);
 };
 
 
 
 template <typename number>
 template <typename somenumber>
-void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
-                                            const Vector<somenumber> &src) const 
+void SparseILU<number>::vmult (Vector<somenumber>       &dst,
+                               const Vector<somenumber> &src) const 
 {
+  SparseLUDecomposition<number>::vmult (dst, src);
+
   Assert (dst.size() == src.size(), ExcSizeMismatch(dst.size(), src.size()));
   Assert (dst.size() == this->m(), ExcSizeMismatch(dst.size(), this->m()));
   
@@ -298,10 +184,7 @@ void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
       const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
                                       // find the position where the part
                                       // right of the diagonal starts
-      const unsigned int * const first_after_diagonal
-       = std::lower_bound (rowstart,
-                           &column_numbers[rowstart_indices[row+1]],
-                           row);
+      const unsigned int * const first_after_diagonal = prebuilt_lower_bound[row];
       
       for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
        dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
@@ -321,10 +204,7 @@ void SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
       const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
                                       // find the position where the part
                                       // right of the diagonal starts
-      const unsigned int * const first_after_diagonal
-       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
-                           &column_numbers[rowstart_indices[row+1]],
-                           static_cast<unsigned int>(row));
+      const unsigned int * const first_after_diagonal = prebuilt_lower_bound[row];
       
       for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
        dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
@@ -342,7 +222,7 @@ template <typename number>
 unsigned int
 SparseILU<number>::memory_consumption () const
 {
-  return SparseMatrix<number>::memory_consumption ();
+  return SparseLUDecomposition<number>::memory_consumption ();
 };
 
 
diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h.x b/deal.II/lac/include/lac/sparse_ilu.templates.h.x
new file mode 100644 (file)
index 0000000..608d889
--- /dev/null
@@ -0,0 +1,353 @@
+//----------------------------  sparse_ilu.templates.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_ilu.templates.h  ---------------------------
+#ifndef __deal2__sparse_ilu_templates_h
+#define __deal2__sparse_ilu_templates_h
+
+
+
+#include <base/config.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 SparsityPattern &sparsity) :
+               SparseMatrix<number> (sparsity)
+{};
+
+
+
+template <typename number>
+void SparseILU<number>::reinit ()
+{
+  SparseMatrix<number>::reinit ();
+};
+
+
+
+template <typename number>
+void SparseILU<number>::reinit (const SparsityPattern &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 (this->m()==this->n(),   ExcMatrixNotSquare ());
+  Assert (matrix.m()==this->m(),  ExcSizeMismatch(matrix.m(), this->m()));
+  
+  Assert (strengthen_diagonal>=0, ExcInvalidStrengthening (strengthen_diagonal));
+
+
+                                  // first thing: copy over all elements
+                                  // of @p{matrix} to the present object
+                                  //
+                                  // note that some elements in this
+                                  // matrix may not be in @p{matrix},
+                                  // so we need to preset our matrix
+                                  // by zeroes.
+  if (true)
+    {
+                                      // preset the elements
+      std::fill_n (&this->global_entry(0),
+                  this->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 unsigned int * const column_numbers
+       = matrix.get_sparsity_pattern().get_column_numbers();
+      
+      for (unsigned int row=0; row<this->m(); ++row)
+       for (const unsigned 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<this->m(); ++row)
+      {
+                                        // get the length of the row
+                                        // (without the diagonal element)
+       const unsigned int
+         rowlength = (this->get_sparsity_pattern().get_rowstart_indices()[row+1]
+                      -
+                      this->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
+         = this->get_sparsity_pattern().get_rowstart_indices()[row] + 1;
+       number * const diagonal_element = &this->global_entry(rowstart-1);
+
+       number rowsum = 0;
+       for (unsigned int global_index=rowstart;
+            global_index<rowstart+rowlength; ++global_index)
+         rowsum += std::fabs(this->global_entry(global_index));
+
+       *diagonal_element += strengthen_diagonal * rowsum;
+      };
+
+
+                                  // now work only on this
+                                  // matrix
+  const SparsityPattern             &sparsity = this->get_sparsity_pattern();
+  const unsigned int * 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]
+*/
+
+
+                                  // i := row
+  for (unsigned int row=1; row<this->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
+      AssertThrow((this->global_entry(rowstart_indices[row-1]) !=0),
+                 ExcDivideByZero());
+      
+      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
+       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
+                           &column_numbers[rowstart_indices[row+1]],
+                           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
+         const int global_index_ki = sparsity(*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)
+           {
+//TODO:[WB] make code faster by using the following comment          
+                                              // 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.
+                     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))
+               this->global_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);
+
+/*
+  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() == this->m(), ExcSizeMismatch(dst.size(), this->m()));
+  
+  const unsigned int N=dst.size();
+  const unsigned int * const rowstart_indices
+    = this->get_sparsity_pattern().get_rowstart_indices();
+  const unsigned int * const column_numbers
+    = this->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 unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
+                                      // find the position where the part
+                                      // right of the diagonal starts
+      const unsigned int * const first_after_diagonal
+       = std::lower_bound (rowstart,
+                           &column_numbers[rowstart_indices[row+1]],
+                           row);
+      
+      for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col)
+       dst(row) -= this->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 unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
+                                      // find the position where the part
+                                      // right of the diagonal starts
+      const unsigned int * const first_after_diagonal
+       = std::lower_bound (&column_numbers[rowstart_indices[row]+1],
+                           &column_numbers[rowstart_indices[row+1]],
+                           static_cast<unsigned int>(row));
+      
+      for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
+       dst(row) -= this->global_entry (col-column_numbers) * dst(*col);
+
+                                      // scale by the diagonal element.
+                                      // note that the diagonal element
+                                      // was stored inverted
+      dst(row) *= this->diag_element(row);
+    };
+};
+
+
+
+template <typename number>
+unsigned int
+SparseILU<number>::memory_consumption () const
+{
+  return SparseMatrix<number>::memory_consumption ();
+};
+
+
+
+/*----------------------------   sparse_ilu.templates.h     ---------------------------*/
+
+#endif
+/*----------------------------   sparse_ilu.templates.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/sparse_mic.h b/deal.II/lac/include/lac/sparse_mic.h
new file mode 100644 (file)
index 0000000..edcdb87
--- /dev/null
@@ -0,0 +1,201 @@
+//----------------------------  sparse_mic.h  ---------------------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_mic.h  ---------------------------
+#ifndef __deal2__sparse_mic_h
+#define __deal2__sparse_mic_h
+
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_decomposition.h>
+
+
+/**
+ * Modified incomplete Cholesky (MIC(0)) preconditioner.  This class
+ * conforms to the state and usage specification in
+ * @ref{SparseLUDecomposition}.
+ *
+ * 
+ * @sect2{The decomposition}
+ * 
+ * Let a sparse matrix A is in the form A = - L - U + D, where -L and
+ * -U are strictly lower and upper triangular matrices. The MIC(0)
+ * decomposition of the matrix A is defined by B = (X-L)X^(-1)(X-U),
+ * where X is a diagonal matrix, defined by the condition rowsum(A) =
+ * rowsum(B).
+ * 
+ * @author Stephen "Cheffo" Kolaroff
+ */
+template <typename number>
+class SparseMIC : public SparseLUDecomposition<number>
+{
+  public:
+                                     /**
+                                      * Constructor. Does nothing, so
+                                      * you have to call @p{reinit}
+                                      * sometimes afterwards.
+                                      */
+    SparseMIC ();
+
+                                     /**
+                                      * Constructor. Initialize the
+                                      * sparsity pattern of this
+                                      * object with the given
+                                      * argument.
+                                      */
+    SparseMIC (const SparsityPattern &sparsity);
+
+
+                                    /**
+                                     * Reinitialize the object but
+                                     * keep to the sparsity pattern
+                                     * previously used.  This may be
+                                     * necessary if you @p{reinit}'d
+                                     * the sparsity structure and
+                                     * want to update the size of the
+                                     * matrix.
+                                     *
+                                     * After this method is invoked,
+                                     * this object is out of synch
+                                     * (not decomposed state).
+                                     *
+                                     * This function only releases
+                                     * some memory and calls the
+                                     * respective function of the
+                                     * base class.
+                                     */
+    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 only releases
+                                     * some memory and calls the
+                                     * respective function of the
+                                     * base class.
+                                     */
+    void reinit (const SparsityPattern &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.
+                                     *
+                                     * If @p{strengthen_diagonal}
+                                     * parameter is greater than
+                                     * zero, this method invokes
+                                     * @p{get_strengthen_diagonal_impl
+                                     * ()}.
+                                     *
+                                     * Refer to
+                                     * @ref{SparseLUDecomposition}
+                                     * documentation for state
+                                     * management.
+                                     */
+    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$.
+                                     *
+                                     * Refer to
+                                     * @ref{SparseLUDecomposition}
+                                     * documentation for state
+                                     * management.
+                                     */
+    template <typename somenumber>
+    void vmult (Vector<somenumber>       &dst,
+                const Vector<somenumber> &src) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     */
+    unsigned int memory_consumption () const;
+
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcMatrixNotSquare);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcInternal);
+                                     /**
+                                      * 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!");
+                                     /**
+                                      * Exception
+                                      */
+    DeclException2(ExcDecompositionNotStable, int, double,
+                  << "The diagonal element (" <<arg1<<","<<arg1<<") is "
+                  << arg2 <<", but must be positive");
+    
+  private:
+                                     /**
+                                      * Values of the computed
+                                      * diagonal.
+                                      */
+    std::vector<number> diag;
+    
+                                     /**
+                                      * Inverses of the the diagonal:
+                                      * precomputed for faster vmult.
+                                      */
+    std::vector<number> inv_diag;
+
+                                     /**
+                                      * Values of the computed "inner
+                                      * sums", i.e. per-row sums of
+                                      * the elements laying on the
+                                      * right side of the diagonal.
+                                      */
+    std::vector<number> inner_sums;
+    
+                                     /**
+                                      * Compute the row-th "inner
+                                      * sum".
+                                      */
+    number get_rowsum (const unsigned int row) const;
+};
+
+
+
+#endif  // __deal2__
diff --git a/deal.II/lac/include/lac/sparse_mic.templates.h b/deal.II/lac/include/lac/sparse_mic.templates.h
new file mode 100644 (file)
index 0000000..b8fe533
--- /dev/null
@@ -0,0 +1,234 @@
+//----------------------------  sparse_mic.templates.h  ---------------------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_mic.templates.h  ---------------------------
+#ifndef __deal2__sparse_mic_templates_h
+#define __deal2__sparse_mic_templates_h
+
+
+#include <base/memory_consumption.h>
+#include <lac/sparse_mic.h>
+#include <lac/vector.h>
+
+
+template <typename number>
+SparseMIC<number>::SparseMIC ()
+                :
+                diag(0),
+                inv_diag(0),
+                inner_sums(0)
+{};
+
+
+
+template <typename number>
+SparseMIC<number>::SparseMIC (const SparsityPattern &sparsity)
+                :
+                SparseLUDecomposition<number> (sparsity),
+                diag(0),
+                inv_diag(0),
+                inner_sums(0)
+{};
+
+
+
+template <typename number>
+void
+SparseMIC<number>::reinit ()
+{
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (diag);
+    };
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (inv_diag);
+    };
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (inner_sums);
+    };
+  
+  SparseLUDecomposition<number>::reinit ();
+}
+
+
+
+template <typename number>
+void SparseMIC<number>::reinit (const SparsityPattern& sparsity)
+{
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (diag);
+    };
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (inv_diag);
+    };
+  if (true)
+    {
+      std::vector<number> tmp;
+      tmp.swap (inner_sums);
+    };
+  SparseLUDecomposition<number>::reinit (sparsity);
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void SparseMIC<number>::decompose (const SparseMatrix<somenumber> &matrix,
+                                  const double                    strengthen_diagonal)
+{
+
+  SparseLUDecomposition<number>::decompose(matrix, 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));
+
+  if (strengthen_diagonal > 0)
+    strengthen_diagonal_impl ();
+
+                                   // MIC implementation: (S. Margenov lectures)
+                                   // x[i] = a[i][i] - sum(k=1, i-1,
+                                   //              a[i][k]/x[k]*sum(j=k+1, N, a[k][j]))
+       
+                                   // TODO: for sake of siplicity,
+                                   // those are placed here A better
+                                   // implementation would store this
+                                   // values in the underlying sparse
+                                   // matrix itself.
+  diag.resize (m());
+  inv_diag.resize (m());
+  inner_sums.resize (m());
+
+                                   // precalc sum(j=k+1, N, a[k][j]))
+  for(unsigned int row=0; row<m(); row++) {
+    inner_sums[row] = get_rowsum(row);
+  }
+
+  const unsigned int* const col_nums = get_sparsity_pattern().get_column_numbers();
+  const unsigned int* const rowstarts = get_sparsity_pattern().get_rowstart_indices();
+
+  for(unsigned int row=0; row<m(); row++) {
+    number temp = diag_element(row);
+    number temp1 = 0;
+    const unsigned int * const first_after_diagonal = prebuilt_lower_bound[row];
+
+    unsigned int k = 0;
+    for (const unsigned int * col=&col_nums[rowstarts[row]+1];
+         col<first_after_diagonal; ++col, k++)
+      temp1 += matrix.global_entry (col-col_nums)/diag[k]*inner_sums[k];
+
+    diag[row] = temp - temp1;
+    inv_diag[row] = 1.0/diag[row];
+    Assert(diag[row]>0, ExcInternal());
+  }
+};
+
+
+
+template <typename number>
+inline number
+SparseMIC<number>::get_rowsum (const unsigned int row) const
+{
+  Assert(m()==n(), ExcMatrixNotSquare());
+                                   // get start of this row. skip the
+                                   // diagonal element
+  const unsigned int * const column_numbers = get_sparsity_pattern().get_column_numbers();
+  const unsigned int * const rowstart_indices = get_sparsity_pattern().get_rowstart_indices();
+  const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
+
+                                   // find the position where the part
+                                   // right of the diagonal starts
+  const unsigned int * const first_after_diagonal = prebuilt_lower_bound[row];
+  number rowsum =  0;
+  for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col)
+    rowsum += global_entry (col-column_numbers);
+
+  return rowsum;       
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMIC<number>::vmult (Vector<somenumber>       &dst,
+                          const Vector<somenumber> &src) const
+{
+  SparseLUDecomposition<number>::vmult (dst, src);
+  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 unsigned int * const column_numbers   = get_sparsity_pattern().get_column_numbers();
+                                   // We assume the underlying matrix A is:
+                                   // A = X - L - U, where -L and -U are
+                                   // strictly lower- and upper- diagonal
+                                   // parts of the system.
+                                   // 
+                                   // Solve (X-L)X{-1}(X-U) x = b
+                                   // in 3 steps:
+  dst = src;
+  for (unsigned int row=0; row<N; ++row)
+    {
+                                       // Now: (X-L)u = b
+
+                                       // get start of this row. skip
+                                       // the diagonal element
+      const unsigned int * const rowstart = &column_numbers[rowstart_indices[row]+1];
+      const unsigned int * const fad = prebuilt_lower_bound[row];
+      for (const unsigned int * col=rowstart; col!=fad; ++col)
+        dst(row) -= global_entry (col-column_numbers) * dst(*col);
+      
+      dst(row) *= inv_diag[row];
+    };
+
+                                   // Now: v = Xu
+  for(unsigned int row=0; row<N; row++)
+    dst(row) *= diag[row];
+
+                                   // x = (X-U)v
+  for (int row=N-1; row>=0; --row)
+    {
+                                      // get end of this row
+      const unsigned int * const rowend = &column_numbers[rowstart_indices[row+1]];
+      const  unsigned int * const fad = prebuilt_lower_bound[row];
+      for (const unsigned int * col=fad; col!=rowend; ++col)
+        dst(row) -= global_entry (col-column_numbers) * dst(*col);
+
+      dst(row) *= inv_diag[row];
+    };
+};
+
+
+
+template <typename number>
+unsigned int
+SparseMIC<number>::memory_consumption () const
+{
+  return (SparseLUDecomposition<number>::memory_consumption () +
+          MemoryConsumption::memory_consumption(diag) +
+          MemoryConsumption::memory_consumption(inv_diag) +
+          MemoryConsumption::memory_consumption(inner_sums));
+};
+
+
+
+#endif // __deal2__sparse_mic_templates_h
diff --git a/deal.II/lac/source/sparse_decomposition.cc b/deal.II/lac/source/sparse_decomposition.cc
new file mode 100644 (file)
index 0000000..49b4399
--- /dev/null
@@ -0,0 +1,25 @@
+//----------------------------  sparse_decomposition.cc  ---------------------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_decomposition.h  ---------------------------
+
+#include <lac/sparse_decomposition.templates.h>
+
+
+template class SparseLUDecomposition<double>;
+template void SparseLUDecomposition<double>::decompose<double> (const SparseMatrix<double> &,
+                                                                const double);
+template void SparseLUDecomposition<double>::decompose<float> (const SparseMatrix<float> &,
+                                                               const double);
+
+template class SparseLUDecomposition<float>;
+template void SparseLUDecomposition<float>::decompose<double> (const SparseMatrix<double> &,
+                                                               const double);
+template void SparseLUDecomposition<float>::decompose<float> (const SparseMatrix<float> &,
+                                                              const double);
index 771c729e57b25d8cf79b434e11284b1e8e285f4b..c29b5fd34f38a626d75bbb60f2ee9b16b945e1b8 100644 (file)
@@ -1,7 +1,4 @@
 //----------------------------  sparse_ilu.cc  ---------------------------
-//    $Id$
-//    Version: $Name$
-//
 //    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
@@ -12,6 +9,7 @@
 //----------------------------  sparse_ilu.cc  ---------------------------
 
 
+#include <lac/sparse_ilu.h>
 #include <lac/sparse_ilu.templates.h>
 
 
 template class SparseILU<double>;
 template void SparseILU<double>::decompose<double> (const SparseMatrix<double> &,
                                                    const double);
-template void SparseILU<double>::apply_decomposition<double> (Vector<double> &,
-                                                             const Vector<double> &) const;
+template void SparseILU<double>::vmult <double> (Vector<double> &,
+                                                 const Vector<double> &) const;
 template void SparseILU<double>::decompose<float> (const SparseMatrix<float> &,
                                                   const double);
-template void SparseILU<double>::apply_decomposition<float> (Vector<float> &,
-                                                            const Vector<float> &) const;
+template void SparseILU<double>::vmult<float> (Vector<float> &,
+                                               const Vector<float> &) const;
 
 
 template class SparseILU<float>;
 template void SparseILU<float>::decompose<double> (const SparseMatrix<double> &,
                                                   const double);
-template void SparseILU<float>::apply_decomposition<double> (Vector<double> &,
-                                                            const Vector<double> &) const;
+template void SparseILU<float>::vmult<double> (Vector<double> &,
+                                               const Vector<double> &) const;
 template void SparseILU<float>::decompose<float> (const SparseMatrix<float> &,
                                                  const double);
-template void SparseILU<float>::apply_decomposition<float> (Vector<float> &,
-                                                           const Vector<float> &) const;
+template void SparseILU<float>::vmult<float> (Vector<float> &,
+                                              const Vector<float> &) const;
diff --git a/deal.II/lac/source/sparse_ilu.cc.x b/deal.II/lac/source/sparse_ilu.cc.x
new file mode 100644 (file)
index 0000000..771c729
--- /dev/null
@@ -0,0 +1,38 @@
+//----------------------------  sparse_ilu.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_ilu.cc  ---------------------------
+
+
+#include <lac/sparse_ilu.templates.h>
+
+
+// explicit instantiations
+template class SparseILU<double>;
+template void SparseILU<double>::decompose<double> (const SparseMatrix<double> &,
+                                                   const double);
+template void SparseILU<double>::apply_decomposition<double> (Vector<double> &,
+                                                             const Vector<double> &) const;
+template void SparseILU<double>::decompose<float> (const SparseMatrix<float> &,
+                                                  const double);
+template void SparseILU<double>::apply_decomposition<float> (Vector<float> &,
+                                                            const Vector<float> &) const;
+
+
+template class SparseILU<float>;
+template void SparseILU<float>::decompose<double> (const SparseMatrix<double> &,
+                                                  const double);
+template void SparseILU<float>::apply_decomposition<double> (Vector<double> &,
+                                                            const Vector<double> &) const;
+template void SparseILU<float>::decompose<float> (const SparseMatrix<float> &,
+                                                 const double);
+template void SparseILU<float>::apply_decomposition<float> (Vector<float> &,
+                                                           const Vector<float> &) const;
diff --git a/deal.II/lac/source/sparse_mic.cc b/deal.II/lac/source/sparse_mic.cc
new file mode 100644 (file)
index 0000000..35b2a41
--- /dev/null
@@ -0,0 +1,35 @@
+//----------------------------  sparse_mic.cc  ---------------------------
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002
+//    by the deal.II authors and Stephen "Cheffo" Kolaroff
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparse_mic.cc  ---------------------------
+
+#include <lac/sparse_mic.templates.h>
+
+
+// explicit instantiations for double and float matrices
+template class SparseMIC<double>;
+template void SparseMIC<double>::decompose<double> (const SparseMatrix<double> &,
+                                                    const double);
+template void SparseMIC<double>::vmult<double> (Vector<double> &,
+                                                const Vector<double> &) const;
+template void SparseMIC<double>::decompose<float> (const SparseMatrix<float> &,
+                                                   const double);
+template void SparseMIC<double>::vmult<float> (Vector<float> &,
+                                               const Vector<float> &) const;
+
+template class SparseMIC<float>;
+template void SparseMIC<float>::decompose<double> (const SparseMatrix<double> &,
+                                                   const double);
+template void SparseMIC<float>::vmult<double> (Vector<double> &,
+                                               const Vector<double> &) const;
+template void SparseMIC<float>::decompose<float> (const SparseMatrix<float> &,
+                                                  const double);
+template void SparseMIC<float>::vmult<float> (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.