From: guido Date: Wed, 4 Sep 2002 19:44:40 +0000 (+0000) Subject: SparseMatrixEZ reprogrammed X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=95c4d1552a78d32ca2ca179658ee235b16593cda;p=dealii-svn.git SparseMatrixEZ reprogrammed git-svn-id: https://svn.dealii.org/trunk@6359 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index a7956b0e9a..a3eb58fa1c 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -27,11 +27,16 @@ template class FullMatrix; /** * Sparse matrix without sparsity pattern. * - * Documentation follows. + * Instead of using a pre-assembled sparsity pattern, this matrix + * 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 + * for rows. * * 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 otheraliens will appreciate it. By the way, the suffix makes + * some other aliens will appreciate it. By the way, the suffix makes * sense by pronouncing it the American way. * * @author Guido Kanschat, 2002 @@ -39,6 +44,51 @@ template class FullMatrix; template class SparseMatrixEZ : public Subscriptor { + private: + /** + * The class for storing the + * column number of an entry + * together with its value. + */ + struct Entry + { + /** + * Standard constructor. Sets + * @p{column} to + * @p{invalid_entry}. + */ + Entry(); + + /** + * Constructor. Fills column + * and value. + */ + Entry(unsigned int column, + const number& value); + + /** + * The column number. + */ + unsigned int column; + /** + * The value there. + */ + number value; + /** + * Comparison operator for finding. + */ + bool operator==(const Entry&) const; + + /** + * Less than operator for sorting. + */ + bool operator < (const Entry&) const; + /** + * Non-existent column number. + */ + static const unsigned int invalid_entry = static_cast(-1); + }; + public: /** * Type of matrix entries. In analogy to @@ -70,10 +120,12 @@ class SparseMatrixEZ : public Subscriptor /** * Constructor. Generates a * matrix of the given size, - * ready to be filled. + * ready to be filled. The optional parameters */ explicit SparseMatrixEZ (unsigned int n_rows, - unsigned int n_columns = n_rows); + unsigned int n_columns = n_rows, + unsigned int default_row_length = Entry::invalid_entry, + unsigned int default_increment = Entry::invalid_entry); /** * Destructor. Free all memory, but do not @@ -95,7 +147,9 @@ class SparseMatrixEZ : public Subscriptor * entries at this point. */ virtual void reinit (unsigned int n_rows, - unsigned int n_columns = n_rows); + unsigned int n_columns = n_rows, + unsigned int default_row_length = Entry::invalid_entry, + unsigned int default_increment = Entry::invalid_entry); /** * Release all memory and return @@ -320,26 +374,6 @@ class SparseMatrixEZ : public Subscriptor number operator () (const unsigned int i, const unsigned int j) const; - /** - * This function is mostly like - * @p{operator()} in that it - * returns the value of the - * matrix entry @p{(i,j)}. The only - * difference is that if this - * entry does not exist in the - * sparsity pattern, then instead - * of raising an exception, zero - * is returned. While this may be - * convenient in some cases, note - * that it is simple to write - * algorithms that are slow - * compared to an optimal - * solution, since the sparsity - * of the matrix is not used. - */ - 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 @@ -650,113 +684,21 @@ class SparseMatrixEZ : public Subscriptor << "The iterators denote a range of " << arg1 << " elements, but the given number of rows was " << arg2); + private: + /** - * The class for storing the - * column number of an entry - * together with its value. + * Find an entry. Return a + * zero-pointer if the entry does + * not exist. */ - struct Entry - { - /** - * Standard constructor. Sets - * @p{column} to - * @p{invalid_entry}. - */ - Entry(); - - /** - * Constructor. Fills column - * and value. - */ - Entry(unsigned int column, - const number& value); - - /** - * The column number. - */ - unsigned int column; - /** - * The value there. - */ - number value; - /** - * Comparison operator for finding. - */ - bool operator==(const Entry&) const; - - /** - * Less than operator for sorting. - */ - bool operator < (const Entry&) const; - /** - * Non-existent column number. - */ - static const unsigned int invalid_entry = static_cast(-1); - }; - + const Entry* locate (unsigned int row, + unsigned int col) const; /** - * The class for storing each row. + * Find an entry or generate it. */ - class Row - { - public: - /** - * Set an entry to a value. - */ - void set(unsigned int column, - const number& value); - /** - * Add value to an entry. - */ - void add(unsigned int column, - const number& value); - /* - * Access to value. - */ - number& operator() (unsigned int column); - - /** - * Read-only access to value. - */ - const number& operator() (unsigned int column) const; - - /** - * Number of entries. - */ - unsigned int size() const; - - /** - * Start of entry list. - */ - typename std::vector::iterator begin(); - - /** - * Start of constant entry list. - */ - typename std::vector::const_iterator begin() const; - - /** - * End of entry list. - */ - typename std::vector::iterator end(); - - /** - * End of constant entry list. - */ - typename std::vector::const_iterator end() const; - - - private: - /** - * Actual data storage. This - * vector contains the - * entries of a row ordered - * by column number. - */ - std::vector values; - }; - + Entry* allocate (unsigned int row, + unsigned int col); /** * Version of @p{vmult} which only @@ -826,7 +768,6 @@ class SparseMatrixEZ : public Subscriptor const std::pair interval, somenumber *partial_norm) const; - /** * Number of columns. This is * used to check vector @@ -834,10 +775,21 @@ class SparseMatrixEZ : public Subscriptor */ unsigned int n_columns; + /** + * Start of indices rows. Points + * into the data field. + */ + std::vector row_start; + /** * Data storage. */ - std::vector rows; + std::vector data; + + /** + * Increment when a row grows. + */ + unsigned int increment; // make all other sparse matrices // friends @@ -886,101 +838,126 @@ SparseMatrixEZ::Entry::operator<(const Entry& e) const +//----------------------------------------------------------------------// template inline -const number& -SparseMatrixEZ::Row::operator()(unsigned int column) const +unsigned int SparseMatrixEZ::m () const { - // find entry - // return its value - Assert(false, ExcNotImplemented()); - return values[0].value; -} + return row_start.size() - 1; +}; template inline -number& -SparseMatrixEZ::Row::operator()(unsigned int column) +unsigned int SparseMatrixEZ::n () const { - // find entry - // return its value - Assert(false, ExcNotImplemented()); - return values[0].value; -} + return n_columns; +}; template inline -typename std::vector::Entry>::iterator -SparseMatrixEZ::Row::begin() +const SparseMatrixEZ::Entry* SparseMatrixEZ::locate ( + const unsigned int row, + const unsigned int col) const { - return values.begin(); -} - + Assert (row -inline -unsigned int -SparseMatrixEZ::Row::size() const -{ - return values.size(); + const unsigned int end = row_start[row+1]; + for (unsigned int i=row_start[row];icolumn == col) + return entry; + if (entry->column == Entry::invalid_entry) + return 0; + } + return 0; } -template -inline -typename std::vector::Entry>::const_iterator -SparseMatrixEZ::Row::begin() const -{ - return values.begin(); -} - template inline -typename std::vector::Entry>::iterator -SparseMatrixEZ::Row::end() +SparseMatrixEZ::Entry* SparseMatrixEZ::allocate ( + const unsigned int row, + const unsigned int col) { - return values.end(); -} - + Assert (row -inline -typename std::vector::Entry>::const_iterator -SparseMatrixEZ::Row::end() const -{ - return values.end(); + const unsigned int end = row_start[row+1]; + for (unsigned int i=row_start[row];icolumn == col) + return entry; + // entry does not exist, + // create it + if (entry->column > col) + { + // Save original entry + Entry temp = data[i]; + // Insert new entry here to + // make sure all entries + // are ordered by column + // index + entry->column = col; + entry->value = 0; + + // Move all entries in this + // row up by one + for (unsigned int j = i+1;j -inline -unsigned int SparseMatrixEZ::m () const -{ - return rows.size(); -}; - - -template -inline -unsigned int SparseMatrixEZ::n () const -{ - return n_columns; -}; - template inline void SparseMatrixEZ::set (const unsigned int i, - const unsigned int j, - const number value) + const unsigned int j, + const number value) { Assert (ivalue = value; }; @@ -988,39 +965,15 @@ void SparseMatrixEZ::set (const unsigned int i, template inline void SparseMatrixEZ::add (const unsigned int i, - const unsigned int j, - const number value) + const unsigned int j, + const number value) { Assert (ivalue += value; }; - -template -inline -number SparseMatrixEZ::diag_element (const unsigned int i) const -{ - Assert (i -inline -number & SparseMatrixEZ::diag_element (const unsigned int i) -{ - Assert (i -template -void -SparseMatrixEZ::Row::set(unsigned int column, - const number& value) -{ - // Store end of vector - const typename std::vector::iterator end_col = end(); - - // Create Entry for inserting - const Entry e(column, value); - - // Find position for inserting - // should return first Entry with - // higher column number. - typename std::vector::iterator col = lower_bound(begin(), end_col, e); - - // Insert Entry if column did not exist. - // Edit existing entry otherwise. - if (col==end_col) - values.push_back(e); - else - if (col->column == column) - col->value = value; - else - values.insert(col, e); -} - - -template -void -SparseMatrixEZ::Row::add(unsigned int column, - const number& value) -{ - const typename std::vector::const_iterator end_col = end(); - - // Create Entry for inserting - const typename SparseMatrixEZ::Entry e(column, value); - - // Find position for inserting - // should return first Entry with - // higher column number. - typename std::vector::iterator col = lower_bound(begin(), end, e); - - // Insert Entry if column did not exist. - // Edit existing entry otherwise. - if (col==end_col) - values.push_back(Entry(column, value)); - else - if (col->column == column) - col->values += value; - else - values.insert(col, Entry(column, value)); -} - - //----------------------------------------------------------------------// template @@ -75,9 +20,11 @@ SparseMatrixEZ::SparseMatrixEZ(const SparseMatrixEZ&) template SparseMatrixEZ::SparseMatrixEZ(unsigned int n_rows, - unsigned int n_cols) + unsigned int n_cols, + unsigned int default_row_length, + unsigned int default_increment) { - reinit(n_rows, n_cols); + reinit(n_rows, n_cols, default_row_length, default_increment); } @@ -98,10 +45,25 @@ SparseMatrixEZ::operator= (const SparseMatrixEZ&) template void SparseMatrixEZ::reinit(unsigned int n_rows, - unsigned int n_cols) + unsigned int n_cols, + unsigned int default_row_length, + unsigned int default_increment) { + if (default_row_length == Entry::invalid_entry) + default_row_length = 5; + if (default_increment == Entry::invalid_entry) + default_increment = 4; + if (default_increment == 0) + default_increment = 4; + increment = default_increment; + n_columns = n_cols; - rows.resize(n_rows); + row_start.resize(n_rows+1); + data.reserve(default_row_length * n_rows + n_rows * increment); + data.resize(default_row_length * n_rows); + + for (unsigned int i=0;i<=n_rows;++i) + row_start[i] = i * default_row_length; } @@ -110,7 +72,8 @@ void SparseMatrixEZ::clear() { n_columns = 0; - rows.resize(0); + row_start.resize(0); + data.resize(0); } @@ -118,7 +81,7 @@ template bool SparseMatrixEZ::empty() const { - return ((n_columns == 0) && (rows.size()==0)); + return ((n_columns == 0) && (row_start.size()==0)); } @@ -131,16 +94,17 @@ SparseMatrixEZ::vmult (Vector& dst, Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); - typename std::vector::const_iterator row = rows.begin(); - const typename std::vector::const_iterator end_row = rows.end(); - for (unsigned int i=0; row != end_row; ++i, ++row) + const unsigned int end_row = row_start.size() - 1; + for (unsigned int i = 0; i < end_row;++i) { - typename std::vector::const_iterator col = row->begin(); - const typename std::vector::const_iterator end_col = row->end(); - + unsigned int index = row_start[i]; + unsigned int end = row_start[i+1]; double s = 0.; - for (;col != end_col; ++col) - s += col->value * src(col->column); + for (;index != end && data[index].column != Entry::invalid_entry;++index) + { + const Entry& entry = data[index]; + s += entry.value * src(entry.column); + } dst(i) = s; } } @@ -166,16 +130,17 @@ SparseMatrixEZ::vmult_add (Vector& dst, Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); - typename std::vector::const_iterator row = rows.begin(); - const typename std::vector::const_iterator end_row = rows.end(); - for (unsigned int i=0; row != end_row; ++i, ++row) + const unsigned int end_row = row_start.size() - 1; + for (unsigned int i = 0; i < end_row;++i) { - typename std::vector::const_iterator col = row->begin(); - const typename std::vector::const_iterator end_col = row->end(); - + unsigned int index = row_start[i]; + unsigned int end = row_start[i+1]; double s = 0.; - for (;col != end_col; ++col) - s += col->value * src(col->column); + for (;index != end && data[index].column != Entry::invalid_entry;++index) + { + const Entry& entry = data[index]; + s += entry.value * src(entry.column); + } dst(i) += s; } } @@ -189,15 +154,16 @@ SparseMatrixEZ::Tvmult_add (Vector& dst, Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size())); Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size())); - typename std::vector::const_iterator row = rows.begin(); - const typename std::vector::const_iterator end_row = rows.end(); - for (unsigned int i=0; row != end_row; ++i, ++row) + const unsigned int end_row = row_start.size() - 1; + for (unsigned int i = 0; i < end_row;++i) { - typename std::vector::const_iterator col = row->begin(); - const typename std::vector::const_iterator end_col = row->end(); - - for (;col != end_col; ++col) - dst(col->column) += col->value * src(i); + unsigned int index = row_start[i]; + unsigned int end = row_start[i+1]; + for (;index != end && data[index].column != Entry::invalid_entry;++index) + { + const Entry& entry = data[index]; + dst(entry.column) += entry.value * src(i); + } } } @@ -205,12 +171,9 @@ template unsigned int SparseMatrixEZ::memory_consumption() const { - unsigned int result = sizeof (*this) - + sizeof(Row) * sizeof (rows); - - for (typename std::vector::const_iterator r = rows.begin(); - r != rows.end(); ++r) - result += r->size() * sizeof(Entry); - + unsigned int result = + sizeof (*this) + + sizeof(unsigned int) * row_start.capacity(); + + sizeof(SparseMatrixEZ::Entry) * data.capacity(); return result; } diff --git a/tests/lac/sparse_matrices.cc b/tests/lac/sparse_matrices.cc index 3400b898a4..df6e0b4b08 100644 --- a/tests/lac/sparse_matrices.cc +++ b/tests/lac/sparse_matrices.cc @@ -37,7 +37,7 @@ check_vmult_quadratic(const MATRIX& A, Vector f(A.m()); GrowingVectorMemory<> mem; - SolverControl control(50, 1.e-3, false); + SolverControl control(20, 1.e-3, false); SolverRichardson<> rich(control, mem, .01); PreconditionIdentity prec; @@ -73,36 +73,44 @@ int main() // and benchmark #ifdef DEBUG deallog.depth_console(0); - const unsigned int size = 100; + const unsigned int size = 10; #else + deallog.depth_console(1000); deallog.log_execution_time(true); deallog.log_time_differences(true); - const unsigned int size = 1000; + const unsigned int size = 500; #endif FDMatrix testproblem (size, size); unsigned int dim = (size-1)*(size-1); - + + deallog << "Structure" << std::endl; SparsityPattern structure(dim, dim, 5); testproblem.five_point_structure(structure); structure.compress(); SparseMatrix A(structure); + deallog << "Assemble" << std::endl; testproblem.five_point(A); check_vmult_quadratic(A, "5-SparseMatrix"); - SparseMatrixEZ E(dim,dim); + SparseMatrixEZ E(dim,dim,5,1); + deallog << "Assemble" << std::endl; testproblem.five_point(E); check_vmult_quadratic(E, "5-SparseMatrixEZ"); A.clear(); + deallog << "Structure" << std::endl; structure.reinit(dim, dim, 9); testproblem.nine_point_structure(structure); structure.compress(); A.reinit(structure); + deallog << "Assemble" << std::endl; + testproblem.nine_point(A); check_vmult_quadratic(A, "9-SparseMatrix"); E.clear(); - E.reinit(dim,dim); + E.reinit(dim,dim,9,2); + deallog << "Assemble" << std::endl; testproblem.nine_point(E); check_vmult_quadratic(E, "9-SparseMatrixEZ"); diff --git a/tests/lac/sparse_matrices.checked b/tests/lac/sparse_matrices.checked index 2b52c9596b..de3ecad61d 100644 --- a/tests/lac/sparse_matrices.checked +++ b/tests/lac/sparse_matrices.checked @@ -1,29 +1,35 @@ +DEAL::Structure +DEAL::Assemble DEAL:5-SparseMatrix:Richardson::Starting -DEAL:5-SparseMatrix:Richardson::Failure step 50 +DEAL:5-SparseMatrix:Richardson::Failure step 20 DEAL:5-SparseMatrix::Transpose DEAL:5-SparseMatrix:RichardsonT::Starting -DEAL:5-SparseMatrix:RichardsonT::Failure step 50 +DEAL:5-SparseMatrix:RichardsonT::Failure step 20 DEAL::GrowingVectorMemory:Overall allocated vectors: 4 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2 +DEAL::Assemble DEAL:5-SparseMatrixEZ:Richardson::Starting -DEAL:5-SparseMatrixEZ:Richardson::Failure step 50 +DEAL:5-SparseMatrixEZ:Richardson::Failure step 20 DEAL:5-SparseMatrixEZ::Transpose DEAL:5-SparseMatrixEZ:RichardsonT::Starting -DEAL:5-SparseMatrixEZ:RichardsonT::Failure step 50 +DEAL:5-SparseMatrixEZ:RichardsonT::Failure step 20 DEAL::GrowingVectorMemory:Overall allocated vectors: 4 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2 +DEAL::Structure +DEAL::Assemble DEAL:9-SparseMatrix:Richardson::Starting -DEAL:9-SparseMatrix:Richardson::Failure step 50 +DEAL:9-SparseMatrix:Richardson::Failure step 20 DEAL:9-SparseMatrix::Transpose DEAL:9-SparseMatrix:RichardsonT::Starting -DEAL:9-SparseMatrix:RichardsonT::Failure step 50 +DEAL:9-SparseMatrix:RichardsonT::Failure step 20 DEAL::GrowingVectorMemory:Overall allocated vectors: 4 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2 +DEAL::Assemble DEAL:9-SparseMatrixEZ:Richardson::Starting -DEAL:9-SparseMatrixEZ:Richardson::Failure step 50 +DEAL:9-SparseMatrixEZ:Richardson::Failure step 20 DEAL:9-SparseMatrixEZ::Transpose DEAL:9-SparseMatrixEZ:RichardsonT::Starting -DEAL:9-SparseMatrixEZ:RichardsonT::Failure step 50 +DEAL:9-SparseMatrixEZ:RichardsonT::Failure step 20 DEAL::GrowingVectorMemory:Overall allocated vectors: 4 DEAL::GrowingVectorMemory:Maximum allocated vectors: 2 diff --git a/tests/lac/testmatrix.cc b/tests/lac/testmatrix.cc index 1f922f2e54..5b13362163 100644 --- a/tests/lac/testmatrix.cc +++ b/tests/lac/testmatrix.cc @@ -178,7 +178,6 @@ FDMatrix::five_point(MATRIX& A) const { // Number of the row to be entered unsigned int row = j+(nx-1)*i; - A.set(row, row, 4.); if (j>0) {