]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
avoid diagonal in quadratic matrices
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 18 Nov 2003 14:28:28 +0000 (14:28 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 18 Nov 2003 14:28:28 +0000 (14:28 +0000)
git-svn-id: https://svn.dealii.org/trunk@8188 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix.templates.h
deal.II/lac/include/lac/sparsity_pattern.h
deal.II/lac/source/sparsity_pattern.cc

index 8dc4ed8f59d3cdb5d79002a25d2bfae1530d34f0..7aedbe586785664c000255d56cd88f455b58d810 100644 (file)
@@ -836,7 +836,9 @@ SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>       &dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
   Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
 
@@ -879,7 +881,9 @@ SparseMatrix<number>::precondition_SSOR (Vector<somenumber>       &dst,
                                   // which is much clearer...
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
   Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
 
@@ -945,7 +949,9 @@ SparseMatrix<number>::precondition_SOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
 
   dst = src;
   SOR(dst,om);
@@ -961,7 +967,9 @@ SparseMatrix<number>::precondition_TSOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
 
   dst = src;
   TSOR(dst,om);
@@ -976,7 +984,9 @@ SparseMatrix<number>::SOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
 
   for (unsigned int row=0; row<m(); ++row)
@@ -1002,7 +1012,9 @@ SparseMatrix<number>::TSOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
 
   for (unsigned int row=m(); row!=0;)
@@ -1028,7 +1040,9 @@ SparseMatrix<number>::PSOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
          ExcDimensionMismatch(m(), permutation.size()));
@@ -1064,7 +1078,9 @@ SparseMatrix<number>::TPSOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
          ExcDimensionMismatch(m(), permutation.size()));
