From: guido Date: Tue, 18 Nov 2003 14:28:28 +0000 (+0000) Subject: avoid diagonal in quadratic matrices X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ffb6a67ab69578ebdf8df2452751255b473e8735;p=dealii-svn.git avoid diagonal in quadratic matrices git-svn-id: https://svn.dealii.org/trunk@8188 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 8dc4ed8f59..7aedbe5867 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -836,7 +836,9 @@ SparseMatrix::precondition_Jacobi (Vector &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::precondition_SSOR (Vector &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::precondition_SOR (Vector& 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::precondition_TSOR (Vector& 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::SOR (Vector& 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::TSOR (Vector& 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::PSOR (Vector& 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::TPSOR (Vector& 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::SOR_step (Vector &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::TSOR_step (Vector &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::SSOR (Vector& 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(); diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h index e63692be9e..ad95028994 100644 --- a/deal.II/lac/include/lac/sparsity_pattern.h +++ b/deal.II/lac/include/lac/sparsity_pattern.h @@ -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 &row_lengths); + const std::vector &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 &row_lengths); + const std::vector &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 &row_lengths); + const std::vector &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 - void copy_from (const FullMatrix &matrix); + void copy_from (const FullMatrix &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(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 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 // diff --git a/deal.II/lac/source/sparsity_pattern.cc b/deal.II/lac/source/sparsity_pattern.cc index 5fb54c8615..dc0992e623 100644 --- a/deal.II/lac/source/sparsity_pattern.cc +++ b/deal.II/lac/source/sparsity_pattern.cc @@ -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 &row_lengths) + const std::vector &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 &row_lengths) + const std::vector &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 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 &row_lengths) + const std::vector &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 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 -void SparsityPattern::copy_from (const FullMatrix &matrix) +void SparsityPattern::copy_from (const FullMatrix &matrix, + const bool optimize_diag) { // first init with the number of // entries per row @@ -506,7 +523,7 @@ void SparsityPattern::copy_from (const FullMatrix &matrix) for (unsigned int col=0; col (const FullMatrix &); -template void SparsityPattern::copy_from (const FullMatrix &); +template void SparsityPattern::copy_from (const FullMatrix &, bool); +template void SparsityPattern::copy_from (const FullMatrix &, bool);