From: guido Date: Thu, 5 Sep 2002 17:43:24 +0000 (+0000) Subject: preconditioners for SparseMatrixEZ X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c39e7e0b5450b42b50ee7c2474b3fb43cdef3c16;p=dealii-svn.git preconditioners for SparseMatrixEZ git-svn-id: https://svn.dealii.org/trunk@6363 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 a3eb58fa1c..b74ab7af60 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -55,7 +55,7 @@ class SparseMatrixEZ : public Subscriptor /** * Standard constructor. Sets * @p{column} to - * @p{invalid_entry}. + * @p{invalid}. */ Entry(); @@ -77,18 +77,55 @@ class SparseMatrixEZ : public Subscriptor /** * Comparison operator for finding. */ - bool operator==(const Entry&) const; +// bool operator==(const Entry&) const; /** * Less than operator for sorting. */ - bool operator < (const Entry&) const; +// bool operator < (const Entry&) const; /** * Non-existent column number. */ - static const unsigned int invalid_entry = static_cast(-1); + static const unsigned int invalid = static_cast(-1); }; + /** + * Structure for storing + * information on a matrix + * row. One object for each row + * will be stored in the matrix. + */ + struct RowInfo + { + /** + * Constructor. + */ + RowInfo (unsigned int start = Entry::invalid); + + /** + * Index of first entry of + * the row in the data field. + */ + unsigned int start; + /** + * Number of entries in this + * row. + */ + unsigned short length; + /** + * Position of the diagonal + * element relative tor the + * start index. + */ + unsigned short diagonal; + /** + * Value for non-existing diagonal. + */ + static const unsigned short + invalid_diagonal = static_cast(-1); + }; + + public: /** * Type of matrix entries. In analogy to @@ -120,12 +157,20 @@ class SparseMatrixEZ : public Subscriptor /** * Constructor. Generates a * matrix of the given size, - * ready to be filled. The optional parameters + * ready to be filled. The + * optional parameters + * @p{default_row_length} and + * @p{default_increment} allow + * for preallocating + * memory. Providing these + * properly is essential for an + * efficient assembling of the + * matrix. */ explicit SparseMatrixEZ (unsigned int n_rows, unsigned int n_columns = n_rows, - unsigned int default_row_length = Entry::invalid_entry, - unsigned int default_increment = Entry::invalid_entry); + unsigned int default_row_length = Entry::invalid, + unsigned int default_increment = Entry::invalid); /** * Destructor. Free all memory, but do not @@ -144,12 +189,20 @@ class SparseMatrixEZ : public Subscriptor * Reinitialize the sparse matrix * to the dimensions provided. * The matrix will have no - * entries at this point. + * entries at this point. The + * optional parameters + * @p{default_row_length} and + * @p{default_increment} allow + * for preallocating + * memory. Providing these + * properly is essential for an + * efficient assembling of the + * matrix. */ virtual void reinit (unsigned int n_rows, unsigned int n_columns = n_rows, - unsigned int default_row_length = Entry::invalid_entry, - unsigned int default_increment = Entry::invalid_entry); + unsigned int default_row_length = Entry::invalid, + unsigned int default_increment = Entry::invalid); /** * Release all memory and return @@ -540,67 +593,6 @@ class SparseMatrixEZ : public Subscriptor const Vector &src, const number om = 1.) const; - /** - * Perform SSOR preconditioning - * in-place. Apply the - * preconditioner matrix without - * copying to a second vector. - * @p{omega} is the relaxation - * parameter. - */ - template - void SSOR (Vector &v, - const number omega = 1.) const; - - /** - * Perform an SOR preconditioning in-place. - * The result is $v = (\omega D - L)^{-1} v$. - * @p{omega} is the damping parameter. - */ - template - void SOR (Vector &v, - const number om = 1.) const; - - /** - * Perform a transpose SOR preconditioning in-place. - * The result is $v = (\omega D - L)^{-1} v$. - * @p{omega} is the damping parameter. - */ - template - void TSOR (Vector &v, - const number om = 1.) const; - - /** - * Do one SOR step on @p{v}. - * Performs a direct SOR step - * with right hand side @p{b}. - */ - template - void SOR_step (Vector &v, - const Vector &b, - const number om = 1.) const; - - /** - * Do one adjoint SOR step on - * @p{v}. Performs a direct TSOR - * step with right hand side @p{b}. - */ - template - void TSOR_step (Vector &v, - const Vector &b, - const number om = 1.) const; - - /** - * Do one adjoint SSOR step on - * @p{v}. Performs a direct SSOR - * step with right hand side @p{b} - * by performing TSOR after SOR. - */ - template - void SSOR_step (Vector &v, - const Vector &b, - const number om = 1.) const; - /** * Print the matrix to the given * stream, using the format @@ -662,29 +654,19 @@ class SparseMatrixEZ : public Subscriptor * of this object. */ unsigned int memory_consumption () const; - - /** - * Exception - */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The entry with index <" << arg1 << ',' << arg2 - << "> does not exist."); /** - * Exception + * Exception for applying + * inverse-type operators to + * rectangular matrices. */ - DeclException0 (ExcMatrixNotSquare); - + DeclException0(ExcNoSquareMatrix); + /** - * Exception + * Exception for missing diagonal entry. */ - DeclException2 (ExcIteratorRange, - int, int, - << "The iterators denote a range of " << arg1 - << " elements, but the given number of rows was " << arg2); + DeclException0(ExcNoDiagonal); - private: /** @@ -694,6 +676,7 @@ class SparseMatrixEZ : public Subscriptor */ const Entry* locate (unsigned int row, unsigned int col) const; + /** * Find an entry or generate it. */ @@ -776,10 +759,9 @@ class SparseMatrixEZ : public Subscriptor unsigned int n_columns; /** - * Start of indices rows. Points - * into the data field. + * Info structure for each row. */ - std::vector row_start; + std::vector row_info; /** * Data storage. @@ -814,28 +796,34 @@ template inline SparseMatrixEZ::Entry::Entry() : - column(invalid_entry), + column(invalid), value(0) {} -template -inline -bool -SparseMatrixEZ::Entry::operator==(const Entry& e) const -{ - return column == e.column; -} +// template +// inline +// bool +// SparseMatrixEZ::Entry::operator==(const Entry& e) const +// { +// return column == e.column; +// } +// template +// inline +// bool +// SparseMatrixEZ::Entry::operator<(const Entry& e) const +// { +// return column < e.column; +// } + template inline -bool -SparseMatrixEZ::Entry::operator<(const Entry& e) const -{ - return column < e.column; -} - +SparseMatrixEZ::RowInfo::RowInfo(unsigned int start) + : + start(start), length(0), diagonal(invalid_diagonal) +{} //----------------------------------------------------------------------// @@ -843,7 +831,7 @@ template inline unsigned int SparseMatrixEZ::m () const { - return row_start.size() - 1; + return row_info.size(); }; @@ -864,13 +852,14 @@ const SparseMatrixEZ::Entry* SparseMatrixEZ::locate ( Assert (rowcolumn == col) return entry; - if (entry->column == Entry::invalid_entry) + if (entry->column == Entry::invalid) return 0; } return 0; @@ -887,63 +876,82 @@ SparseMatrixEZ::Entry* SparseMatrixEZ::allocate ( Assert (row= row) + i += r.diagonal; + // Find position of entry + while (icolumn == col) + return entry; + + // Now, we must insert the new + // entry and move all successive + // entries back. + + // If no more space is available + // for this row, insert new + // elements into the vector. + if (row != row_info.size()-1) { - Entry* entry = &data[i]; - // entry found - if (entry->column == col) - return entry; - // entry does not exist, - // create it - if (entry->column > col) + if (end >= row_info[row+1].start) + { + // Insert new entries + data.insert(data.begin()+end, increment, Entry()); + entry = &data[i]; + // Update starts of + // following rows + for (unsigned int rn=row+1;rn= data.size()) { - // 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;jcolumn = col; + entry->value = 0; + // Update row_info + ++r.length; + if (col == row) + r.diagonal = i - r.start; + else if (col::reinit(unsigned int n_rows, unsigned int default_row_length, unsigned int default_increment) { - if (default_row_length == Entry::invalid_entry) + clear(); + if (default_row_length == Entry::invalid) default_row_length = 5; - if (default_increment == Entry::invalid_entry) + if (default_increment == Entry::invalid) default_increment = 4; if (default_increment == 0) default_increment = 4; increment = default_increment; n_columns = n_cols; - row_start.resize(n_rows+1); + row_info.resize(n_rows); 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; + for (unsigned int i=0;i::clear() { n_columns = 0; - row_start.resize(0); + row_info.resize(0); data.resize(0); } @@ -94,18 +95,20 @@ SparseMatrixEZ::vmult (Vector& dst, Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); - const unsigned int end_row = row_start.size() - 1; - for (unsigned int i = 0; i < end_row;++i) + const unsigned int end_row = row_info.size(); + for (unsigned int row = 0; row < end_row; ++row) { - unsigned int index = row_start[i]; - unsigned int end = row_start[i+1]; + const RowInfo& ri = row_info[row]; + typename std::vector::const_iterator + entry = data.begin() + ri.start; double s = 0.; - for (;index != end && data[index].column != Entry::invalid_entry;++index) + for (unsigned short i=0;icolumn != Entry::invalid, + ExcInternalError()); + s += entry->value * src(entry->column); } - dst(i) = s; + dst(row) = s; } } @@ -130,18 +133,20 @@ SparseMatrixEZ::vmult_add (Vector& dst, Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); - const unsigned int end_row = row_start.size() - 1; - for (unsigned int i = 0; i < end_row;++i) + const unsigned int end_row = row_info.size(); + for (unsigned int row = 0; row < end_row; ++row) { - unsigned int index = row_start[i]; - unsigned int end = row_start[i+1]; + const RowInfo& ri = row_info[i]; + typename std::vector::const_iterator + entry = data.begin() + ri.start; double s = 0.; - for (;index != end && data[index].column != Entry::invalid_entry;++index) + for (unsigned short i=0;icolumn != Entry::invalid, + ExcInternalError()); + s += entry->value * src(entry->column); } - dst(i) += s; + dst(row) += s; } } @@ -154,19 +159,157 @@ SparseMatrixEZ::Tvmult_add (Vector& dst, Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size())); Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size())); - const unsigned int end_row = row_start.size() - 1; - for (unsigned int i = 0; i < end_row;++i) + const unsigned int end_row = row_info.size(); + for (unsigned int row = 0; row < end_row; ++row) { - unsigned int index = row_start[i]; - unsigned int end = row_start[i+1]; - for (;index != end && data[index].column != Entry::invalid_entry;++index) + const RowInfo& ri = row_info[row]; + typename std::vector::const_iterator + entry = data.begin() + ri.start; + for (unsigned short i=0;icolumn != Entry::invalid, + ExcInternalError()); + dst(entry->column) += entry->value * src(row); } } } + +template +template +void +SparseMatrixEZ::precondition_Jacobi (Vector &dst, + const Vector &src, + const number om) const +{ + Assert (m() == n(), ExcNoSquareMatrix()); + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + somenumber *dst_ptr = dst.begin(); + const somenumber *src_ptr = src.begin(); + typename std::vector::const_iterator ri = row_info.begin(); + const typename std::vector::const_iterator end = row_info.end(); + + for (; ri != end; ++dst_ptr, ++src_ptr, ++ri) + { + Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal()); + *dst_ptr = om * *src_ptr / data[ri->start + ri->diagonal].value; + } +}; + + + +template +template +void +SparseMatrixEZ::precondition_SOR (Vector &dst, + const Vector &src, + const number om) const +{ + Assert (m() == n(), ExcNoSquareMatrix()); + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + somenumber *dst_ptr = dst.begin(); + const somenumber *src_ptr = src.begin(); + typename std::vector::const_iterator ri = row_info.begin(); + const typename std::vector::const_iterator end = row_info.end(); + + for (; ri != end; ++dst_ptr, ++src_ptr, ++ri) + { + Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal()); + number s = *src_ptr; + const unsigned int end_row = ri->start + ri->diagonal; + for (unsigned int i=ri->start;istart + ri->diagonal].value; + } +}; + + + +template +template +void +SparseMatrixEZ::precondition_TSOR (Vector &dst, + const Vector &src, + const number om) const +{ + Assert (m() == n(), ExcNoSquareMatrix()); + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + somenumber *dst_ptr = dst.begin()+dst.size()-1; + const somenumber *src_ptr = src.begin()+src.size()-1; + typename std::vector::const_reverse_iterator + ri = row_info.rbegin(); + const typename std::vector::const_reverse_iterator + end = row_info.rend(); + + for (; ri != end; --dst_ptr, --src_ptr, ++ri) + { + Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal()); + number s = *src_ptr; + const unsigned int end_row = ri->start + ri->length; + for (unsigned int i=ri->start+ri->diagonal+1;istart + ri->diagonal].value; + } +}; + + + +template +template +void +SparseMatrixEZ::precondition_SSOR (Vector &dst, + const Vector &src, + const number om) const +{ + Assert (m() == n(), ExcNoSquareMatrix()); + Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); + Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); + + somenumber *dst_ptr = dst.begin(); + const somenumber *src_ptr = src.begin(); + typename std::vector::const_iterator ri; + const typename std::vector::const_iterator end = row_info.end(); + + // Forward + for (ri = row_info.begin(); ri != end; ++dst_ptr, ++src_ptr, ++ri) + { + Assert (ri->diagonal != RowInfo::invalid_diagonal, ExcNoDiagonal()); + number s = *src_ptr; + const unsigned int end_row = ri->start + ri->diagonal; + for (unsigned int i=ri->start;istart + ri->diagonal].value; + } + // Diagonal + dst_ptr = dst.begin(); + for (ri = row_info.begin(); ri != end; ++dst_ptr, ++ri) + *dst_ptr *= (2.-om) * data[ri->start + ri->diagonal].value; + + // Backward + typename std::vector::const_reverse_iterator rri; + const typename std::vector::const_reverse_iterator + rend = row_info.rend(); + dst_ptr = dst.begin()+dst.size()-1; + for (rri = row_info.rbegin(); rri != rend; --dst_ptr, ++rri) + { + const unsigned int end_row = rri->start + rri->length; + for (unsigned int i=rri->start+rri->diagonal+1;istart + rri->diagonal].value; + } +}; + + + template unsigned int SparseMatrixEZ::memory_consumption() const diff --git a/tests/lac/sparse_matrices.cc b/tests/lac/sparse_matrices.cc index df6e0b4b08..0775afe138 100644 --- a/tests/lac/sparse_matrices.cc +++ b/tests/lac/sparse_matrices.cc @@ -37,28 +37,31 @@ check_vmult_quadratic(const MATRIX& A, Vector f(A.m()); GrowingVectorMemory<> mem; - SolverControl control(20, 1.e-3, false); + SolverControl control(10, 1.e-13, false); SolverRichardson<> rich(control, mem, .01); + SolverRichardson<> prich(control, mem, 1.); PreconditionIdentity prec; - + PreconditionJacobi jacobi; + jacobi.initialize(A, .5); + PreconditionSOR sor; + sor.initialize(A, 1.2); + PreconditionSSOR ssor; + ssor.initialize(A, 1.2); + u = 0.; f = 1.; - try - { - rich.solve(A, u, f, prec); - } - catch (...) - { - } + 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 (...) {} + + u = 0.; deallog << "Transpose" << std::endl; - try - { - rich.Tsolve(A, u, f, prec); - } - catch (...) - { - } + 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 (...) {} deallog.pop(); } @@ -72,13 +75,15 @@ int main() // Switch between regression test // and benchmark #ifdef DEBUG - deallog.depth_console(0); - const unsigned int size = 10; + deallog.depth_console(3); + const unsigned int size = 5; + const unsigned int row_length = 3; #else deallog.depth_console(1000); deallog.log_execution_time(true); deallog.log_time_differences(true); const unsigned int size = 500; + const unsigned int row_length = 9; #endif FDMatrix testproblem (size, size); @@ -90,12 +95,12 @@ int main() structure.compress(); SparseMatrix A(structure); deallog << "Assemble" << std::endl; - testproblem.five_point(A); + testproblem.five_point(A, true); check_vmult_quadratic(A, "5-SparseMatrix"); - SparseMatrixEZ E(dim,dim,5,1); + SparseMatrixEZ E(dim,dim,row_length,2); deallog << "Assemble" << std::endl; - testproblem.five_point(E); + testproblem.five_point(E, true); check_vmult_quadratic(E, "5-SparseMatrixEZ"); A.clear(); @@ -109,7 +114,7 @@ int main() check_vmult_quadratic(A, "9-SparseMatrix"); E.clear(); - E.reinit(dim,dim,9,2); + E.reinit(dim,dim,row_length,2); deallog << "Assemble" << std::endl; testproblem.nine_point(E); check_vmult_quadratic(E, "9-SparseMatrixEZ"); diff --git a/tests/lac/testmatrix.cc b/tests/lac/testmatrix.cc index 5b13362163..7ca6bdb4c3 100644 --- a/tests/lac/testmatrix.cc +++ b/tests/lac/testmatrix.cc @@ -114,7 +114,7 @@ FDMatrix::nine_point_structure(SparsityPattern& structure) const template void -FDMatrix::nine_point(MATRIX& A) const +FDMatrix::nine_point(MATRIX& A, bool nonsymmetric) const { for(unsigned int i=0;i<=ny-2;i++) { @@ -170,7 +170,7 @@ FDMatrix::nine_point(MATRIX& A) const template void -FDMatrix::five_point(MATRIX& A) const +FDMatrix::five_point(MATRIX& A, bool nonsymmetric) const { for(unsigned int i=0;i<=ny-2;i++) { @@ -178,10 +178,16 @@ 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 (nonsymmetric) + A.set(row, row, 5.); + else + A.set(row, row, 4.); if (j>0) { - A.set(row-1, row, -1.); + if (nonsymmetric) + A.set(row-1, row, -2.); + else + A.set(row-1, row, -1.); A.set(row, row-1, -1.); } if (j - void five_point(MATRIX&) const; + void five_point(MATRIX&, bool nonsymmetric = false) const; /** * Fill the matrix with values. */ template - void nine_point(MATRIX&) const; + void nine_point(MATRIX&, bool nonsymmetric = false) const; template void gnuplot_print(std::ostream&, const Vector&) const; @@ -52,4 +52,3 @@ class FDMatrix */ unsigned int ny; }; -