@@ -1098,7 +1114,9 @@ SparseMatrix<number>::SOR_step (Vector<somenumber> &v,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
   Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
 
@@ -1124,7 +1142,9 @@ SparseMatrix<number>::TSOR_step (Vector<somenumber> &v,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
   Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
 
@@ -1162,7 +1182,9 @@ SparseMatrix<number>::SSOR (Vector<somenumber>& dst,
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
-  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (cols->optimize_diagonal(),
+         typename SparsityPattern::ExcNotQuadratic());
+  
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
 
   const unsigned int  n = dst.size();
index e63692be9ec20e09a3c646f7533eb25b481c425f..ad950289943fbd5a3438a7144cf63f89922c01ab 100644 (file)
@@ -181,10 +181,18 @@ class SparsityPattern : public Subscriptor
                                      * @p{m} rows and @p{n} columns.
                                      * The matrix may contain at most @p{max_per_row}
                                      * nonzero entries per row.
+                                     *
+                                     * If the matrix is quadratic,
+                                     * then the last parameter
+                                     * controls optimized storage of
+                                     * diagonal elements for
+                                     * relaxation methods of
+                                     * SparseMatrix.
                                      */
     SparsityPattern (const unsigned int m,
                     const unsigned int n,
-                    const unsigned int max_per_row);
+                    const unsigned int max_per_row,
+                    const bool optimize_diagonal = true);
 
                                     /**
                                      * Initialize a rectangular
@@ -196,18 +204,26 @@ class SparsityPattern : public Subscriptor
                                      */
     SparsityPattern (const unsigned int               m,
                     const unsigned int               n,
-                    const std::vector<unsigned int> &row_lengths);
+                    const std::vector<unsigned int> &row_lengths,
+                    const bool optimize_diagonal = true);
     
                                     /**
-                                     * Initialize a square matrix of dimension
+                                     * Initialize a quadratic matrix of dimension
                                      * @p{n} with at most @p{max_per_row}
                                      * nonzero entries per row.
+                                     *
+                                     * This constructor automatically
+                                     * enables optimized storage of
+                                     * diagonal elements. To avoid
+                                     * this, use the constructor
+                                     * taking row and column numbers
+                                     * separately.
                                      */
     SparsityPattern (const unsigned int n,
                     const unsigned int max_per_row);
 
                                     /**
-                                     * Initialize a square
+                                     * Initialize a quadratic
                                      * matrix with @p{m} rows and @p{m}
                                      * columns.  The maximal number
                                      * of nonzero entries for each
@@ -215,7 +231,8 @@ class SparsityPattern : public Subscriptor
                                      * @p{row_lengths} array.
                                      */
     SparsityPattern (const unsigned int               m,
-                    const std::vector<unsigned int> &row_lengths);
+                    const std::vector<unsigned int> &row_lengths,
+                    const bool optimize_diagonal = true);
 
                                     /**
                                      * Make a copy with extra off-diagonals.
@@ -246,7 +263,7 @@ class SparsityPattern : public Subscriptor
                                      * than their being on side-diagonals.
                                      *
                                      * This function requires that @p{original}
-                                     * refer to a square matrix structure.
+                                     * refer to a quadratic matrix structure.
                                      * It shall be compressed. The matrix 
                                      * structure is not compressed
                                      * after this function finishes.
@@ -283,7 +300,8 @@ class SparsityPattern : public Subscriptor
                                      */
     void reinit (const unsigned int m,
                 const unsigned int n,
-                const unsigned int max_per_row);
+                const unsigned int max_per_row,
+                const bool optimize_diagonal = true);
 
                                     /**
                                      * Reallocate memory for a matrix
@@ -301,10 +319,19 @@ class SparsityPattern : public Subscriptor
                                      * size extends the old one. This is done
                                      * to save time and to avoid fragmentation
                                      * of the heap.
+                                     *
+                                     * IF the number of rows equals
+                                     * the number of columns and the
+                                     * last parameter is true,
+                                     * diagonals elements are stored
+                                     * first in each row to allow
+                                     * optimized access in relaxation
+                                     * methods of SparseMatrix.
                                      */
     void reinit (const unsigned int               m,
                 const unsigned int               n,
-                const std::vector<unsigned int> &row_lengths);
+                const std::vector<unsigned int> &row_lengths,
+                const bool optimize_diagonal = true);
     
                                     /**
                                      * This function compresses the sparsity
@@ -314,7 +341,7 @@ class SparsityPattern : public Subscriptor
                                      * ones to allow faster access by usage
                                      * of binary search algorithms. A special
                                      * sorting scheme is used for the diagonal
-                                     * entry of square matrices, which is
+                                     * entry of quadratic matrices, which is
                                      * always the first entry of each row.
                                      *
                                      * The memory which is no more
@@ -470,7 +497,8 @@ class SparsityPattern : public Subscriptor
     void copy_from (const unsigned int    n_rows,
                    const unsigned int    n_cols,
                    const ForwardIterator begin,
-                   const ForwardIterator end);
+                   const ForwardIterator end,
+                   const bool optimize_diagonal = true);
 
                                     /**
                                      * Copy data from an object of
@@ -481,7 +509,8 @@ class SparsityPattern : public Subscriptor
                                      * sparsity pattern is in
                                      * compressed mode afterwards.
                                      */
-    void copy_from (const CompressedSparsityPattern &csp);
+    void copy_from (const CompressedSparsityPattern &csp,
+                   const bool optimize_diagonal = true);
 
                                     /**
                                      * Take a full matrix and use its
@@ -495,7 +524,8 @@ class SparsityPattern : public Subscriptor
                                      * compressed mode afterwards.
                                      */
     template <typename number>
-    void copy_from (const FullMatrix<number> &matrix);
+    void copy_from (const FullMatrix<number> &matrix,
+                   const bool optimize_diagonal = true);
     
                                     /**
                                      * Return whether the object is empty. It
@@ -586,7 +616,7 @@ class SparsityPattern : public Subscriptor
                                      * This function throws an
                                      * exception if the sparsity
                                      * pattern does not represent a
-                                     * square matrix.
+                                     * quadratic matrix.
                                      */
     void symmetrize ();
     
@@ -646,10 +676,11 @@ class SparsityPattern : public Subscriptor
                                      * Access to column number field.
                                      * Return the column number of
                                      * the @p{index}th entry in
