From edbd520bbc8e3fc355661c0fc0726260bfa63640 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 6 Sep 2002 13:36:49 +0000 Subject: [PATCH] SparseMatrixEZ const_iterator git-svn-id: https://svn.dealii.org/trunk@6370 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix_ez.h | 455 ++++++++++++++++++--- tests/lac/sparse_matrices.cc | 70 +++- tests/lac/sparse_matrices.checked | 69 ++++ 3 files changed, 526 insertions(+), 68 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index b74ab7af60..a04a120a9d 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -31,9 +31,50 @@ template class FullMatrix; * builds the pattern on the fly. Filling the matrix may consume more * time as with @p{SparseMatrix}, since large memory movements may be * involved. To help optimizing things, an expected row-length may be - * provided to the constructor, as well as a mininmum size increment + * provided to the constructor, as well as an increment size * for rows. * + * The storage structure: like with the usual sparse matrix, it is + * attempted to store only non-zero elements. these are stored in a + * single data array @p{data}. They are ordered by row and inside each + * row by column number. Each row is described by its starting point + * in the data array and its length. These are stored in the + * @p{row_info} array, together with additional useful information. + * + * Due to the structure, gaps may occur between rows. Whenever a new + * entry must be created, an attempt is made to use the gap in its + * row. If the gap is full, the row must be extended and all + * subsequent rows must be shifted backwards. This is a very expensive + * operation and should be avoided as much as possible. + * + * This is, where the optimization parameters, provided to the + * constructor or to the function @p{reinit} come + * in. @p{default_row_length} is the amount of entries that will be + * allocated for each row on initialization (the actual length of the + * rows is still zero). This means, that @p{default_row_length} + * entries can be added to this row without shifting other rows. If + * less entries are added, the additional memory will be wasted. + * + * If the space for a row is not sufficient, then it is enlarged by + * @p{default_increment} entries. This way, the subsequent rows are + * not shifted by single entries very often. + * + * Finally, the @p{default_reserve} allocates extra space at the end + * of the data array. This space is used whenever a row must be + * enlarged. Since @p{std::vector} doubles the capacity everytime it + * must increase it, this value should allow for all the growth needed. + * + * Suggested settings: @p{default_row_length} should be the length of + * a typical row, for instance the size of the stencil in regular + * parts of the grid. Then, @p{default_increment} may be the expected + * amount of entries added to the row by having one hanging node. This + * way, a good compromise between memory consumption and speed should + * be achieved. @p{default_reserve} should then be an estimate for the + * number of hanging nodes times @p{default_increment}. + * + * If the rows are expected to be filled more or less from first to + * last, using a @p{default_row_length} may not be such a bad idea. + * * The name of this matrix is in reverence to a publication of the * Internal Revenue Service of the United States of America. I hope * some other aliens will appreciate it. By the way, the suffix makes @@ -124,9 +165,124 @@ class SparseMatrixEZ : public Subscriptor static const unsigned short invalid_diagonal = static_cast(-1); }; - - + public: + /** + * Accessor class for iterators + */ + class Accessor + { + public: + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const SparseMatrixEZ*, + unsigned int row, + unsigned short index); + + /** + * Row number of the element + * represented by this + * object. + */ + unsigned int row() const; + + /** + * Index in row of the element + * represented by this + * object. + */ + unsigned short index() const; + + /** + * Column number of the + * element represented by + * this object. + */ + unsigned int column() const; + + /** + * Value of this matrix entry. + */ + number value() const; + + protected: + /** + * The matrix accessed. + */ + const SparseMatrixEZ* matrix; + + /** + * Current row number. + */ + unsigned int a_row; + + /** + * Current index in row. + */ + unsigned short a_index; + }; + + /** + * STL conforming iterator. + */ + class const_iterator : private Accessor + { + public: + /** + * Constructor. + */ + const_iterator(const SparseMatrixEZ*, + unsigned int row, + unsigned short index); + + /** + * Prefix increment. + */ + const_iterator& operator++ (); + + /** + * Postfix increment. + */ + const_iterator& operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor& operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor* operator-> () const; + + /** + * Comparison. True, if + * both iterators point to + * the same matrix + * position. + */ + bool operator == (const const_iterator&) const; + /** + * Inverse of @p{==}. + */ + bool operator != (const const_iterator&) const; + + /** + * Comparison + * operator. Result is true + * if either the first row + * number is smaller or if + * the row numbers are + * equal and the first + * index is smaller. + */ + bool operator < (const const_iterator&) const; + }; + /** * Type of matrix entries. In analogy to * the STL container classes. @@ -237,31 +393,6 @@ class SparseMatrixEZ : public Subscriptor */ unsigned int n () const; - /** - * Return the number of nonzero - * elements of this - * matrix. Actually, it returns - * the number of entries in the - * sparsity pattern; if any of - * the entries should happen to - * be zero, it is counted anyway. - */ -// unsigned int n_nonzero_elements () const; - - /** - * Return the number of actually - * nonzero elements of this - * matrix. - * - * Note, that this function does - * (in contrary to the - * @p{n_nonzero_elements}) NOT - * count all entries of the - * sparsity pattern but only the - * ones that are nonzero. - */ -// unsigned int n_actually_nonzero_elements () const; - /** * Set the element @p{(i,j)} to * @p{value}. Allocates the entry @@ -427,6 +558,14 @@ class SparseMatrixEZ : public Subscriptor number operator () (const unsigned int i, const unsigned int j) const; + /** + * Return the value of the entry + * (i,j). Returns zero for all + * non-existing entries. + */ + number el (const unsigned int i, + const unsigned int j) const; + /** * Return the main diagonal element in * the @p{i}th row. This function throws an @@ -441,14 +580,14 @@ class SparseMatrixEZ : public Subscriptor * involve searching for the * right column number. */ - number diag_element (const unsigned int i) const; +// number diag_element (const unsigned int i) const; /** * Same as above, but return a * writeable reference. You're * sure you know what you do? */ - number & diag_element (const unsigned int i); +// number & diag_element (const unsigned int i); /** * Matrix-vector multiplication: @@ -592,6 +731,42 @@ class SparseMatrixEZ : public Subscriptor void precondition_TSOR (Vector &dst, const Vector &src, const number om = 1.) const; + + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * STL-like iterator with the + * first entry of row @p{r}. + */ + const_iterator begin (unsigned int r) const; + + /** + * Final iterator of row @p{r}. + */ + const_iterator end (unsigned int r) const; + + /** + * Return the number of nonzero + * elements of this + * matrix. + */ +// unsigned int n_nonzero_elements () const; + + /** + * Return the number of actually + * nonzero elements of this + * matrix. + */ +// unsigned int n_actually_nonzero_elements () const; /** * Print the matrix to the given @@ -667,8 +842,15 @@ class SparseMatrixEZ : public Subscriptor */ DeclException0(ExcNoDiagonal); - private: + /** + * Exception + */ + DeclException2 (ExcInvalidIndex, + int, int, + << "The entry with index <" << arg1 << ',' << arg2 + << "> does not exist."); + private: /** * Find an entry. Return a * zero-pointer if the entry does @@ -772,10 +954,9 @@ class SparseMatrixEZ : public Subscriptor * Increment when a row grows. */ unsigned int increment; - - // make all other sparse matrices - // friends - template friend class SparseMatrix; + + friend class Accessor; + friend class const_iterator; }; @@ -801,31 +982,142 @@ SparseMatrixEZ::Entry::Entry() {} -// template -// inline -// bool -// SparseMatrixEZ::Entry::operator==(const Entry& e) const -// { -// return column == e.column; -// } +template +inline +SparseMatrixEZ::RowInfo::RowInfo(unsigned int start) + : + start(start), length(0), diagonal(invalid_diagonal) +{} + + +//----------------------------------------------------------------------// +template +inline +SparseMatrixEZ::Accessor::Accessor ( + const SparseMatrixEZ* matrix, + unsigned int r, + unsigned short i) + : + matrix(matrix), + a_row(r), + a_index(i) +{} + + +template +inline +unsigned int +SparseMatrixEZ::Accessor::row() const +{ + return a_row; +} + + +template +inline +unsigned int +SparseMatrixEZ::Accessor::column() const +{ + return matrix->data[matrix->row_info[a_row].start+a_index].column; +} + + +template +inline +unsigned short +SparseMatrixEZ::Accessor::index() const +{ + return a_index; +} -// template -// inline -// bool -// SparseMatrixEZ::Entry::operator<(const Entry& e) const -// { -// return column < e.column; -// } template inline -SparseMatrixEZ::RowInfo::RowInfo(unsigned int start) +number +SparseMatrixEZ::Accessor::value() const +{ + return matrix->data[matrix->row_info[a_row].start+a_index].value; +} + + +template +inline +SparseMatrixEZ::const_iterator::const_iterator( + const SparseMatrixEZ* matrix, + unsigned int r, + unsigned short i) : - start(start), length(0), diagonal(invalid_diagonal) + Accessor(matrix, r, i) {} +template +inline +SparseMatrixEZ::const_iterator& +SparseMatrixEZ::const_iterator::operator++ () +{ + Assert (a_row < matrix->m(), ExcIteratorPastEnd()); + + ++a_index; + if (a_index >= matrix->row_info[a_row].length) + { + a_index = 0; + a_row++; + } + return *this; +} + + +template +inline +const SparseMatrixEZ::Accessor& +SparseMatrixEZ::const_iterator::operator* () const +{ + return *this; +} + + +template +inline +const SparseMatrixEZ::Accessor* +SparseMatrixEZ::const_iterator::operator-> () const +{ + return this; +} + + +template +inline +bool +SparseMatrixEZ::const_iterator::operator == ( + const const_iterator& other) const +{ + return (row() == other->row() && index() == other->index()); +} + + +template +inline +bool +SparseMatrixEZ::const_iterator::operator != ( + const const_iterator& other) const +{ + return ! (*this == other); +} + + +template +inline +bool +SparseMatrixEZ::const_iterator::operator < ( + const const_iterator& other) const +{ + return (row() < other->row() || + (row() == other->row() && index() < other->index())); +} + + //----------------------------------------------------------------------// template inline @@ -983,5 +1275,66 @@ void SparseMatrixEZ::add (const unsigned int i, }; + +template +inline +number SparseMatrixEZ::el (const unsigned int i, + const unsigned int j) const +{ + Entry* entry = locate(i,j); + if (entry) + return entry->value; + return 0.; +} + + + +template +inline +number SparseMatrixEZ::operator() (const unsigned int i, + const unsigned int j) const +{ + Entry* entry = locate(i,j); + if (entry) + return entry->value; + Assert(false, ExcInvalidEntry(i,j)); + return 0.; +} + + +template +inline +SparseMatrixEZ::const_iterator +SparseMatrixEZ::begin () const +{ + return const_iterator(this, 0, 0); +} + +template +inline +SparseMatrixEZ::const_iterator +SparseMatrixEZ::end () const +{ + return const_iterator(this, m(), 0); +} + +template +inline +SparseMatrixEZ::const_iterator +SparseMatrixEZ::begin (unsigned int r) const +{ + Assert (r +inline +SparseMatrixEZ::const_iterator +SparseMatrixEZ::end (unsigned int r) const +{ + Assert (r +#define PREC_CHECK(solver, method, precond) try \ + { solver.method (A, u, f, precond); } catch (...) {} \ + residuals.push_back(control.last_value()) + template void -check_vmult_quadratic(const MATRIX& A, +check_vmult_quadratic(std::vector& residuals, + const MATRIX& A, const char* prefix) { deallog.push(prefix); @@ -40,7 +45,7 @@ check_vmult_quadratic(const MATRIX& A, SolverControl control(10, 1.e-13, false); SolverRichardson<> rich(control, mem, .01); SolverRichardson<> prich(control, mem, 1.); - PreconditionIdentity prec; + PreconditionIdentity identity; PreconditionJacobi jacobi; jacobi.initialize(A, .5); PreconditionSOR sor; @@ -51,20 +56,41 @@ check_vmult_quadratic(const MATRIX& A, u = 0.; f = 1.; - try { rich.solve(A, u, f, prec); } catch (...) {} - try { prich.solve(A, u, f, jacobi); } catch (...) {} - try { prich.solve(A, u, f, ssor); } catch (...) {} - try { prich.solve(A, u, f, sor); } catch (...) {} - + PREC_CHECK(rich, solve, identity); + PREC_CHECK(prich, solve, jacobi); + PREC_CHECK(prich, solve, ssor); + PREC_CHECK(prich, solve, sor); + u = 0.; deallog << "Transpose" << std::endl; - try { rich.Tsolve(A, u, f, prec); } catch (...) {} - try { prich.Tsolve(A, u, f, jacobi); } catch (...) {} - try { prich.Tsolve(A, u, f, ssor); } catch (...) {} - try { prich.Tsolve(A, u, f, sor); } catch (...) {} + PREC_CHECK(rich, Tsolve, identity); + PREC_CHECK(prich, Tsolve, jacobi); + PREC_CHECK(prich, Tsolve, ssor); + PREC_CHECK(prich, Tsolve, sor); deallog.pop(); } + +template +void +check_iterator (const MATRIX& A) +{ + for (typename MATRIX::const_iterator i = A.begin(); i!= A.end(); ++i) + deallog << '\t' << i->row() + << '\t' << i->column() + << '\t' << i->index() + << '\t' << i->value() + << std::endl; + deallog << "Repeat row 2" << std::endl; + for (typename MATRIX::const_iterator i = A.begin(2); i!= A.end(2); ++i) + deallog << '\t' << i->row() + << '\t' << i->column() + << '\t' << i->index() + << '\t' << i->value() + << std::endl; +} + + int main() { std::ofstream logfile("sparse_matrices.output"); @@ -89,6 +115,9 @@ int main() FDMatrix testproblem (size, size); unsigned int dim = (size-1)*(size-1); + std::vector A_res; + std::vector E_res; + deallog << "Structure" << std::endl; SparsityPattern structure(dim, dim, 5); testproblem.five_point_structure(structure); @@ -96,13 +125,16 @@ int main() SparseMatrix A(structure); deallog << "Assemble" << std::endl; testproblem.five_point(A, true); - check_vmult_quadratic(A, "5-SparseMatrix"); + check_vmult_quadratic(A_res, A, "5-SparseMatrix"); SparseMatrixEZ E(dim,dim,row_length,2); deallog << "Assemble" << std::endl; testproblem.five_point(E, true); - check_vmult_quadratic(E, "5-SparseMatrixEZ"); - + check_vmult_quadratic(E_res, E, "5-SparseMatrixEZ"); +#ifdef DEBUG + check_iterator(E); +#endif + A.clear(); deallog << "Structure" << std::endl; structure.reinit(dim, dim, 9); @@ -111,12 +143,16 @@ int main() A.reinit(structure); deallog << "Assemble" << std::endl; testproblem.nine_point(A); - check_vmult_quadratic(A, "9-SparseMatrix"); + check_vmult_quadratic(A_res, A, "9-SparseMatrix"); E.clear(); E.reinit(dim,dim,row_length,2); deallog << "Assemble" << std::endl; testproblem.nine_point(E); - check_vmult_quadratic(E, "9-SparseMatrixEZ"); - + check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ"); + + for (unsigned int i=0;i 1.e-14) + deallog << "SparseMatrix and SparseMatrixEZ differ!!!" + << std::endl; } diff --git a/tests/lac/sparse_matrices.checked b/tests/lac/sparse_matrices.checked index 15304de6ce..7b9387f15b 100644 --- a/tests/lac/sparse_matrices.checked +++ b/tests/lac/sparse_matrices.checked @@ -40,6 +40,75 @@ DEAL:5-SparseMatrixEZ:RichardsonT::Starting DEAL:5-SparseMatrixEZ:RichardsonT::Failure step 10 DEAL::GrowingVectorMemory:Overall allocated vectors: 16 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2 +DEAL:: 0 0 0 5.00 +DEAL:: 0 1 1 -2.00 +DEAL:: 0 4 2 -1.00 +DEAL:: 1 0 0 -1.00 +DEAL:: 1 1 1 5.00 +DEAL:: 1 2 2 -2.00 +DEAL:: 1 5 3 -1.00 +DEAL:: 2 1 0 -1.00 +DEAL:: 2 2 1 5.00 +DEAL:: 2 3 2 -2.00 +DEAL:: 2 6 3 -1.00 +DEAL:: 3 2 0 -1.00 +DEAL:: 3 3 1 5.00 +DEAL:: 3 7 2 -1.00 +DEAL:: 4 0 0 -1.00 +DEAL:: 4 4 1 5.00 +DEAL:: 4 5 2 -2.00 +DEAL:: 4 8 3 -1.00 +DEAL:: 5 1 0 -1.00 +DEAL:: 5 4 1 -1.00 +DEAL:: 5 5 2 5.00 +DEAL:: 5 6 3 -2.00 +DEAL:: 5 9 4 -1.00 +DEAL:: 6 2 0 -1.00 +DEAL:: 6 5 1 -1.00 +DEAL:: 6 6 2 5.00 +DEAL:: 6 7 3 -2.00 +DEAL:: 6 10 4 -1.00 +DEAL:: 7 3 0 -1.00 +DEAL:: 7 6 1 -1.00 +DEAL:: 7 7 2 5.00 +DEAL:: 7 11 3 -1.00 +DEAL:: 8 4 0 -1.00 +DEAL:: 8 8 1 5.00 +DEAL:: 8 9 2 -2.00 +DEAL:: 8 12 3 -1.00 +DEAL:: 9 5 0 -1.00 +DEAL:: 9 8 1 -1.00 +DEAL:: 9 9 2 5.00 +DEAL:: 9 10 3 -2.00 +DEAL:: 9 13 4 -1.00 +DEAL:: 10 6 0 -1.00 +DEAL:: 10 9 1 -1.00 +DEAL:: 10 10 2 5.00 +DEAL:: 10 11 3 -2.00 +DEAL:: 10 14 4 -1.00 +DEAL:: 11 7 0 -1.00 +DEAL:: 11 10 1 -1.00 +DEAL:: 11 11 2 5.00 +DEAL:: 11 15 3 -1.00 +DEAL:: 12 8 0 -1.00 +DEAL:: 12 12 1 5.00 +DEAL:: 12 13 2 -2.00 +DEAL:: 13 9 0 -1.00 +DEAL:: 13 12 1 -1.00 +DEAL:: 13 13 2 5.00 +DEAL:: 13 14 3 -2.00 +DEAL:: 14 10 0 -1.00 +DEAL:: 14 13 1 -1.00 +DEAL:: 14 14 2 5.00 +DEAL:: 14 15 3 -2.00 +DEAL:: 15 11 0 -1.00 +DEAL:: 15 14 1 -1.00 +DEAL:: 15 15 2 5.00 +DEAL::Repeat row 2 +DEAL:: 2 1 0 -1.00 +DEAL:: 2 2 1 5.00 +DEAL:: 2 3 2 -2.00 +DEAL:: 2 6 3 -1.00 DEAL::Structure DEAL::Assemble DEAL:9-SparseMatrix:Richardson::Starting -- 2.39.5