]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add initialize and clear functions in order to unify the interface of preconditioners...
authorhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Feb 2003 17:39:03 +0000 (17:39 +0000)
committerhartmann <hartmann@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Feb 2003 17:39:03 +0000 (17:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@7183 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_decomposition.h
deal.II/lac/include/lac/sparse_decomposition.templates.h
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.templates.h
deal.II/lac/include/lac/sparse_mic.h
deal.II/lac/include/lac/sparse_mic.templates.h

index 04f09bf192a46f320dd3d3aaafcb6893574d3533..8877120d9e6f811147b278ba2fdf2cb170c26a16 100644 (file)
  * 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<double> lu (lu_sparsity);
- *   lu.decompose (global_matrix);
+ * SparseLUImplementation<double> lu;
+ * lu.initialize(matrix, SparseLUImplementation<double>::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
  * 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
  * 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 <typename number>
 class SparseLUDecomposition : protected SparseMatrix<number>
@@ -96,41 +124,20 @@ 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.
+                                     * 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<number>
                                       * 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<number>
     {
       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 <typename somenumber>
+    void initialize (const SparseMatrix<somenumber> &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 <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}.
+                                     * 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<number>
                                       * 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<number>::is_decomposed () const
 
 
 template <typename number>
-inline
 SparseLUDecomposition<number>::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)
 {}
 
 
index d72f899a7c2517731406877fbf1a7457dec01868..2e6af4b8440202c0eed3e4c3d4400a6f91afc5db 100644 (file)
@@ -18,7 +18,8 @@ template<typename number>
 SparseLUDecomposition<number>::SparseLUDecomposition()
                 :
                 SparseMatrix<number>(),
-                decomposed(false)
+                decomposed(false),
+                own_sparsity(0)
 {}
 
 
@@ -27,17 +28,94 @@ template<typename number>
 SparseLUDecomposition<number>::
 SparseLUDecomposition (const SparsityPattern& sparsity) :
                 SparseMatrix<number>(sparsity),
-                decomposed(false)
+                decomposed(false),
+                own_sparsity(0)
 {}
 
 
 
 template<typename number>
 SparseLUDecomposition<number>::~SparseLUDecomposition()
-{}
+{
+  clear();
+}
+
+
+template<typename number>
+void SparseLUDecomposition<number>::clear()
+{
+  decomposed = false;
+  
+  std::vector<const unsigned int*> tmp;
+  tmp.swap (prebuilt_lower_bound);
+
+  SparseMatrix<number>::clear();
+  
+  if (own_sparsity)
+    {
+      delete own_sparsity;
+      own_sparsity=0;
+    }
+}
 
 
 
+template<typename number>
+template <typename somenumber>
+void SparseLUDecomposition<number>::initialize (
+  const SparseMatrix<somenumber> &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<number>::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<typename number>
 template<typename somenumber>
 void
index 3b2471dc92a7ce1b6511fd6807760f748166e107..78fc47993c15b325d754663a50f4b12ffa6bace0 100644 (file)
  * 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 <typename number>
 class SparseILU : public SparseLUDecomposition<number>
 {
   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<number>
                                      * 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 <typename somenumber>
     void initialize (const SparseMatrix<somenumber> &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 <typename somenumber>
     void decompose (const SparseMatrix<somenumber> &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 <typename somenumber>
     void apply_decomposition (Vector<somenumber>       &dst,
@@ -147,10 +154,8 @@ class SparseILU : public SparseLUDecomposition<number>
                                      * 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 <typename somenumber>
     void vmult (Vector<somenumber>       &dst,
@@ -197,14 +202,5 @@ SparseILU<number>::apply_decomposition (Vector<somenumber>       &dst,
 
 
 
-template <typename number>
-template <typename somenumber>
-inline
-void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
-                                  const double strengthen_diagonal)
-{
-  initialize(matrix, AdditionalData(strengthen_diagonal));
-}
-
 
 #endif // __deal2__sparse_ilu_h
index 6728c4dba5d18ffb52a4bc9599a3b20cbc993b2e..f4109f6781950229e405215eae744b93d4f2dfe7 100644 (file)
@@ -34,22 +34,36 @@ SparseILU<number>::SparseILU (const SparsityPattern &sparsity) :
 
 
 
+
 template <typename number>
 template <typename somenumber>
 void SparseILU<number>::initialize (const SparseMatrix<somenumber> &matrix,
                                    const AdditionalData data)
 {
-  SparseLUDecomposition<number>::decompose (matrix, data.strengthen_diagonal);
+  SparseLUDecomposition<number>::initialize(matrix, data);
+      
+  decompose(matrix, data.strengthen_diagonal);
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+inline
+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 (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<number>::initialize (const SparseMatrix<somenumber> &matrix,
 
 
 
+
 template <typename number>
 template <typename somenumber>
 void SparseILU<number>::vmult (Vector<somenumber>       &dst,
index 9298f29abde0ef9d64908d5d8dc7b3e4ba2fed2a..f5b7c345880c9805f140b9e5737d4bd29acda29e 100644 (file)
@@ -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 <typename number>
 class SparseMIC : public SparseLUDecomposition<number>
@@ -49,7 +50,20 @@ class SparseMIC : public SparseLUDecomposition<number>
                                       * 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<number>
     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<number>
                     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 <typename somenumber>
     void decompose (const SparseMatrix<somenumber> &matrix,
@@ -139,10 +112,8 @@ class SparseMIC : public SparseLUDecomposition<number>
                                      * 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 <typename somenumber>
     void vmult (Vector<somenumber>       &dst,
@@ -214,15 +185,4 @@ class SparseMIC : public SparseLUDecomposition<number>
 
 
 
-template <typename number>
-template <typename somenumber>
-inline
-void SparseMIC<number>::initialize (const SparseMatrix<somenumber> &matrix,
-                                   const AdditionalData data)
-{
-  decompose(matrix, data.strengthen_diagonal);
-}
-
-
-
 #endif  // __deal2__
index a60b9944fbffa8acb5582201d9619d98a43669ce..a2f66bd22dc93cf32dea6bebbbec154d5c12bb9d 100644 (file)
@@ -37,6 +37,48 @@ SparseMIC<number>::SparseMIC (const SparsityPattern &sparsity)
 {}
 
 
+template <typename number>
+SparseMIC<number>::~SparseMIC()
+{
+  clear();
+}
+
+
+template <typename number>
+void SparseMIC<number>::clear()
+{
+  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>::clear();
+}
+
+
+template <typename number>
+template <typename somenumber>
+inline
+void SparseMIC<number>::initialize (const SparseMatrix<somenumber> &matrix,
+                                   const AdditionalData data)
+{
+  SparseLUDecomposition<number>::initialize(matrix, data);
+
+  decompose(matrix, data.strengthen_diagonal);
+}
+
+
 
 template <typename number>
 void

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.