-                                     * @p{row}. Note that the if the
-                                     * matrix is square, then the
-                                     * first element in each row is
-                                     * the diagonal element,
+                                     * @p{row}. Note that the if
+                                     * diagonal elements are
+                                     * optimized, the first element
+                                     * in each row is the diagonal
+                                     * element,
                                      * i.e. @p{column_number(row,0)==row}.
                                      *
                                      * If the sparsity pattern is
@@ -692,6 +723,15 @@ class SparsityPattern : public Subscriptor
                                      * compressed or not.
                                      */
     bool is_compressed () const;
+
+                                    /**
+                                     * Determine whether the matrix
+                                     * uses special convention for
+                                     * quadratic matrices,
+                                     * i.e. diagonal element first in
+                                     * row.
+                                     */
+    bool optimize_diagonal () const;
     
                                     /**
                                      * This is kind of an expert mode. Get 
@@ -842,7 +882,11 @@ class SparsityPattern : public Subscriptor
                                     /**
                                      * Exception
                                      */
-    DeclException0 (ExcNotSquare);
+    DeclException0 (ExcNotQuadratic);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcDiagonalNotOptimized);
                                     /**
                                      * Exception
                                      */
@@ -943,31 +987,31 @@ class SparsityPattern : public Subscriptor
                                      * matrix will also be at position @p{p}
                                      * of the values array of that class.
                                      *
