--- /dev/null
+//---------------------- 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
--- /dev/null
+//--------------------- 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;
+}
+
+
//---------------------------- 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
* 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,
<< " 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
--- /dev/null
+//---------------------------- 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
//---------------------------- 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
// 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();
(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]
*/
// 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);
{
//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
// 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()));
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);
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);
unsigned int
SparseILU<number>::memory_consumption () const
{
- return SparseMatrix<number>::memory_consumption ();
+ return SparseLUDecomposition<number>::memory_consumption ();
};
--- /dev/null
+//---------------------------- 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 ---------------------------*/
--- /dev/null
+//---------------------------- 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__
--- /dev/null
+//---------------------------- 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
--- /dev/null
+//---------------------------- 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);
//---------------------------- 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
//---------------------------- 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;
--- /dev/null
+//---------------------------- 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;
--- /dev/null
+//---------------------------- 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;
+