]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Doc updates and small formatting things.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 May 1999 12:09:45 +0000 (12:09 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 May 1999 12:09:45 +0000 (12:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@1287 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/Todo
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/include/lac/ivector.h
deal.II/lac/include/lac/precondition_block.h
deal.II/lac/include/lac/precondition_block.templates.h
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/vector.h
deal.II/lac/include/lac/vector_memory.h
deal.II/lac/source/Makefile
deal.II/lac/source/sparse_matrix.cc

index 21b97904c7539ab828126304a8d14a418fcab863..a706b31c4867e60bcc1c647ddbb87e0b8fcb4ff6 100644 (file)
@@ -29,4 +29,14 @@ Make ILU decomposition faster by following the comment in the
   innermost loop (about the better use of the index j).
 
 Change in-class defined SolverSelector functions out-of-class to
-  conform to general style guide. (Ralf?)
\ No newline at end of file
+  conform to general style guide. (Ralf?)
+
+Eliminate #vec_len# from the SparseMatrixStruct class. This should
+  be rather simple since as far as I can see it is only used once
+  (at all other places, it is only written) and also its value
+  should be equal to rowstart[row].
+
+Why not replace Vector::scale by Vector::operator *= ?
+
+Use the commented-out version in PreconditionBlock::invert_diagblocks
+  using the try-catch clauses
\ No newline at end of file
index 732efc0b345c3b29627097deeed507927a63f2d7..523f6c6708662e76d93e4fc90680564fc069568c 100644 (file)
@@ -26,18 +26,20 @@ class iVector;
  * Rectangular/quadratic full matrix.
  *
  * Implementation of a classical rectangular scheme of numbers. The
- * data type of the entries is provided in the template argument #number#.
- * The interface is quite fat and in fact has grown every time a new
- * feature was needed. So, a lot of functions are provided.
+ * data type of the entries is provided in the template argument
+ * #number#.  The interface is quite fat and in fact has grown every
+ * time a new feature was needed. So, a lot of functions are provided.
  *
  * Since the instantiation of this template is quite an effort,
- * standard versions are precompiled into the library. These include all
- * combinations of 'float' and 'double' for matrices and vectors. If you need more
- * data types, the implementation of non-inline functions is in
- * "fullmatrix.templates.h". Driver files are in the source tree.
+ * standard versions are precompiled into the library. These include
+ * all combinations of #float# and #double# for matrices and
+ * vectors. If you need more data types, the implementation of
+ * non-inline functions is in #fullmatrix.templates.h#. Driver files
+ * are in the source tree.
  *
- * Internal calculations are usually done with the accuracy of the vector argument to
- * functions. If there is no argument with a number type, the matrix number type is used.
+ * Internal calculations are usually done with the accuracy of the
+ * vector argument to functions. If there is no argument with a number
+ * type, the matrix number type is used.
  *
  * <TABLE BORDER=1>
  * <TR><TH ALIGN=CENTER><B>this</B><TH ALIGN=CENTER><B>other
@@ -52,67 +54,27 @@ class iVector;
  * <CAPTION>Instantiations provided in the library</CAPTION>
  * </TABLE>
  *
- *  CONVENTIONS for used 'equations' : <p>
- *  - THIS matrix is always named 'A' <p>
- *  - matrices are always uppercase , vectors and scalars are lowercase <p>
- *  - Transp(A) used for transpose of matrix A
+ * In the documentation of member functions, the following conventions are adopted:
+ * \begin{itemize}
+ *  \item THIS matrix is always named #A#.
+ *  \item Matrices are always uppercase, vectors and scalars are lowercase.
+ *  \item #Transp(A)# denotes the transpose of matrix A.
+ * \end{itemize}
  *
  * @author Guido Kanschat, Franz-Theo Suttmeier, Wolfgang Bangerth
  */
 template<typename number>
 class FullMatrix : public Subscriptor
 {
-  private:
-                                    /**
-                                     * Component-array.
-                                     */
-    number* val;
-                                    /** 
-                                     * Dimension of range. Actual number of Columns
-                                     */
-    unsigned int dim_range;
-                                    /**
-                                     * Dimension of image. Actual number of Rows
-                                     */
-    unsigned int dim_image;
-                                    /**
-                                     * Dimension. Determines amount of reserved memory
-                                     */
-    unsigned int val_size;
-    
-                                    /**
-                                     * Initialization. Initialize
-                                     * memory for a #FullMatrix#
-                                     * of #m# rows and #n#
-                                     * columns to zero.
-                                     */
-    void init (const unsigned int m, const unsigned int n);
-    
-                                    /**
-                                     * Return a read-write reference to the
-                                     * element #(i,j)#.
-                                     *
-                                     * This function does no bounds
-                                     * checking and is only to be used
-                                     * internally and in functions
-                                     * already checked.
-                                     */
-    number& el (const unsigned int i, const unsigned int j);
-    
-                                    /**
-                                     * Return the value of the element #(i,j)#.
-                                     *
-                                     * This function does no bounds checking and is only to be used
-                                     * internally and in functions
-                                     * already checked.
-                                     */
-    number el (const unsigned int i, const unsigned int j) const;
-    
-    
   public:
                                     /**
                                      * Constructor. Initialize the matrix as
                                      * a square matrix with dimension #n#.
+                                     *
+                                     * In order to avoid the implicit
+                                     * conversion of integers and other types
+                                     * to a matrix, this constructor is
+                                     * declared #explicit#.
                                      */
     explicit FullMatrix (const unsigned int n = 1);
     
@@ -172,38 +134,40 @@ class FullMatrix : public Subscriptor
                                      * Fill rectangular block.
                                      *
                                      * The matrix #src# is copied
-                                     into the target. The optional
-                                     values #i# and #j# determine the
-                                     upper left corner of the image
-                                     of #src#.
+                                     into the target. The optional
+                                     values #i# and #j# determine the
+                                     upper left corner of the image
+                                     of #src#.
                                      *
                                      * This function requires that
-                                     #i+src.m()<=m()# and
-                                     #j+src.n()<=n()#, that is, the
-                                     image fits into the space of #this#.
+                                     #i+src.m()<=m()# and
+                                     #j+src.n()<=n()#, that is, the
+                                     image fits into the space of #this#.
                                      */
     template<typename number2>
-    void fill (const FullMatrix<number2>& src,
-              const unsigned int i=0, const unsigned int j=0);
+    void fill (const FullMatrix<number2> &src,
+              const unsigned int         i=0,
+              const unsigned int         j=0);
     
                                     /**
-                                     * Change  Dimensions.
-                                     * Set dimension to (m,n) <p>
-                                     * ( reinit rectangular matrix )
+                                     * Set dimension to $m\times n$ and
+                                     * allocate memory if necessary. Forget
+                                     * the previous content of the matrix.
                                      */
-    void reinit (const unsigned int m, const unsigned int n);
+    void reinit (const unsigned int m,
+                const unsigned int n);
     
                                     /**
-                                     * Change  Dimensions.
-                                     * Set dimension to (n,n) <p>
-                                     * ( reinit quadratic matrix )
+                                     * Set dimension to $n\times n$ and
+                                     * allocate memory if necessary. Forget
+                                     * the previous content of the matrix.
                                      */
     void reinit (const unsigned int n);
     
                                     /**
-                                     * Adjust  Dimension.
-                                     * Set dimension to ( m(B),n(B) ) <p>
-                                     * ( adjust to dimensions of another matrix B )
+                                     * Set dimension to $m(B)\times n(B)$ and
+                                     * allocate memory if necessary. Forget
+                                     * the previous content of the matrix.
                                      */
     template<typename number2>
     void reinit (const FullMatrix<number2> &B);
@@ -226,66 +190,64 @@ class FullMatrix : public Subscriptor
                                      * Return whether the matrix contains only
                                      * elements with value zero. This function
                                      * is mainly for internal consistency
-                                     * check and should seldomly be used when
+                                     * checks and should seldomly be used when
                                      * not in debug mode since it uses quite
                                      * some time.
                                      */
     bool all_zero () const;
 
-                                    /*
-                                     *   Access Elements. returns element at relative 'address' i <p>
-                                     *   ( -> access to A(i/n , i mod n) )
-                                     */
-//    number el (const unsigned int i) const;
-    
                                     /**
                                      * Return the value of the element #(i,j)#.
-                                     * Does the same as the private #el(i,j)# function
-                                     * but does bounds checking in
+                                     * Does the same as the private #el(i,j)#
+                                     * function but does bounds checking in
                                      * debug mode.
                                      */
-    number operator() (const unsigned int i, const unsigned int j) const;
+    number operator() (const unsigned int i,
+                      const unsigned int j) const;
     
                                     /**
                                      * Return a read-write reference to
                                      * the element #(i,j)#.
-                                     * Does the same as the private #el(i,j)# function
-                                     * but does bounds checking in
+                                     * Does the same as the private #el(i,j)#
+                                     * function but does bounds checking in
                                      * debug mode.
                                      */
-    number& operator() (const unsigned int i, const unsigned int j);
+    number& operator() (const unsigned int i,
+                       const unsigned int j);
     
                                     /**
                                      * Set all entries in the matrix to
-                                     * zero.
+                                     * zero. Do not resize the matrix.
                                      */
     void clear ();
 
                                     /**
-                                     *  Weighted addition. The matrix
-                                     #s*B# is added to #this#.
+                                     * Weighted addition. The matrix
+                                     #s*B# is added to #this#.
                                      *
                                      * $A += sB$
                                      */
     template<typename number2>
-    void add (const number s, const FullMatrix<number2>& B);
+    void add (const number               s,
+             const FullMatrix<number2> &B);
 
                                     /**
                                      * Weighted addition of the
-                                     transpose of #B# to #this#.
+                                     transpose of #B# to #this#.
                                      *
                                      * $A += s B^T$
                                      */
     template<typename number2>
-    void Tadd (const number s, const FullMatrix<number2>& B);
+    void Tadd (const number               s,
+              const FullMatrix<number2> &B);
     
                                     /**
                                      * Matrix-matrix-multiplication.
                                      * $C=A*B$.
                                      */
     template<typename number2>
-    void mmult (FullMatrix<number2>& C, const FullMatrix<number2>& B) const;
+    void mmult (FullMatrix<number2>       &C,
+               const FullMatrix<number2> &B) const;
     
                                     /**
                                      * Matrix-matrix-multiplication using
@@ -293,30 +255,36 @@ class FullMatrix : public Subscriptor
                                      * $C=A^T*B$.
                                      */
     template<typename number2>
-    void Tmmult (FullMatrix<number2>& C, const FullMatrix<number2>& B) const;
+    void Tmmult (FullMatrix<number2>       &C,
+                const FullMatrix<number2> &B) const;
     
                                     /**
                                      * Matrix-vector-multiplication.
                                      *
                                      * The optional parameter
                                      * #adding# determines, whether the
-                                     * result is stored in #w# or addet
+                                     * result is stored in #w# or added
                                      * to #w#.
                                      *
                                      * if (adding)
-                                     *  w += A*v
+                                     *  $w += A*v$
                                      *
                                      * if (!adding)
-                                     *  w = A*v
+                                     *  $w = A*v$
                                      */
     template<typename number2>
-    void vmult (Vector<number2>& w, const Vector<number2>& v, const bool adding=false) const;
+    void vmult (Vector<number2>       &w,
+               const Vector<number2> &v,
+               const bool             adding=false) const;
     
                                     /**
-                                     * Transpose matrix-vector-multiplication. See #vmult# above.
+                                     * Transpose matrix-vector-multiplication.
+                                     * See #vmult# above.
                                      */
     template<typename number2>
-    void Tvmult (Vector<number2>& w, const Vector<number2>& v, const bool adding=false) const;
+    void Tvmult (Vector<number2>       &w,
+                const Vector<number2> &v,
+                const bool             adding=false) const;
 
                                     /**
                                      * Return the norm of the vector #v# with
@@ -350,38 +318,48 @@ class FullMatrix : public Subscriptor
                                      * the finite element context.
                                      */
     template<typename number2>
-    double matrix_scalar_product (const Vector<number2> &u, const Vector<number2> &v) const;
+    double matrix_scalar_product (const Vector<number2> &u,
+                                 const Vector<number2> &v) const;
 
                                     /**
-                                     * Return the l1-norm of the matrix, i.e.
+                                     * Return the $l_1$-norm of the matrix, i.e.
                                      * $|M|_1=max_{all columns j}\sum_{all 
                                      * rows i} |M_ij|$,
-                                     * (max. sum of columns).
-                                     * This is the
+                                     * (max. sum of columns). This is the
                                      * natural matrix norm that is compatible
-                                     * to the l1-norm for vectors, i.e.
+                                     * to the $l_1$-norm for vectors, i.e.
                                      * $|Mv|_1\leq |M|_1 |v|_1$.
                                      * (cf. Rannacher Numerik0)
                                      */
     number l1_norm () const;
 
                                     /**
-                                     * Return the linfty-norm of the
+                                     * Return the $l_\infty$-norm of the
                                      * matrix, i.e.
-                                     * $|M|_infty=max_{all rows i}\sum_{all 
-                                     * columns j} |M_ij|$,
+                                     * $|M|_\infty=\max_{all rows i}\sum_{all 
+                                     * columns j} |M_{ij}|$,
                                      * (max. sum of rows).
                                      * This is the
                                      * natural matrix norm that is compatible
-                                     * to the linfty-norm of vectors, i.e.
-                                     * $|Mv|_infty \leq |M|_infty |v|_infty$.
+                                     * to the $l_\infty$-norm of vectors, i.e.
+                                     * $|Mv|_\infty \leq |M|_\infty |v|_\infty$.
                                      * (cf. Rannacher Numerik0)
                                      */
     number linfty_norm () const;
     
+                                    /**
+                                     * Compute the quadratic matrix norm.
+                                     * Return value is the root of the square
+                                     * sum of all matrix entries.
+                                     */
+    number norm2 () const;
+    
                                     /**
                                      * A=Inverse(A). Inversion of this by
-                                     * Gauss-Jordan-algorithm
+                                     * Gauss-Jordan-algorithm. Note that this
+                                     * is a rather expensive operation, so
+                                     * you may not want to use it for
+                                     * larger matrices if not necessary.
                                      */
     void gauss_jordan ();
 
@@ -395,45 +373,59 @@ class FullMatrix : public Subscriptor
                                       */
     double determinant () const;
 
-                                    /**
-                                     * Compute the quadratic matrix norm.
-                                     * Return value is the root of the square
-                                     * sum of all matrix entries.
-                                     */
-    double norm2 () const;
-
                                     /**
                                      * Assign the inverse of the given
                                      * matrix to #*this#. This function is
                                      * only implemented (hardcoded) for
                                      * square matrices of dimension one,
-                                     * two and three.
+                                     * two, three and four, since the
+                                     * amount of code needed grows quickly.
+                                     * The implementation does not use
+                                     * an elimination method like the
+                                     * Gauss-Jordan one, but rather sets
+                                     * the element directly; their values
+                                     * are precomputed symbolically using
+                                     * Maple. This way, we can avoid the
+                                     * overhead of loops and local variables
+                                     * but the number of lines of code
+                                     * grows rapidly.
+                                     *
+                                     * For all other sizes than the ones given
+                                     * above, an exception of type
+                                     * #ExcNotImplemented(dim_range)# is
+                                     * thrown, which you can catch and use
+                                     * some other method to invert the matrix,
+                                     * e.g. the #gauss_jordan# function.
                                      */
     void invert (const FullMatrix<number> &M);
 
                                     /**
-                                     *  A(i,1-n)+=s*A(j,1-n).
+                                     * $A(i,1-n)+=s*A(j,1-n)$.
                                      * Simple addition of rows of this
                                      */
-    void add_row (const unsigned int i, const number s, const unsigned int j);
+    void add_row (const unsigned int i,
+                 const number       s,
+                 const unsigned int j);
 
                                     /**
-                                     *  A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n).
-                                     *  Multiple addition of rows of this
+                                     * $A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n)$.
+                                     * Multiple addition of rows of this.
                                      */
     void add_row (const unsigned int i,
                  const number s, const unsigned int j,
                  const number t, const unsigned int k);
 
                                     /**
-                                     *  A(1-n,i)+=s*A(1-n,j).
-                                     *  Simple addition of columns of this
+                                     * $A(1-n,i)+=s*A(1-n,j)$.
+                                     *  Simple addition of columns of this.
                                      */
-    void add_col (const unsigned int i, const number s, const unsigned int j);
+    void add_col (const unsigned int i,
+                 const number       s,
+                 const unsigned int j);
 
                                     /**
-                                     *  A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k).
-                                     *  Multiple addition of columns of this
+                                     * $A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k)$.
+                                     *  Multiple addition of columns of this.
                                      */
     void add_col (const unsigned int i,
                  const number s, const unsigned int j,
@@ -452,11 +444,14 @@ class FullMatrix : public Subscriptor
     void swap_col (const unsigned int i, const unsigned int j);
 
                                     /**
-                                     *  w=b-A*v.
-                                     *  Residual calculation , returns |w|
+                                     * $w=b-A*v$.
+                                     * Residual calculation , returns
+                                     * the $l_2$-norm $|w|$
                                      */
     template<typename number2, typename number3>
-    double residual (Vector<number2>& w, const Vector<number2>& v, const Vector<number3>& b) const;
+    double residual (Vector<number2>      & w,
+                    const Vector<number2>& v,
+                    const Vector<number3>& b) const;
 
                                     /**
                                      * Forward elimination of lower triangle.
@@ -470,31 +465,37 @@ class FullMatrix : public Subscriptor
                                      * is considered
                                      */
     template<typename number2>
-    void forward (Vector<number2>& dst, const Vector<number2>& src) const;
+    void forward (Vector<number2>       &dst,
+                 const Vector<number2> &src) const;
 
                                     /**
                                      * Backward elimination of upper triangle.
                                      * @see forward
                                      */
     template<typename number2>
-    void backward (Vector<number2>& dst, const Vector<number2>& src) const;
+    void backward (Vector<number2>       &dst,
+                  const Vector<number2> &src) const;
 
                                     /**
-                                     * QR - factorization of a matrix.
+                                     * QR-factorization of a matrix.
                                      * The orthogonal transformation Q is
-                                     * applied to the vector y and this matrix. <p>
+                                     * applied to the vector y and this matrix.
+                                     *
                                      * After execution of householder, the upper
-                                     *  triangle contains the resulting matrix R, <p>
-                                     * the lower the incomplete factorization matrices.
+                                     * triangle contains the resulting matrix R,
+                                     * the lower the incomplete factorization
+                                     * matrices.
                                      */
     template<typename number2>
-    void householder (Vector<number2>y);
+    void householder (Vector<number2> &y);
 
                                     /**
-                                     * Least - Squares - Approximation by QR-factorization.
+                                     * Least-Squares-Approximation by
+                                     * QR-factorization.
                                      */
     template<typename number2>
-    double least_squares (Vector<number2>& dst, Vector<number2>& src);
+    double least_squares (Vector<number2> &dst,
+                         Vector<number2> &src);
 
                                     /**
                                      *  A(i,i)+=B(i,1-n). Addition of complete
@@ -502,7 +503,8 @@ class FullMatrix : public Subscriptor
                                      *  ( i = 1 ... m )
                                      */
     template<typename number2>
-    void add_diag (const number s, const FullMatrix<number2>& B);
+    void add_diag (const number               s,
+                  const FullMatrix<number2> &B);
 
                                     /**
                                      *  A(i,i)+=s  i=1-m.
@@ -592,6 +594,70 @@ class FullMatrix : public Subscriptor
                                      * Exception
                                      */
     DeclException0 (ExcIO);
+
+  private:
+                                    /**
+                                     * Component-array.
+                                     */
+    number *val;
+    
+                                    /** 
+                                     * Dimension of range.
+                                     * Actual number of Columns
+                                     */
+    unsigned int dim_range;
+    
+                                    /**
+                                     * Dimension of image. Actual number of Rows
+                                     */
+    unsigned int dim_image;
+    
+                                    /**
+                                     * Dimension of the array
+                                     * holding the values of the
+                                     * matrix elements. Determines
+                                     * amount of reserved memory.
+                                     *
+                                     * Actually, the allocated array may
+                                     * not have a size equal to the number
+                                     * of elements of this matrix, since
+                                     * reallocation only happens when the
+                                     * size of the matrix is increased.
+                                     * Therefore, if the matrix size was
+                                     * decreased somewhen in the past,
+                                     * #val_size# will be larger than
+                                     * #dim_range * dim_image#.
+                                     */
+    unsigned int val_size;
+    
+                                    /**
+                                     * Initialization. Initialize
+                                     * memory for a #FullMatrix#
+                                     * of #m# rows and #n#
+                                     * columns to zero.
+                                     */
+    void init (const unsigned int m, const unsigned int n);
+    
+                                    /**
+                                     * Return a read-write reference to the
+                                     * element #(i,j)#.
+                                     *
+                                     * This function does no bounds
+                                     * checking and is only to be used
+                                     * internally and in functions
+                                     * already checked.
+                                     */
+    number & el (const unsigned int i, const unsigned int j);
+    
+                                    /**
+                                     * Return the value of the element #(i,j)#.
+                                     *
+                                     * This function does no bounds checking
+                                     * and is only to be used
+                                     * internally and in functions
+                                     * already checked.
+                                     */
+    number el (const unsigned int i, const unsigned int j) const;    
 };
 
 
index f6193a393d328948d08833c4cc48844db911b15a..998cf0ad771900f3f308e061694c09bd01776980 100644 (file)
@@ -71,7 +71,8 @@ FullMatrix<number>::all_zero () const
 
 template <typename number>
 void
-FullMatrix<number>::reinit (const unsigned int mm, const unsigned int nn)
+FullMatrix<number>::reinit (const unsigned int mm,
+                           const unsigned int nn)
 {
   if (val_size<nn*mm)
     {
@@ -405,14 +406,13 @@ double FullMatrix<number>::residual (Vector<number2>& dst,
   Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
   Assert(right.size() == m(), ExcDimensionMismatch(right.size(), m()));
 
-  unsigned int i,j;
   number2 s, res = 0.;
   const unsigned int size_m = m(),
                     size_n = n();
-  for (i=0; i<size_n; ++i)
+  for (unsigned int i=0; i<size_n; ++i)
     {
       s = right(i);
-      for (j=0; j<size_m; ++j)
+      for (unsigned int j=0; j<size_m; ++j)
        s -= src(j) * el(i,j);
       dst(i) = s;
       res += s*s;
@@ -1185,7 +1185,7 @@ FullMatrix<number>::determinant () const
 
 
 template <typename number>
-double
+number
 FullMatrix<number>::norm2 () const
 {
   number s = 0.;
index cc2d22ee48bba0beba28dad5444d2f9291291c06..d263cc2b006a57172267bd21769e8c861d08d1bd 100644 (file)
@@ -26,7 +26,7 @@ class iVector
 {
   friend class dFMatrix;
 
-protected:
+  protected:
 
                                     /// Dimension. Actual number of components
     unsigned int dim;
index 26e9e0f8ad0749e657753ad0a300e26958f20585..8eb6353cc013dc876fc60a0581e0a86780398b7a 100644 (file)
  * Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ...
  * This class assumes the #SparseMatrix<number># consisting of invertible blocks 
  * of #blocksize# on the diagonal and provides the inversion of the diagonal blocks
- * of the matrix. NOT only diagonal block matrices are allowed but all
- * matrices of arbitrary structure with the minimal property of having got
+ * of the matrix. NOT only block diagonal matrices are allowed but all
+ * matrices of arbitrary structure with the minimal property of having
  * invertible blocks on the diagonal!
  *
  * This block matrix structure is given e.g. for the DG method
  * for the transport equation. For a downstream numbering the matrices
- * even have got a lower block diagonal matrix structure, i.e. the matrices
+ * even have got a block lower left matrix structure, i.e. the matrices
  * are empty above the diagonal blocks.
  *
  * For all matrices that are empty above and below the diagonal
  * blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner
  * is a direct solver. For all matrices that are empty only above the diagonal blocks
- * (e.g. the matrices one get by the DG method with downstream numbering) the
+ * (e.g. the matrices one gets by the DG method with downstream numbering) the
  * #BlockSOR# is a direct solver.
  * 
  * This first implementation of the #PreconditionBlock# assumes the
- * matrix having blocks each of the same block size. Varying
+ * matrix has blocks each of the same block size. Varying
  * block sizes within the matrix must still be implemented if needed.
  *
  * The first template parameter denotes the type of number representation in
  * the sparse matrix, the second denotes the type of number representation in
- * which the inverse diagonal block
- * matrices are stored by #invert_diagblocks()#.
+ * which the inverted diagonal block matrices are stored within this class
+ * by #invert_diagblocks()#. If you don't want to use the block inversion as
+ * an exact solver, but rather as a preconditioner, you may probably want to
+ * store the inverted blocks with less accuracy than the original matrix;
+ * for example, #number==double, inverse_type=float# might be a viable choice.
  */
 template<typename number, typename inverse_type>
 class PreconditionBlock: public Subscriptor
@@ -66,7 +69,10 @@ class PreconditionBlock: public Subscriptor
 
                                     /**
                                      * Takes the matrix that should be used
-                                     * for the preconditioning.
+                                     * for the preconditioning. A reference
+                                     * to it is stored within this class,
+                                     * but ownership of the matrix remains
+                                     * with the caller of this function.
                                      */
     void use_matrix(const SparseMatrix<number> &M);
     
@@ -85,7 +91,7 @@ class PreconditionBlock: public Subscriptor
 
                                     /**
                                      * Stores the inverse of
-                                     * the diagonal blocks matrices
+                                     * the diagonal blocks
                                      * in #inverse#. This costs some 
                                      * additional memory - for DG
                                      * methods about 1/3 (for double inverses) 
@@ -95,14 +101,22 @@ class PreconditionBlock: public Subscriptor
                                      *
                                      * It is not allowed to call this function
                                      * twice (will produce an error) before
-                                     * a call of #reinit(..)#
+                                     * a call of #clear(..)#
                                      * because at the second time there already
                                      * exist the inverse matrices.
+                                     *
+                                     * After this function is called, the
+                                     * lock on the matrix given through the
+                                     * #use_matrix# function is released,
+                                     * i.e. you may overwrite of delete it.
+                                     * You may want to do this in case
+                                     * you use this matrix to precondition
+                                     * another matrix.
                                      */
     void invert_diagblocks();
 
                                     /**
-                                     * Gives back the size of the blocks.
+                                     * Return the size of the blocks.
                                      */
     unsigned int block_size () const;
 
@@ -191,7 +205,7 @@ class PreconditionBlockSOR : public PreconditionBlock<number, inverse_type>
                                      * Constructor. Takes the damping
                                      * Parameter as 
                                      */
-    PreconditionBlockSOR(number omega=1.);
+    PreconditionBlockSOR (const number omega=1.);
     
                                     /**
                                      * Destructor.
@@ -220,7 +234,10 @@ class PreconditionBlockSOR : public PreconditionBlock<number, inverse_type>
     void operator() (Vector<number2>&, const Vector<number2>&) const;
 
   private:
-    number omega;
+                                    /**
+                                     * Damping parameter.
+                                     */
+    const number omega;
 };
 
 
@@ -229,3 +246,4 @@ class PreconditionBlockSOR : public PreconditionBlock<number, inverse_type>
 /* end of #ifndef __precondition_block_H */
 #endif
 /*------------------------   precondition_block.h     ---------------------------*/
+
index a15886ff343e3f2f7dd84b224bdf26c52b84c7de..bb747393efca41834d99ce73e1c7ef918f850d0d 100644 (file)
@@ -34,8 +34,7 @@ template <typename number, typename inverse_type>
 void PreconditionBlock<number, inverse_type>::use_matrix(
   const SparseMatrix<number> &M)
 {
-  A=0;
-  A=&M;
+  A = &M;
 }
 
 
@@ -62,7 +61,7 @@ void PreconditionBlock<number, inverse_type>::invert_diagblocks()
   Assert (blocksize!=0, ExcBlockSizeNotSet());
   Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
 
-  unsigned int n_cells=M.m()/blocksize;
+  const unsigned int n_cells = M.m()/blocksize;
 
                                   // cell_row, cell_column are the
                                   // numbering of the blocks (cells).
@@ -72,8 +71,21 @@ void PreconditionBlock<number, inverse_type>::invert_diagblocks()
                                   // row, column are the global numbering
                                   // of the unkowns.
 
-  inverse = vector<FullMatrix<inverse_type> > (
-    n_cells, FullMatrix<inverse_type>(blocksize));
+                                  // set the #inverse# array to the right
+                                  // size. we could do it like this:
+                                  // inverse = vector<>(n_cells,FullMatrix<>())
+                                  // but this would involve copying many
+                                  // FullMatrix objects.
+                                  //
+                                  // the following is a neat trick which
+                                  // avoids copying
+  if (true)
+    {
+      vector<FullMatrix<inverse_type> > tmp(n_cells,
+                                           FullMatrix<inverse_type>(blocksize));
+      inverse.swap (tmp);
+    };
+  
   FullMatrix<inverse_type> M_cell(blocksize);
   
   for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
@@ -114,14 +126,15 @@ PreconditionBlockSOR<number, inverse_type>::~PreconditionBlockSOR(){}
 
 template <typename number, typename inverse_type>
 template <typename number2>
-void PreconditionBlockSOR<number, inverse_type>::operator() (Vector<number2> &dst, const Vector<number2> &src) const
+void PreconditionBlockSOR<number, inverse_type>::operator() (Vector<number2> &dst,
+                                                            const Vector<number2> &src) const
 {
   Assert(A!=0, ExcNoMatrixGivenToUse());
   const SparseMatrix<number> &M=*A;
   Assert (M.m() == M.n(), ExcMatrixNotSquare());
   Assert (blocksize!=0, ExcBlockSizeNotSet());
   Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
-  unsigned int n_cells=M.m()/blocksize;
+  const unsigned int n_cells=M.m()/blocksize;
   Assert (inverse.size()==0 || inverse.size()==n_cells,
          ExcWrongNumberOfInverses(inverse.size(), n_cells));
 
index 06ff1a3bb7a849b627c1d86f2ebca01aec6ee7df..88f0b6e956f98c72a17ed53b6910de6c2429d6ec 100644 (file)
@@ -26,9 +26,10 @@ class ostream;
 
 
 /*
-
-   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
-   */
+ * Structure representing the sparsity pattern of a sparse matrix.
+ *
+ * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
+ */
 class SparseMatrixStruct : public Subscriptor
 {
   public:
@@ -337,7 +338,7 @@ class SparseMatrixStruct : public Subscriptor
                                      * readonly.
                                      *
                                      * Though the return value is declared
-                                     * #const#, you shoudl be aware that it
+                                     * #const#, you should be aware that it
                                      * may change if you call any nonconstant
                                      * function of objects which operate on
                                      * it.
@@ -407,12 +408,126 @@ class SparseMatrixStruct : public Subscriptor
     DeclException0 (ExcNotSquare);
     
   private:
+                                    /**
+                                     * Maximum number of rows that can
+                                     * be stored in the #row_start# array.
+                                     * Since reallocation of that array
+                                     * only happens if the present one is
+                                     * too small, but never when the size
+                                     * of this matrix structure shrinks,
+                                     * #max_dim# might be larger than
+                                     * #rows# and in this case #row_start#
+                                     * has more elements than are used.
+                                     */
     unsigned int max_dim;
-    unsigned int rows, cols;
-    unsigned int vec_len, max_vec_len;
+
+                                    /**
+                                     * Number of rows that this sparsity
+                                     * structure shall represent.
+                                     */
+    unsigned int rows;
+
+                                    /**
+                                     * Number of columns that this sparsity
+                                     * structure shall represent.
+                                     */
+    unsigned int cols;
+
+                                    /**
+                                     * Size of the used part of the
+                                     * #colnums# array. Might be lower than
+                                     * #max_vec_len# if the size was reduced
+                                     * somewhen in the past.
+                                     */
+    unsigned int vec_len;
+
+                                    /**
+                                     * Size of the actually allocated array
+                                     * #colnums#. Here, the same applies as
+                                     * for the #rowstart# array, i.e. it
+                                     * may be larger than the actually used
+                                     * part of the array.
+                                     */
+    unsigned int max_vec_len;
+
+                                    /**
+                                     * Maximum number of elements per
+                                     * row. This is set to the value given
+                                     * to the #reinit# function (or to the
+                                     * constructor). Its value is more
+                                     * or less meaningsless after #compress()#
+                                     * has been called.
+                                     */
     unsigned int max_row_len;
-    unsigned int* rowstart;
-    int* colnums;
+
+                                    /**
+                                     * Array which hold for each row which
+                                     * is the first element in #colnums#
+                                     * belonging to that row. Note that
+                                     * the size of the array is one larger
+                                     * than the number of rows, because
+                                     * the last element is used for
+                                     * #row=rows#, i.e. the row past the
+                                     * last used one. The value of
+                                     * #rowstart[rows]# equals the index
+                                     * of the element past the end in
+                                     * #colnums#; this way, we are able to
+                                     * write loops like
+                                     * #for (i=rowstart[k]; i<rowstart[k+1]; ++i)#
+                                     * also for the last row.
+                                     *
+                                     * Note that the actual size of the
+                                     * allocated memory may be larger than
+                                     * the region that is used. The actual
+                                     * number of elements that was allocated
+                                     * is stored in #max_dim#.
+                                     */
+    unsigned int *rowstart;
+
+                                    /**
+                                     * Array of column numbers. In this array,
+                                     * we store for each non-zero element its
+                                     * column number. The column numbers for
+                                     * the elements in row #r# are stored
+                                     * within the index range
+                                     * #rowstart[r]...rowstart[r+1]#. Therefore
+                                     * to find out whether a given element
+                                     * #(r,c)# exists, we have to check
+                                     * whether the column number #c# exists in
+                                     * the abovementioned range within this
+                                     * array. If it exists, say at position
+                                     * #p# within this array, the value of
+                                     * the respective element in the sparse
+                                     * matrix will also be at position #p#
+                                     * of the values array of that class.
+                                     *
+                                     * At the beginning, all elements of
+                                     * this array are set to #-1# indicating
+                                     * invalid (unused) column numbers
+                                     * (however, note that if this object
+                                     * refers to a square matrix, the diagonal
+                                     * elements are preset, see below). Now, if
+                                     * nonzero elements are added, one column
+                                     * number in the row's respective range
+                                     * after the other is set to the column
+                                     * number of the added element. When
+                                     * compress is called, unused elements
+                                     * (indicated by column numbers #-1#)
+                                     * are eliminated by copying the column
+                                     * number of subsequent rows and the
+                                     * column numbers within each row (with
+                                     * the exception of the diagonal element)
+                                     * are sorted, such that finding whether
+                                     * an element exists and determining its
+                                     * position can be done by a binary search.
+                                     *
+                                     * If this object represents a square
+                                     * matrix, the first element in each
+                                     * row always denotes the diagonal
+                                     * element, i.e.
+                                     * #colnums[rowstart[r]]==r#.
+                                     */
+    int *colnums;
 
                                     /**
                                      * Store whether the #compress# function
@@ -427,17 +542,17 @@ class SparseMatrixStruct : public Subscriptor
 
 
 /*
-CLASS
-   SparseMatrix
-
-   @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
-   */
+ * Sparse matrix.
+ *
+ * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
+ */
 template <typename number>
 class SparseMatrix : public Subscriptor
 {
   public:
                                     /**
-                                     * Type of matrix entries.
+                                     * Type of matrix entries. In analogy to
+                                     * the STL container classes.
                                      */
     typedef number value_type;
     
@@ -765,26 +880,42 @@ class SparseMatrix : public Subscriptor
                                      * #dst#.
                                      */
     template <typename somenumber>
-    somenumber residual (Vector<somenumber>& dst, const Vector<somenumber>& x,
-                    const Vector<somenumber>& b) const;
+    somenumber residual (Vector<somenumber>       &dst,
+                        const Vector<somenumber> &x,
+                        const Vector<somenumber> &b) const;
                                     //
     template <typename somenumber>
-    void precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
-                             const number om = 1.) const;
+    void precondition_Jacobi (Vector<somenumber>       &dst,
+                             const Vector<somenumber> &src,
+                             const number              om = 1.) const;
                                     //
     template <typename somenumber>
-    void precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
-                           const number om = 1.) const;
+    void precondition_SSOR (Vector<somenumber>       &dst,
+                           const Vector<somenumber> &src,
+                           const number              om = 1.) const;
                                     //
     template <typename somenumber>
-    void precondition_SOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
-                          const number om = 1.) const;
-                                    //
+    void precondition_SOR (Vector<somenumber>       &dst,
+                          const Vector<somenumber> &src,
+                          const number              om = 1.) const;
+    
+                                    /**
+                                     * Perform an SSOR step in-place, i.e.
+                                     * without copying to a second vector.
+                                     * #omega# is the damping parameter.
+                                     */
     template <typename somenumber>
-    void SSOR (Vector<somenumber>& dst, const number om = 1.) const;
-                                    //
+    void SSOR (Vector<somenumber> &dst,
+              const number        omega = 1.) const;
+
+                                    /**
+                                     * Perform an SOR step in-place, i.e.
+                                     * without copying to a second vector.
+                                     * #omega# is the damping parameter.
+                                     */
     template <typename somenumber>
-    void SOR (Vector<somenumber>& dst, const number om = 1.) const;
+    void SOR (Vector<somenumber> &dst,
+             const number        om = 1.) const;
 
                                     /**
                                      * Return a (constant) reference to the
@@ -792,7 +923,7 @@ class SparseMatrix : public Subscriptor
                                      * matrix.
                                      *
                                      * Though the return value is declared
-                                     * #const#, you shoudl be aware that it
+                                     * #const#, you should be aware that it
                                      * may change if you call any nonconstant
                                      * function of objects which operate on
                                      * it.
@@ -881,8 +1012,36 @@ class SparseMatrix : public Subscriptor
     DeclException0 (ExcInvalidConstructorCall);
     
   private:
+                                    /**
+                                     * Pointer to the sparsity pattern used
+                                     * for this matrix. In order to guarantee
+                                     * that it is not deleted while still in
+                                     * use, we subscribe to it using the
+                                     * #SmartPointer# class.
+                                     */
     SmartPointer<const SparseMatrixStruct> cols;
-    number* val;
+
+                                    /**
+                                     * Array of values for all the nonzero
+                                     * entries. The position within the matrix,
+                                     * i.e. the row and column number for a
+                                     * given entry can only be deduced using
+                                     * the sparsity pattern. The same holds
+                                     * for the more common operation of
+                                     * finding an entry by its coordinates.
+                                     */
+    number *val;
+
+                                    /**
+                                     * Allocated size of #val#. This can
+                                     * be larger than the actually used part
+                                     * if the size of the matrix was
+                                     * reduced somewhen in the past by
+                                     * associating a sparsity pattern
+                                     * with a smaller size to this object
+                                     * somewhen, using the #reinit#
+                                     * function.
+                                     */
     unsigned int max_len;
     
 
index 5e82d24093bb5323b7c566c3fb853278fb51d7b7..4b51edef3f6bb9f510f089f8d9e1b7f158fa46b5 100644 (file)
  */
 template <typename Number>
 class Vector {
-  protected:
-
-                                    /**
-                                     * Dimension. Actual number of components
-                                     * contained in the vector.
-                                     * Get this number by calling #size()#.
-                                     */
-    unsigned int dim;
-
-                                    /**
-                                     * Amount of memory actually reserved for
-                                     * this vector. This number may be greater
-                                     * than #dim# if a #reinit# was called with
-                                     * less memory requirements than the vector
-                                     * needed last time. At present #reinit#
-                                     * does not free memory when the number of
-                                     * needed elements is reduced.
-                                     */
-    unsigned int maxdim;
-
-                                    /**
-                                     * Pointer to the array of components.
-                                     */
-    Number *val;
-
   public:
-
-                                    /**
-                                     * Declare iterator types just like those
-                                     * for the C++ standard library:
-                                     * Data type stored by this container.
-                                     */
-    typedef Number value_type;
-
                                     /**
                                      * Declare standard types used in all
-                                     * containers.
+                                     * containers. These types parallel
+                                     * those in the #C++# standard libraries
+                                     * #vector<...># class.
                                      */
+    typedef Number value_type;
     typedef value_type* pointer;
     typedef const value_type* const_pointer;
     typedef value_type* iterator;
@@ -85,8 +55,9 @@ class Vector {
     Vector ();
     
                                     /**
-                                     *   Copy-Constructor. Dimension set to that of V , <p>
-                                     *                     all components are copied from V
+                                     * Copy-Constructor. Dimension set to
+                                     * that of V, all components are copied
+                                     * from V
                                      */
     Vector (const Vector<Number>& V);
 
@@ -94,7 +65,7 @@ class Vector {
 // note: I disabled this function for the time being, since egcs1.1.2
 // does not respect the "explicit" keyword for template constructors.
 // this leads to unwanted conversions and in some places to automatically
-// generated temporaries, where this is not a good idea    
+// generated temporaries, where this is not a good idea. [WB]
 //                                  /**
 //                                   * Copy constructor taking a vector of
 //                                   * another data type. This will fail if
@@ -115,7 +86,7 @@ class Vector {
     Vector (const unsigned int n);
     
                                     /**
-                                     *         Destructor. Clears memory
+                                     * Destructor. Clears memory
                                      */
     ~Vector ();
 
@@ -129,28 +100,28 @@ class Vector {
     void clear ();
     
                                     /**
-                                     *  U(0-N) = s       . Fill all components
+                                     * $U(0-N) = s$: fill all components.
                                      */
-    Vector<Number>& operator= (const Number s);
+    Vector<Number> & operator= (const Number s);
     
                                     /**
-                                     *  U = V            . Copy all components
+                                     *  $U = V$: copy all components.
                                      */
-    Vector<Number>& operator= (const Vector<Number>& V);
+    Vector<Number> & operator= (const Vector<Number>& V);
 
                                     /**
-                                     * U = V for different types.
+                                     * $U = V$ for different types.
                                      */
     template<typename Number2>
-    Vector<Number>& operator= (const Vector<Number2>& V);
+    Vector<Number> & operator= (const Vector<Number2>& V);
     
                                     /**
-                                     *  U = U * V        . Scalar Produkt
+                                     * $U = U * V$: scalar product.
                                      */
     Number operator* (const Vector<Number>& V) const;
 
                                     /**
-                                     * Return square of the l2-norm.
+                                     * Return square of the $l_2$-norm.
                                      */
     Number norm_sqr () const;
 
@@ -161,13 +132,13 @@ class Vector {
     Number mean_value () const;
 
                                     /**
-                                     * Return the l1-norm of the vector, i.e.
+                                     * Return the $l_1$-norm of the vector, i.e.
                                      * the sum of the absolute values.
                                      */
     Number l1_norm () const;
 
                                     /**
-                                     * Return the l2-norm of the vector, i.e.
+                                     * Return the $l_2$-norm of the vector, i.e.
                                      * the square root of the sum of the
                                      * squares of the elements.
                                      */
@@ -175,7 +146,8 @@ class Vector {
 
                                     /**
                                      * Return the maximum absolute value of the
-                                     * elements of this vector.
+                                     * elements of this vector, which is the
+                                     * $l_\infty$-norm of a vector.
                                      */
     Number linfty_norm () const;
     
@@ -197,7 +169,8 @@ class Vector {
                                      * On #fast==false#, the vector is filled by
                                      * zeros.
                                      */ 
-    void reinit (const unsigned int N, const bool fast=false);
+    void reinit (const unsigned int N,
+                const bool         fast=false);
     
                                     /**
                                      * Change the dimension to that of the
@@ -208,7 +181,8 @@ class Vector {
                                      * this function is the same as calling
                                      * #reinit (V.size(), fast)#.
                                      */
-    void reinit (const Vector<Number>& V, const bool fast=false);
+    void reinit (const Vector<Number> &V,
+                const bool            fast=false);
     
                                     /**
                                      * Return dimension of the vector. This
@@ -263,14 +237,13 @@ class Vector {
                                      */
                                     //@{
                                     /**
-                                     *  Access Components. returns U(i) , 
-                                     *             INLINE
+                                     * Access components, returns U(i).
                                      */
     Number operator() (const unsigned int i) const;
     
                                     /**
-                                     *  Access Components. returns U(i) , 
-                                     *             INLINE
+                                     * Access components, returns U(i)
+                                     * as a writeable reference.
                                      */
     Number& operator() (const unsigned int i);
                                     //@}
@@ -356,7 +329,7 @@ class Vector {
     void scale (const Number factor);
     
                                     /**
-                                     *  U=a*V. Replacing
+                                     *  U=a*V. Replacing.
                                      */
     void equ (const Number a, const Vector<Number>& V);
     
@@ -382,7 +355,8 @@ class Vector {
                                      * is undefined. No attempt is made to
                                      * catch such situations.
                                      */
-    void ratio (const Vector<Number> &a, const Vector<Number> &b);
+    void ratio (const Vector<Number> &a,
+               const Vector<Number> &b);
                                     //@}
     
     
@@ -465,6 +439,31 @@ class Vector {
                                      * Exception
                                      */
     DeclException0 (ExcIO);
+
+  protected:
+
+                                    /**
+                                     * Dimension. Actual number of components
+                                     * contained in the vector.
+                                     * Get this number by calling #size()#.
+                                     */
+    unsigned int dim;
+
+                                    /**
+                                     * Amount of memory actually reserved for
+                                     * this vector. This number may be greater
+                                     * than #dim# if a #reinit# was called with
+                                     * less memory requirements than the vector
+                                     * needed last time. At present #reinit#
+                                     * does not free memory when the number of
+                                     * needed elements is reduced.
+                                     */
+    unsigned int maxdim;
+
+                                    /**
+                                     * Pointer to the array of components.
+                                     */
+    Number *val;
 };
 
 
index 7ccb5028153cbf2eec103848f47119b7d2c3ff57..f6e989c5306b34925851472e1efa92f651abd23e 100644 (file)
@@ -7,6 +7,8 @@
 
 #include <base/subscriptor.h>
 
+
+
 /**
  * Memory management for vectors. This class is used by all
  * iterative methods to allocate space for auxilliary
index d59ead60edbdfbceeaf661a2cce7bcb8e1e5d72a..72204a29f6efe06e0b84d71ca332aef9b3c30307 100644 (file)
@@ -14,7 +14,7 @@ include $D/deal.II/Make.global_options
 
 ../lib/o/%.go :
        @echo ============================ Compiling with debugging information:   $<
-       $(CXX) $(CXXFLAGS.g) -c $< -o $@
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
 ../lib/o/%.o :
        @echo ============================ Compiling with optimization:   $<
        @$(CXX) $(CXXFLAGS) -c $< -o $@
@@ -24,9 +24,11 @@ lib: ../lib/liblac.g.a ../lib/liblac.a
 
 
 ../lib/liblac.a: $(forward-declarations) $(o-files)
+       @echo ======================================== Updating library:   $@
        @ar ruv $@ $(o-files)
 
 ../lib/liblac.g.a: $(forward-declarations) $(go-files)
+       @echo ======================================== Updating library:   $@
        @ar ruv $@ $(go-files)
 
 clean:
index 1bfd73193a06ac2bbffe50c075b840db35474359..1027c32c9d52435f8f0ed8e0fb0ac3c84d02b41d 100644 (file)
@@ -226,7 +226,7 @@ SparseMatrixStruct::compress ()
                    row_length = 0;
 
                                   // reserve temporary storage to
-                                  // store the entries of one wor
+                                  // store the entries of one row
   int *tmp_entries = new int[max_row_len];
 
                                   // Traverse all rows
@@ -324,6 +324,16 @@ SparseMatrixStruct::empty () const {
 unsigned int
 SparseMatrixStruct::max_entries_per_row () const 
 {
+                                  // if compress() has not yet been
+                                  // called, we can get the maximum
+                                  // number of elements per row using
+                                  // the stored value
+  if (!compressed)
+    return max_row_len;
+
+                                  // if compress() was called, we
+                                  // use a better algorithm which
+                                  // gives us a sharp bound
   unsigned int m = 0;
   for (unsigned int i=1; i<rows; ++i)
     m = max (m, rowstart[i]-rowstart[i-1]);

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.