From e97c66a4e6846db1386b205907de489944387255 Mon Sep 17 00:00:00 2001 From: wolf Date: Mon, 2 Dec 2002 17:28:05 +0000 Subject: [PATCH] Add Stephens implementation of sparse decompositions. git-svn-id: https://svn.dealii.org/trunk@6792 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/sparse_decomposition.h | 327 ++++++++++++++++ .../lac/sparse_decomposition.templates.h | 170 +++++++++ deal.II/lac/include/lac/sparse_ilu.h | 170 +++------ deal.II/lac/include/lac/sparse_ilu.h.x | 242 ++++++++++++ .../lac/include/lac/sparse_ilu.templates.h | 176 ++------- .../lac/include/lac/sparse_ilu.templates.h.x | 353 ++++++++++++++++++ deal.II/lac/include/lac/sparse_mic.h | 201 ++++++++++ .../lac/include/lac/sparse_mic.templates.h | 234 ++++++++++++ deal.II/lac/source/sparse_decomposition.cc | 25 ++ deal.II/lac/source/sparse_ilu.cc | 20 +- deal.II/lac/source/sparse_ilu.cc.x | 38 ++ deal.II/lac/source/sparse_mic.cc | 35 ++ 12 files changed, 1715 insertions(+), 276 deletions(-) create mode 100644 deal.II/lac/include/lac/sparse_decomposition.h create mode 100644 deal.II/lac/include/lac/sparse_decomposition.templates.h create mode 100644 deal.II/lac/include/lac/sparse_ilu.h.x create mode 100644 deal.II/lac/include/lac/sparse_ilu.templates.h.x create mode 100644 deal.II/lac/include/lac/sparse_mic.h create mode 100644 deal.II/lac/include/lac/sparse_mic.templates.h create mode 100644 deal.II/lac/source/sparse_decomposition.cc create mode 100644 deal.II/lac/source/sparse_ilu.cc.x create mode 100644 deal.II/lac/source/sparse_mic.cc diff --git a/deal.II/lac/include/lac/sparse_decomposition.h b/deal.II/lac/include/lac/sparse_decomposition.h new file mode 100644 index 0000000000..24b13647e4 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_decomposition.h @@ -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 +#include + +#include + +/** + * 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 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 +class SparseLUDecomposition : protected SparseMatrix{ + 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 + void decompose (const SparseMatrix &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 + void copy_from (const SparseMatrix& 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}.
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 prebuilt_lower_bound; + + private: + /** + * Fills the + * @ref{prebuilt_lower_bound} + * array. + */ + void prebuild_lower_bound (); + +}; + + + +template +inline number +SparseLUDecomposition:: +get_strengthen_diagonal(const number rowsum, const unsigned int row) const +{ + return strengthen_diagonal; +}; + + + +template +inline bool +SparseLUDecomposition::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 index 0000000000..52911d63ff --- /dev/null +++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h @@ -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 +#include +#include + + +template +SparseLUDecomposition::SparseLUDecomposition() + : + SparseMatrix(), + decomposed(false) +{} + + + +template +SparseLUDecomposition:: +SparseLUDecomposition (const SparsityPattern& sparsity) : + SparseMatrix(sparsity), + decomposed(false) +{} + + + +template +SparseLUDecomposition::~SparseLUDecomposition() +{} + + + +template +template +void +SparseLUDecomposition:: +decompose (const SparseMatrix &matrix, + const double strengthen_diagonal) +{ + decomposed = false; + + this->strengthen_diagonal = strengthen_diagonal; + prebuild_lower_bound (); + copy_from (matrix); + decomposed = true; +} + + + +template +void +SparseLUDecomposition::reinit () +{ + decomposed = false; + if (true) + { + std::vector tmp; + tmp.swap (prebuilt_lower_bound); + }; + SparseMatrix::reinit (); +} + + + +template +void SparseLUDecomposition::reinit (const SparsityPattern& sparsity) +{ + decomposed = false; + if (true) + { + std::vector tmp; + tmp.swap (prebuilt_lower_bound); + }; + SparseMatrix::reinit (sparsity); +} + + + +template +void +SparseLUDecomposition::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 +template +void +SparseLUDecomposition::copy_from (const SparseMatrix& 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 +void +SparseLUDecomposition::strengthen_diagonal_impl () +{ + for (unsigned int row=0; row +inline unsigned int +SparseLUDecomposition::memory_consumption () const +{ + unsigned int + res = (SparseMatrix::memory_consumption () + + MemoryConsumption::memory_consumption(prebuilt_lower_bound)); + return res; +} + + diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index 00359dcdc9..e4a3746114 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -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 @@ -10,12 +8,14 @@ // further information on this license. // //---------------------------- sparse_ilu.h --------------------------- + #ifndef __deal2__sparse_ilu_h #define __deal2__sparse_ilu_h #include #include +#include /** @@ -40,50 +40,14 @@ * 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 ilu (ilu_sparsity); - * ilu.decompose (global_matrix); * - * somesolver.solve (A, x, f, - * PreconditionUseMatrix,Vector > - * (ilu,&SparseILU::template apply_decomposition)); - * @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 @@ -98,72 +62,24 @@ * @author Wolfgang Bangerth, 1999, based on a similar implementation by Malte Braack */ template -class SparseILU : protected SparseMatrix +class SparseILU : public SparseLUDecomposition { 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 * 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 void decompose (const SparseMatrix &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 void apply_decomposition (Vector &dst, const Vector &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 void vmult (Vector &dst, @@ -229,14 +163,16 @@ class SparseILU : protected SparseMatrix << " is not greater or equal than zero!"); }; + + template template +inline void -SparseILU::vmult (Vector &dst, - const Vector &src) const +SparseILU::apply_decomposition (Vector &dst, + const Vector &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 index 0000000000..00359dcdc9 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_ilu.h.x @@ -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 +#include + + +/** + * 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 ilu (ilu_sparsity); + * ilu.decompose (global_matrix); + * + * somesolver.solve (A, x, f, + * PreconditionUseMatrix,Vector > + * (ilu,&SparseILU::template apply_decomposition)); + * @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 +class SparseILU : protected SparseMatrix +{ + 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 + void decompose (const SparseMatrix &matrix, + const double strengthen_diagonal=0.); + + /** + * Apply the incomplete decomposition, + * i.e. do one forward-backward step + * $dst=(LU)^{-1}src$. + */ + template + void apply_decomposition (Vector &dst, + const Vector &src) const; + + /** + * Same as + * @p{apply_decomposition}, + * format for LAC. + */ + template + void vmult (Vector &dst, + const Vector &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 +template +void +SparseILU::vmult (Vector &dst, + const Vector &src) const +{ + apply_decomposition(dst, src); +} + + +#endif diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 608d889a0c..617da1df7f 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -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 @@ -24,102 +22,35 @@ template -SparseILU::SparseILU () +SparseILU::SparseILU () {}; template SparseILU::SparseILU (const SparsityPattern &sparsity) : - SparseMatrix (sparsity) + SparseLUDecomposition (sparsity) {}; -template -void SparseILU::reinit () -{ - SparseMatrix::reinit (); -}; - - - -template -void SparseILU::reinit (const SparsityPattern &sparsity) -{ - SparseMatrix::reinit (sparsity); -}; - - - template template void SparseILU::decompose (const SparseMatrix &matrix, const double strengthen_diagonal) { + SparseLUDecomposition::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; rowm(); ++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; rowm(); ++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_indexglobal_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::decompose (const SparseMatrix &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::decompose (const SparseMatrix &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::decompose (const SparseMatrix &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::decompose (const SparseMatrix &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; kdiag_element(this->m()-1) = 1./this->diag_element(this->m()-1); }; template template -void SparseILU::apply_decomposition (Vector &dst, - const Vector &src) const +void SparseILU::vmult (Vector &dst, + const Vector &src) const { + SparseLUDecomposition::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::apply_decomposition (Vector &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::apply_decomposition (Vector &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(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 unsigned int SparseILU::memory_consumption () const { - return SparseMatrix::memory_consumption (); + return SparseLUDecomposition::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 index 0000000000..608d889a0c --- /dev/null +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h.x @@ -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 +#include +#include + +#include +#include + + +template +SparseILU::SparseILU () +{}; + + + +template +SparseILU::SparseILU (const SparsityPattern &sparsity) : + SparseMatrix (sparsity) +{}; + + + +template +void SparseILU::reinit () +{ + SparseMatrix::reinit (); +}; + + + +template +void SparseILU::reinit (const SparsityPattern &sparsity) +{ + SparseMatrix::reinit (sparsity); +}; + + + +template +template +void SparseILU::decompose (const SparseMatrix &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; rowm(); ++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; rowm(); ++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_indexglobal_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; rowm(); ++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 +template +void SparseILU::apply_decomposition (Vector &dst, + const Vector &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; rowglobal_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(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 +unsigned int +SparseILU::memory_consumption () const +{ + return SparseMatrix::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 index 0000000000..edcdb8701f --- /dev/null +++ b/deal.II/lac/include/lac/sparse_mic.h @@ -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 +#include + + +/** + * 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 +class SparseMIC : public SparseLUDecomposition +{ + 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 + void decompose (const SparseMatrix &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 + void vmult (Vector &dst, + const Vector &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 (" < diag; + + /** + * Inverses of the the diagonal: + * precomputed for faster vmult. + */ + std::vector 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 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 index 0000000000..b8fe533948 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_mic.templates.h @@ -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 +#include +#include + + +template +SparseMIC::SparseMIC () + : + diag(0), + inv_diag(0), + inner_sums(0) +{}; + + + +template +SparseMIC::SparseMIC (const SparsityPattern &sparsity) + : + SparseLUDecomposition (sparsity), + diag(0), + inv_diag(0), + inner_sums(0) +{}; + + + +template +void +SparseMIC::reinit () +{ + if (true) + { + std::vector tmp; + tmp.swap (diag); + }; + if (true) + { + std::vector tmp; + tmp.swap (inv_diag); + }; + if (true) + { + std::vector tmp; + tmp.swap (inner_sums); + }; + + SparseLUDecomposition::reinit (); +} + + + +template +void SparseMIC::reinit (const SparsityPattern& sparsity) +{ + if (true) + { + std::vector tmp; + tmp.swap (diag); + }; + if (true) + { + std::vector tmp; + tmp.swap (inv_diag); + }; + if (true) + { + std::vector tmp; + tmp.swap (inner_sums); + }; + SparseLUDecomposition::reinit (sparsity); +} + + + +template +template +void SparseMIC::decompose (const SparseMatrix &matrix, + const double strengthen_diagonal) +{ + + SparseLUDecomposition::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; row0, ExcInternal()); + } +}; + + + +template +inline number +SparseMIC::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 +template +void +SparseMIC::vmult (Vector &dst, + const Vector &src) const +{ + SparseLUDecomposition::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=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 +unsigned int +SparseMIC::memory_consumption () const +{ + return (SparseLUDecomposition::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 index 0000000000..49b43992c1 --- /dev/null +++ b/deal.II/lac/source/sparse_decomposition.cc @@ -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 + + +template class SparseLUDecomposition; +template void SparseLUDecomposition::decompose (const SparseMatrix &, + const double); +template void SparseLUDecomposition::decompose (const SparseMatrix &, + const double); + +template class SparseLUDecomposition; +template void SparseLUDecomposition::decompose (const SparseMatrix &, + const double); +template void SparseLUDecomposition::decompose (const SparseMatrix &, + const double); diff --git a/deal.II/lac/source/sparse_ilu.cc b/deal.II/lac/source/sparse_ilu.cc index 771c729e57..c29b5fd34f 100644 --- a/deal.II/lac/source/sparse_ilu.cc +++ b/deal.II/lac/source/sparse_ilu.cc @@ -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 #include @@ -19,20 +17,20 @@ template class SparseILU; template void SparseILU::decompose (const SparseMatrix &, const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; +template void SparseILU::vmult (Vector &, + const Vector &) const; template void SparseILU::decompose (const SparseMatrix &, const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; +template void SparseILU::vmult (Vector &, + const Vector &) const; template class SparseILU; template void SparseILU::decompose (const SparseMatrix &, const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; +template void SparseILU::vmult (Vector &, + const Vector &) const; template void SparseILU::decompose (const SparseMatrix &, const double); -template void SparseILU::apply_decomposition (Vector &, - const Vector &) const; +template void SparseILU::vmult (Vector &, + const Vector &) 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 index 0000000000..771c729e57 --- /dev/null +++ b/deal.II/lac/source/sparse_ilu.cc.x @@ -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 + + +// explicit instantiations +template class SparseILU; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; + + +template class SparseILU; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; +template void SparseILU::decompose (const SparseMatrix &, + const double); +template void SparseILU::apply_decomposition (Vector &, + const Vector &) const; diff --git a/deal.II/lac/source/sparse_mic.cc b/deal.II/lac/source/sparse_mic.cc new file mode 100644 index 0000000000..35b2a41afb --- /dev/null +++ b/deal.II/lac/source/sparse_mic.cc @@ -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 + + +// explicit instantiations for double and float matrices +template class SparseMIC; +template void SparseMIC::decompose (const SparseMatrix &, + const double); +template void SparseMIC::vmult (Vector &, + const Vector &) const; +template void SparseMIC::decompose (const SparseMatrix &, + const double); +template void SparseMIC::vmult (Vector &, + const Vector &) const; + +template class SparseMIC; +template void SparseMIC::decompose (const SparseMatrix &, + const double); +template void SparseMIC::vmult (Vector &, + const Vector &) const; +template void SparseMIC::decompose (const SparseMatrix &, + const double); +template void SparseMIC::vmult (Vector &, + const Vector &) const; + -- 2.39.5