-                                     * At the beginning, all elements of
-                                     * this array are set to @p{-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 @p{-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.
-                                     * @p{colnums[rowstart[r]]==r}.
+                                     * At the beginning, all elements
+                                     * of this array are set to
+                                     * @p{-1} indicating invalid
+                                     * (unused) column numbers
+                                     * (diagonal elements are preset
+                                     * if optimized storage is
+                                     * requested, though). 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
+                                     * @p{-1}) are eliminated by
+                                     * copying the column number of
+                                     * subsequent rows and the column
+                                     * numbers within each row (with
+                                     * possible 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.
                                      */
     unsigned int *colnums;
 
@@ -977,6 +1021,12 @@ class SparsityPattern : public Subscriptor
                                      */
     bool compressed;
 
+                                    /**
+                                     * Is special treatment of
+                                     * diagonals enabled?
+                                     */
+    bool diagonal_optimized;
+    
                                     /**
                                      * Optimized replacement for
                                      * @p{std::lower_bound} for
@@ -1164,6 +1214,14 @@ SparsityPattern::is_compressed () const
 }
 
 
+inline
+bool
+SparsityPattern::optimize_diagonal () const
+{
+  return (rows == cols);
+}
+
+
 inline
 const unsigned int *
 SparsityPattern::get_rowstart_indices () const
@@ -1247,14 +1305,15 @@ void
 SparsityPattern::copy_from (const unsigned int    n_rows,
                            const unsigned int    n_cols,
                            const ForwardIterator begin,
-                           const ForwardIterator end)
+                           const ForwardIterator end,
+                           const bool optimize_diag)
 {
   Assert (static_cast<unsigned int>(std::distance (begin, end)) == n_rows,
          ExcIteratorRange (std::distance (begin, end), n_rows));
   
                                   // first determine row lengths for
                                   // each row. if the matrix is
-                                  // square, then we might have to
+                                  // quadratic, then we might have to
                                   // add an additional entry for the
                                   // diagonal, if that is not yet
                                   // present. as we have to call
@@ -1262,18 +1321,18 @@ SparsityPattern::copy_from (const unsigned int    n_rows,
                                   // bother to check whether that
                                   // diagonal entry is in a certain
                                   // row or not
-  const bool is_square = (n_rows == n_cols);
+  const bool is_square = optimize_diag && (n_rows == n_cols);
   std::vector<unsigned int> row_lengths;
   row_lengths.reserve(n_rows);
   for (ForwardIterator i=begin; i!=end; ++i)
     row_lengths.push_back (std::distance (i->begin(), i->end())
                           +
                           (is_square ? 1 : 0));
-  reinit (n_rows, n_cols, row_lengths);
+  reinit (n_rows, n_cols, row_lengths, is_square);
 
                                   // now enter all the elements into
                                   // the matrix. note that if the
-                                  // matrix is square, then we
+                                  // matrix is quadratic, then we
                                   // already have the diagonal
                                   // element preallocated
                                   //
index 5fb54c8615052128f7e6d96c380ff5b9c0c2af51..dc0992e62344df4b1f7bcf76edb4eb33bcca04ce 100644 (file)
@@ -32,7 +32,9 @@ SparsityPattern::SparsityPattern () :
                max_dim(0),
                max_vec_len(0),
                rowstart(0),
-               colnums(0)
+               colnums(0),
+               compressed(false),
+               diagonal_optimized(false)
 {
   reinit (0,0,0);
 }
@@ -44,40 +46,46 @@ SparsityPattern::SparsityPattern (const SparsityPattern &s) :
                max_dim(0),
                max_vec_len(0),
                rowstart(0),
-               colnums(0)
+               colnums(0),
+               compressed(false),
+               diagonal_optimized(false)
 {
   Assert (s.rowstart == 0, ExcInvalidConstructorCall());
   Assert (s.colnums == 0, ExcInvalidConstructorCall());
   Assert (s.rows == 0, ExcInvalidConstructorCall());
   Assert (s.cols == 0, ExcInvalidConstructorCall());
   
-  reinit (0,0,0);
+  reinit (0,0,0, false);
 }
 
 
 
 SparsityPattern::SparsityPattern (const unsigned int m,
                                  const unsigned int n,
-                                 const unsigned int max_per_row) 
+                                 const unsigned int max_per_row,
+                                 const bool optimize_diag) 
                : max_dim(0),
                  max_vec_len(0),
                  rowstart(0),
-                 colnums(0)
+                 colnums(0),
+                 compressed(false),
+                 diagonal_optimized(false)
 {
-  reinit (m,n,max_per_row);
+  reinit (m,n,max_per_row, optimize_diag);
 }
 
 
 
 SparsityPattern::SparsityPattern (const unsigned int               m,
                                  const unsigned int               n,
-                                 const std::vector<unsigned int> &row_lengths) 
+                                 const std::vector<unsigned int> &row_lengths,
+                                 const bool optimize_diag) 
                : max_dim(0),
                  max_vec_len(0),
                  rowstart(0),
                  colnums(0)
 {
-  reinit (m, n, row_lengths);
+  reinit (m, n, row_lengths, optimize_diag);
 }
 
 
@@ -89,19 +97,20 @@ SparsityPattern::SparsityPattern (const unsigned int n,
                  rowstart(0),
                  colnums(0)
 {
-  reinit (n,n,max_per_row);
+  reinit (n, n, max_per_row, true);
 }
 
 
 
 SparsityPattern::SparsityPattern (const unsigned int               m,
-                                 const std::vector<unsigned int> &row_lengths) 
+                                 const std::vector<unsigned int> &row_lengths,
+                                 const bool optimize_diag) 
                : max_dim(0),
                  max_vec_len(0),
                  rowstart(0),
                  colnums(0)
 {
-  reinit (m, m, row_lengths);
+  reinit (m, m, row_lengths, optimize_diag);
 }
 
 
@@ -114,10 +123,12 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original,
                  rowstart(0),
                  colnums(0)
 {
-  Assert (original.rows==original.cols, ExcNotSquare());
+  Assert (original.rows==original.cols, ExcNotQuadratic());
+  Assert (original.optimize_diagonal(), ExcNotQuadratic());
   Assert (original.is_compressed(), ExcNotCompressed());
   
-  reinit (original.rows, original.cols, max_per_row);
+  reinit (original.rows, original.cols, max_per_row,
+         original.optimize_diagonal());
 
                                   // now copy the entries from
                                   // the other object
@@ -236,12 +247,13 @@ SparsityPattern::operator = (const SparsityPattern &s)
 void
 SparsityPattern::reinit (const unsigned int m,
                         const unsigned int n,
-                        const unsigned int max_per_row)
+                        const unsigned int max_per_row,
+                        const bool optimize_diag)
 {
                                   // simply map this function to the
                                   // other @p{reinit} function
   const std::vector<unsigned int> row_lengths (m, max_per_row);
-  reinit (m, n, row_lengths);
+  reinit (m, n, row_lengths, optimize_diag);
 }
 
 
@@ -249,7 +261,8 @@ SparsityPattern::reinit (const unsigned int m,
 void
 SparsityPattern::reinit (const unsigned int               m,
                         const unsigned int               n,
-                        const std::vector<unsigned int> &row_lengths)
+                        const std::vector<unsigned int> &row_lengths,
+                        const bool optimize_diag)
 {
   Assert (row_lengths.size() == m, ExcInvalidNumber (m));
          
@@ -272,13 +285,14 @@ SparsityPattern::reinit (const unsigned int               m,
     };
 
                                   // first, if the matrix is
-                                  // quadratic, we will have to make
-                                  // sure that each row has at least
-                                  // one entry for the diagonal
-                                  // element. make this more obvious
-                                  // by having a variable which we
-                                  // can query
-  const bool matrix_is_quadratic = (m == n ? true : false);
+                                  // quadratic and special treatment
+                                  // of diagonals is requested, we
+                                  // will have to make sure that each
+                                  // row has at least one entry for
+                                  // the diagonal element. make this
+                                  // more obvious by having a
+                                  // variable which we can query
+  diagonal_optimized = (m == n ? true : false) && optimize_diag;
 
                                   // find out how many entries we
                                   // need in the @p{colnums} array. if
@@ -291,7 +305,7 @@ SparsityPattern::reinit (const unsigned int               m,
                                   // of columns
   unsigned int vec_len = 0;
   for (unsigned int i=0; i<m; ++i)
-    vec_len += std::min((matrix_is_quadratic ?
+    vec_len += std::min((diagonal_optimized ?
                         std::max(row_lengths[i], 1U) :
                         row_lengths[i]),
                        n);
@@ -311,7 +325,7 @@ SparsityPattern::reinit (const unsigned int               m,
                    std::min (*std::max_element(row_lengths.begin(), row_lengths.end()),
                              n));
 
-  if (matrix_is_quadratic && (max_row_length==0) && (m!=0))
+  if (diagonal_optimized && (max_row_length==0) && (m!=0))
     max_row_length = 1;
 
                                   // allocate memory for the rowstart
@@ -336,7 +350,7 @@ SparsityPattern::reinit (const unsigned int               m,
   rowstart[0] = 0;
   for (unsigned int i=1; i<=rows; ++i)
     rowstart[i] = rowstart[i-1] +
-                 (matrix_is_quadratic ?
+                 (diagonal_optimized ?
                   std::max(std::min(row_lengths[i-1],n),1U) :
                   std::min(row_lengths[i-1],n));
   Assert ((rowstart[rows]==vec_len)
@@ -349,10 +363,10 @@ SparsityPattern::reinit (const unsigned int               m,
                                   // use
   std::fill_n (&colnums[0], vec_len, invalid_entry);
 
-                                  // if the matrix is square: let the
-                                  // first entry in each row be the
-                                  // diagonal value
-  if (rows == cols)
+                                  // if diagonal elements are
+                                  // special: let the first entry in
+                                  // each row be the diagonal value
+  if (diagonal_optimized)
     for (unsigned int i=0;i<rows;i++)
       colnums[rowstart[i]] = i;
 
@@ -404,18 +418,18 @@ SparsityPattern::compress ()
                                       // now @p{rowstart} is
                                       // the number of entries in
                                       // this line
-
-                                      // for square matrices, the
-                                      // first entry in each row
-                                      // is the diagonal one. In
-                                      // this case only sort the
-                                      // remaining entries, otherwise
-                                      // sort all
-
+      
+                                      // Sort only beginning at the
+                                      // second entry, if optimized
+                                      // storage of diagonal entries
+                                      // is on.
+      
                                       // if this line is empty or has
                                       // only one entry, don't sort
       if (row_length > 1)
-       std::sort ((rows==cols) ? tmp_entries.begin()+1 : tmp_entries.begin(),
+       std::sort ((diagonal_optimized)
+                  ? tmp_entries.begin()+1
+                  : tmp_entries.begin(),
                   tmp_entries.begin()+row_length);
       
                                       // insert column numbers
@@ -435,7 +449,7 @@ SparsityPattern::compress ()
                                       // the diagonal element
                                       // (i.e. with column
                                       // index==line number)
-      Assert ((rows!=cols) ||
+      Assert ((!diagonal_optimized) ||
              (new_colnums[rowstart[line]] == line),
              ExcInternalError());
                                       // assert that the first entry
@@ -488,16 +502,19 @@ SparsityPattern::compress ()
 
 
 void
-SparsityPattern::copy_from (const CompressedSparsityPattern &csp) 
+SparsityPattern::copy_from (const CompressedSparsityPattern &csp,
+                           const bool optimize_diag) 
 {
   copy_from (csp.n_rows(), csp.n_cols(),
             csp.column_indices.begin(),
-            csp.column_indices.end());
+            csp.column_indices.end(),
+            optimize_diag);
 }
 
 
 template <typename number>
-void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
+void SparsityPattern::copy_from (const FullMatrix<number> &matrix,
+                                const bool optimize_diag)
 {
                                   // first init with the number of
                                   // entries per row
@@ -506,7 +523,7 @@ void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
     for (unsigned int col=0; col<matrix.n(); ++col)
       if (matrix(row,col) != 0)
        ++entries_per_row[row];
-  reinit (matrix.m(), matrix.n(), entries_per_row);
+  reinit (matrix.m(), matrix.n(), entries_per_row, optimize_diag);
 
                                   // now set entries
   for (unsigned int row=0; row<matrix.m(); ++row)
@@ -582,11 +599,12 @@ SparsityPattern::operator () (const unsigned int i,
                                   // something in this line
   if (rowstart[i] == rowstart[i+1])
     return invalid_entry;
-  
-                                  // check first entry separately, since
-                                  // for square matrices this is
-                                  // the diagonal entry
-  if ((i==j) && (rows==cols))
+
+                                  // If special storage of diagonals
+                                  // was requested, we can get the
+                                  // diagonal element faster by this
+                                  // query.
+  if (diagonal_optimized && (i==j))
     return rowstart[i];
 
                                   // all other entries are sorted, so
@@ -601,7 +619,7 @@ SparsityPattern::operator () (const unsigned int i,
                                   // top of this function, so it may
                                   // not be called for noncompressed
                                   // structures.
-  const unsigned int * const sorted_region_start = (rows==cols ?
+  const unsigned int * const sorted_region_start = (diagonal_optimized ?
                                                    &colnums[rowstart[i]+1] :
                                                    &colnums[rowstart[i]]);
   const unsigned int * const p = optimized_lower_bound (sorted_region_start,
@@ -698,7 +716,10 @@ SparsityPattern::symmetrize ()
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());  
   Assert (compressed==false, ExcMatrixIsCompressed());
-  Assert (rows==cols, ExcNotSquare());
+                                  // Note that we only require a
+                                  // quadratic matrix here, no
+                                  // special treatment of diagonals
+  Assert (rows==cols, ExcNotQuadratic());
 
                                   // loop over all elements presently
                                   // in the sparsity pattern and add
@@ -774,7 +795,7 @@ SparsityPattern::bandwidth () const
 }
 
 
-
+//TODO: Write and read flag diagonal_optimized?
 void
 SparsityPattern::block_write (std::ostream &out) const 
 {
@@ -854,5 +875,5 @@ SparsityPattern::memory_consumption () const
 
 
 // explicit instantiations
-template void SparsityPattern::copy_from<float> (const FullMatrix<float> &);
-template void SparsityPattern::copy_from<double> (const FullMatrix<double> &);
+template void SparsityPattern::copy_from<float> (const FullMatrix<float> &, bool);
+template void SparsityPattern::copy_from<double> (const FullMatrix<double> &, bool);

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.