From 0314fb18a90a0fd6055d3204284a53772f55b60e Mon Sep 17 00:00:00 2001 From: hartmann Date: Wed, 19 Feb 2003 17:39:03 +0000 Subject: [PATCH] Add initialize and clear functions in order to unify the interface of preconditioners. Reinit, decompose, apply_decompose functions are now deprecated. Large parts of doc changed. git-svn-id: https://svn.dealii.org/trunk@7183 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/sparse_decomposition.h | 320 ++++++++++++------ .../lac/sparse_decomposition.templates.h | 84 ++++- deal.II/lac/include/lac/sparse_ilu.h | 94 +++-- .../lac/include/lac/sparse_ilu.templates.h | 23 +- deal.II/lac/include/lac/sparse_mic.h | 100 ++---- .../lac/include/lac/sparse_mic.templates.h | 42 +++ 6 files changed, 432 insertions(+), 231 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_decomposition.h b/deal.II/lac/include/lac/sparse_decomposition.h index 04f09bf192..8877120d9e 100644 --- a/deal.II/lac/include/lac/sparse_decomposition.h +++ b/deal.II/lac/include/lac/sparse_decomposition.h @@ -20,40 +20,67 @@ * 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 + * The decomposition is stored as a sparse matrix 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. + * 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{initialize} + * 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{Unified use of preconditioners} * - * @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): + * An object of this class can be used in the same form as all + * @ref{PreconditionBlock} preconditioners: * @begin{verbatim} - * SparseLUImplementation lu (lu_sparsity); - * lu.decompose (global_matrix); + * SparseLUImplementation lu; + * lu.initialize(matrix, SparseLUImplementation::AdditionalData(...)); * - * somesolver.solve (A, x, f, lu); + * somesolver.solve (A, x, f, lu); * @end{verbatim} * + * Through the @p{AdditionalData} object it is possible to specify + * additional parameters of the LU decomposition. + * + * 1/ The matrix diagonals can be strengthened by adding + * @p{strengthen_diagonal} times the sum of the absolute row entries + * of each row to the respective diagonal entries. By default no + * strengthening is performed. + * + * 2/ By default, each @p{initialize} function call creates its own + * sparsity. For that, it copies the sparsity of @p{matrix} and adds a + * specific number of extra off diagonal entries specified by + * @p{extra_off_diagonals}. + * + * 3/ By setting @p{use_previous_sparsity=true} the sparsity is not + * recreated but the sparsity of the previous @p{initialize} call is + * reused (recycled). This might be useful when several linear + * problems on the same sparsity need to solved, as for example + * several Newton iteration steps on the same triangulation. The + * default is @p{false}. + * + * 4/ It is possible to give a user defined sparsity to + * @p{use_this_sparsity}. Then, no sparsity is created but + * @p[*use_this_sparsity} is used to store the decomposed matrix. For + * restrictions on the sparsity see section `Fill-in' above). + * + * * @sect2{State management} * + * The state management simply requires the @p{initialize} function to + * be called before the object is used as preconditioner. + * + * Obsolete: * 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 @@ -72,10 +99,9 @@ * method. It is legal to apply this decomposition (@p{vmult} method) in * decomposed state. * - * * @sect2{Particular implementations} * - * It is enough to override the @p{decompose} and @p{vmult} methods to + * It is enough to override the @p{initialize} 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 @@ -88,7 +114,9 @@ * 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 + * @author Stephen "Cheffo" Kolaroff, 2002, based on SparseILU + * implementation by Wolfgang Bangerth; unified interface: Ralf + * Hartmann, 2003 */ template class SparseLUDecomposition : protected SparseMatrix @@ -96,41 +124,20 @@ 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. + * Constructor. Does nothing. * - * You have to initialize the - * matrix before usage with - * @p{reinit(SparsityPattern)}. + * Call the @p{initialize} + * function before using this + * object as preconditioner + * (@p{vmult}). */ - 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. + * This method is deprecated, and + * left for backward + * compability. It will be removed + * in later versions. */ SparseLUDecomposition (const SparsityPattern& sparsity); @@ -138,6 +145,14 @@ class SparseLUDecomposition : protected SparseMatrix * Destruction. */ virtual ~SparseLUDecomposition (); + + /** + * Deletes all member + * variables. Leaves the class in + * the state that it had directly + * after calling the constructor + */ + virtual void clear(); /** * Parameters for @@ -147,80 +162,152 @@ class SparseLUDecomposition : protected SparseMatrix { public: /** - * Constructor. + * Constructor. For the + * parameters' description, + * see below. */ - AdditionalData (const double strengthen_diagonal=0); + AdditionalData (const double strengthen_diagonal=0, + const unsigned int extra_off_diagonals=0, + const bool use_previous_sparsity=false, + const SparsityPattern *use_this_sparsity=0); /** - * strengthen_diag parameter. + * @p{strengthen_diag} times + * the sum of absolute row + * entries is added to the + * diagonal entries. + * + * Per default, this value is + * zero, i.e. the diagonal is + * not strengthened. */ double strengthen_diagonal; + + /** + * By default, the + * @p{initialize(matrix, + * data)} function creates + * its own sparsity. This + * sparsity has the same + * @p{SparsityPattern} as + * @p{matrix} with some extra + * off diagonals the number + * of which is specified by + * @p{extra_off_diagonals}. + * + * The user can give a + * @p{SparsityPattern} to + * @p{use_this_sparsity}. Then + * this sparsity is used and + * the + * @p{extra_off_diagonals} + * argument is ignored. + */ + unsigned int extra_off_diagonals; + + /** + * If this flag is true the + * @p{initialize} function uses + * the same sparsity that was + * used during the previous + * @p{initialize} call. + * + * This might be useful when + * several linear problems on + * the same sparsity need to + * solved, as for example + * several Newton iteration + * steps on the same + * triangulation. + */ + bool use_previous_sparsity; + + /** + * When a + * @ref{SparsityPattern} is + * given to this argument, + * the @p{initialize} + * function calls + * @p{reinit(*use_this_sparsity)} + * causing this sparsity to + * be used. + * + * Note that the sparsity + * structures of + * @p{*use_this_sparsity} and + * the matrix passed to the + * initialize 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. + */ + const SparsityPattern *use_this_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 needs to be + * called before an object of + * this class is used as + * preconditioner. * - * After this method is invoked, - * this object is out of synch - * (not decomposed state). + * For more detail about possible + * parameters, see the class + * documentation and the + * documentation of the + * @p{SparseLUDecomposition::AdditionalData} + * class. * - * This function only releases - * some memory and calls the - * respective function of the - * base class. + * According to the + * @p{parameters}, this function + * creates a new SparsityPattern + * or keeps the previous sparsity + * or takes the sparsity given by + * the user to @p{data}. Then, + * this function performs the LU + * decomposition. + * + * After this function is called + * the preconditioner is ready to + * be used (vmult). + */ + template + void initialize (const SparseMatrix &matrix, + const AdditionalData parameters); + + /** + * This method is deprecated, + * and left for backward + * compability. It will be removed + * in later versions. */ 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. + * This method is deprecated, + * and left for backward + * compability. It will be removed + * in later versions. */ void reinit (const SparsityPattern &sparsity); /** - * Perform the sparse LU - * factorization of the given - * matrix. After this method - * invocation, and before - * consecutive reinit invocation - * 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. + * This method is deprecated, + * and left for backward + * compability. It will be removed + * in later versions. */ template void decompose (const SparseMatrix &matrix, const double strengthen_diagonal=0.); /** - * Determines if this object is - * in synch with the underlying - * @ref{SparsityPattern}. + * This method is deprecated, + * and left for backward + * compability. It will be removed + * in later versions. */ virtual bool is_decomposed () const; @@ -321,7 +408,25 @@ class SparseLUDecomposition : protected SparseMatrix * array. */ void prebuild_lower_bound (); - + + /** + * In general this pointer is + * zero except for the case that + * no @p{SparsityPattern} is + * given to this class. Then, a + * @p{SparsityPattern} is created + * and is passed down to the + * @{SparseMatrix} base class. + * + * Nevertheless, the + * @p{SparseLUDecomposition} + * needs to keep ownership of + * this sparsity. It keeps this + * pointer to it enabling it to + * delete this sparsity at + * destruction time. + */ + SparsityPattern *own_sparsity; }; @@ -350,10 +455,15 @@ SparseLUDecomposition::is_decomposed () const template -inline SparseLUDecomposition::AdditionalData::AdditionalData ( - const double strengthen_diagonal): - strengthen_diagonal(strengthen_diagonal) + const double strengthen_diag, + const unsigned int extra_off_diag, + const bool use_prev_sparsity, + const SparsityPattern *use_this_spars): + strengthen_diagonal(strengthen_diag), + extra_off_diagonals(extra_off_diag), + use_previous_sparsity(use_prev_sparsity), + use_this_sparsity(use_this_spars) {} diff --git a/deal.II/lac/include/lac/sparse_decomposition.templates.h b/deal.II/lac/include/lac/sparse_decomposition.templates.h index d72f899a7c..2e6af4b844 100644 --- a/deal.II/lac/include/lac/sparse_decomposition.templates.h +++ b/deal.II/lac/include/lac/sparse_decomposition.templates.h @@ -18,7 +18,8 @@ template SparseLUDecomposition::SparseLUDecomposition() : SparseMatrix(), - decomposed(false) + decomposed(false), + own_sparsity(0) {} @@ -27,17 +28,94 @@ template SparseLUDecomposition:: SparseLUDecomposition (const SparsityPattern& sparsity) : SparseMatrix(sparsity), - decomposed(false) + decomposed(false), + own_sparsity(0) {} template SparseLUDecomposition::~SparseLUDecomposition() -{} +{ + clear(); +} + + +template +void SparseLUDecomposition::clear() +{ + decomposed = false; + + std::vector tmp; + tmp.swap (prebuilt_lower_bound); + + SparseMatrix::clear(); + + if (own_sparsity) + { + delete own_sparsity; + own_sparsity=0; + } +} +template +template +void SparseLUDecomposition::initialize ( + const SparseMatrix &matrix, + const AdditionalData data) +{ + const SparsityPattern &matrix_sparsity=matrix.get_sparsity_pattern(); + + if (data.use_this_sparsity) + reinit(*data.use_this_sparsity); + else if (data.use_previous_sparsity && + !this->empty() && + (this->m()==matrix.m())) + { + // Use the sparsity that was + // previously used. This is + // particularly useful when + // matrix entries change but + // not the sparsity, as for the + // case of several Newton + // iteration steps on an + // unchanged grid. + reinit(this->get_sparsity_pattern()); + } + else if (data.extra_off_diagonals==0) + { + // Use same sparsity as matrix + reinit(matrix_sparsity); + } + else + { + // Create new sparsity + + // for the case that + // own_sparsity wasn't deleted + // before (e.g. by clear()), do + // it here + if (own_sparsity) + { + // release the sparsity + SparseMatrix::clear(); + // delete it + delete own_sparsity; + } + + // and recreate + own_sparsity=new SparsityPattern(matrix_sparsity, + matrix_sparsity.max_entries_per_row() + +2*data.extra_off_diagonals, + data.extra_off_diagonals); + own_sparsity->compress(); + reinit(*own_sparsity); + } +} + + template template void diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index 3b2471dc92..78fc47993c 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -59,24 +59,28 @@ * corresponding ``.templates.h'' file and instantiate the respective * class yourself. * - * @author Wolfgang Bangerth, 1999, based on a similar implementation by Malte Braack + * @author Wolfgang Bangerth, 1999, based on a similar implementation + * by Malte Braack; unified interface: Ralf Hartmann */ template class SparseILU : public SparseLUDecomposition { public: /** - * Constructor. Does nothing, so - * you have to call @p{reinit} - * sometimes afterwards. + * Constructor. Does nothing. + * + * Call the @p{initialize} + * function before using this + * object as preconditioner + * (@p{vmult}). */ SparseILU (); /** - * Constructor. Initialize the - * sparsity pattern of this - * object with the given - * argument. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ SparseILU (const SparsityPattern &sparsity); @@ -95,48 +99,51 @@ class SparseILU : public SparseLUDecomposition * 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. + * This function needs to be + * called before an object of + * this class is used as + * preconditioner. + * + * For more details about + * possible parameters, see the + * class documentation of + * SparseLUDecomposition and the + * documentation of the + * @p{SparseLUDecomposition::AdditionalData} + * class. * - * If the - * @p{AdditionalData::strengthen_diagonal} - * parameter is greater than - * zero, this method invokes - * @p{get_strengthen_diagonal_impl - * ()}. + * According to the + * @p{parameters}, this function + * creates a new SparsityPattern + * or keeps the previous sparsity + * or takes the sparsity given by + * the user to @p{data}. Then, + * this function performs the LU + * decomposition. * - * Refer to - * @ref{SparseLUDecomposition} - * documentation for state - * management. + * After this function is called + * the preconditioner is ready to + * be used. */ template void initialize (const SparseMatrix &matrix, const AdditionalData parameters); /** - * Same as @p{initialize}. This method - * is deprecated, and left for - * backward compability. It may - * be removed in later versions. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ template void decompose (const SparseMatrix &matrix, const double strengthen_diagonal=0.); /** - * Same as @p{vmult}. This method - * is deprecated, and left for - * backward compability. It may - * be removed in later versions. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ template void apply_decomposition (Vector &dst, @@ -147,10 +154,8 @@ class SparseILU : public SparseLUDecomposition * i.e. do one forward-backward step * $dst=(LU)^{-1}src$. * - * Refer to - * @ref{SparseLUDecomposition} - * documentation for state - * management. + * The @p{initialize} function + * needs to be called beforehand. */ template void vmult (Vector &dst, @@ -197,14 +202,5 @@ SparseILU::apply_decomposition (Vector &dst, -template -template -inline -void SparseILU::decompose (const SparseMatrix &matrix, - const double strengthen_diagonal) -{ - initialize(matrix, AdditionalData(strengthen_diagonal)); -} - #endif // __deal2__sparse_ilu_h diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 6728c4dba5..f4109f6781 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -34,22 +34,36 @@ SparseILU::SparseILU (const SparsityPattern &sparsity) : + template template void SparseILU::initialize (const SparseMatrix &matrix, const AdditionalData data) { - SparseLUDecomposition::decompose (matrix, data.strengthen_diagonal); + SparseLUDecomposition::initialize(matrix, data); + + decompose(matrix, data.strengthen_diagonal); +} + + + +template +template +inline +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 (data.strengthen_diagonal>=0, - ExcInvalidStrengthening (data.strengthen_diagonal)); + Assert (strengthen_diagonal>=0, + ExcInvalidStrengthening (strengthen_diagonal)); this->copy_from (matrix); - if (data.strengthen_diagonal>0) + if (strengthen_diagonal>0) this->strengthen_diagonal_impl(); const SparsityPattern &sparsity = this->get_sparsity_pattern(); @@ -150,6 +164,7 @@ void SparseILU::initialize (const SparseMatrix &matrix, + template template void SparseILU::vmult (Vector &dst, diff --git a/deal.II/lac/include/lac/sparse_mic.h b/deal.II/lac/include/lac/sparse_mic.h index 9298f29abd..f5b7c34588 100644 --- a/deal.II/lac/include/lac/sparse_mic.h +++ b/deal.II/lac/include/lac/sparse_mic.h @@ -29,7 +29,8 @@ * where X is a diagonal matrix, defined by the condition rowsum(A) = * rowsum(B). * - * @author Stephen "Cheffo" Kolaroff, 2002. + * @author Stephen "Cheffo" Kolaroff, 2002, unified interface: Ralf + * Hartmann 2003. */ template class SparseMIC : public SparseLUDecomposition @@ -49,7 +50,20 @@ class SparseMIC : public SparseLUDecomposition * argument. */ SparseMIC (const SparsityPattern &sparsity); - + + /** + * Destruction. + */ + virtual ~SparseMIC(); + + /** + * Deletes all member + * variables. Leaves the class in + * the state that it had directly + * after calling the constructor + */ + virtual void clear(); + /** * Make the @p{AdditionalData} * type in the base class @@ -61,38 +75,18 @@ class SparseMIC : public SparseLUDecomposition AdditionalData; /** - * 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. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ 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. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ void reinit (const SparsityPattern &sparsity); @@ -104,31 +98,10 @@ class SparseMIC : public SparseLUDecomposition const AdditionalData parameters); /** - * 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 the - * @p{strengthen_diagonal_impl} - * function of the base class. - * - * Refer to - * @ref{SparseLUDecomposition} - * documentation for state - * management. + * This method is deprecated, and + * left for backward + * compability. It will be + * removed in later versions. */ template void decompose (const SparseMatrix &matrix, @@ -139,10 +112,8 @@ class SparseMIC : public SparseLUDecomposition * i.e. do one forward-backward step * $dst=(LU)^{-1}src$. * - * Refer to - * @ref{SparseLUDecomposition} - * documentation for state - * management. + * Call @p{initialize} before + * calling this function. */ template void vmult (Vector &dst, @@ -214,15 +185,4 @@ class SparseMIC : public SparseLUDecomposition -template -template -inline -void SparseMIC::initialize (const SparseMatrix &matrix, - const AdditionalData data) -{ - decompose(matrix, data.strengthen_diagonal); -} - - - #endif // __deal2__ diff --git a/deal.II/lac/include/lac/sparse_mic.templates.h b/deal.II/lac/include/lac/sparse_mic.templates.h index a60b9944fb..a2f66bd22d 100644 --- a/deal.II/lac/include/lac/sparse_mic.templates.h +++ b/deal.II/lac/include/lac/sparse_mic.templates.h @@ -37,6 +37,48 @@ SparseMIC::SparseMIC (const SparsityPattern &sparsity) {} +template +SparseMIC::~SparseMIC() +{ + clear(); +} + + +template +void SparseMIC::clear() +{ + 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::clear(); +} + + +template +template +inline +void SparseMIC::initialize (const SparseMatrix &matrix, + const AdditionalData data) +{ + SparseLUDecomposition::initialize(matrix, data); + + decompose(matrix, data.strengthen_diagonal); +} + + template void -- 2.39.5