From 91596ab368feb14a252f60e9c66891fc467a0ae0 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 25 Jun 1999 15:56:46 +0000 Subject: [PATCH] Direct implementation of SOR, TSOR, SSOR git-svn-id: https://svn.dealii.org/trunk@1470 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/sparse_matrix.2.templates | 44 +++++ deal.II/lac/include/lac/sparse_matrix.h | 72 +++++++- .../lac/include/lac/sparse_matrix.templates.h | 108 +++++++++++ deal.II/lac/source/sparse_matrix.double.cc | 67 +------ deal.II/lac/source/sparse_matrix.float.cc | 67 +------ .../lac/source/sparse_matrix.long_double.cc | 171 +----------------- 6 files changed, 228 insertions(+), 301 deletions(-) create mode 100644 deal.II/lac/include/lac/sparse_matrix.2.templates diff --git a/deal.II/lac/include/lac/sparse_matrix.2.templates b/deal.II/lac/include/lac/sparse_matrix.2.templates new file mode 100644 index 0000000000..9dbccff468 --- /dev/null +++ b/deal.II/lac/include/lac/sparse_matrix.2.templates @@ -0,0 +1,44 @@ +// $Id$ +// Driver file for SparseMatrix functions with two types. + +// TYPEMAT and TYPE2 are defined in sparsematrix?.cc + +template SparseMatrix & +SparseMatrix::copy_from (const SparseMatrix &); + +template void SparseMatrix::add_scaled (const TYPEMAT, + const SparseMatrix &); + +template void SparseMatrix::vmult (Vector &, + const Vector &) const; +template void SparseMatrix::Tvmult (Vector &, + const Vector &) const; +template void SparseMatrix::vmult_add (Vector &, + const Vector &) const; +template void SparseMatrix::Tvmult_add (Vector &, + const Vector &) const; + +template TYPE2 +SparseMatrix::matrix_norm (const Vector &) const; + +template TYPE2 SparseMatrix::residual (Vector &, + const Vector &, + const Vector &) const; + +template void SparseMatrix::precondition_SSOR (Vector &, + const Vector &, + TYPEMAT) const; + +template void SparseMatrix::precondition_SOR (Vector &, + const Vector &, + TYPEMAT) const; + +template void SparseMatrix::precondition_Jacobi (Vector &, + const Vector &, + TYPEMAT) const; + +template void SparseMatrix::SOR (Vector &, TYPEMAT) const; +template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +template void SparseMatrix::SOR_step (Vector &, const Vector &, TYPEMAT) const; +template void SparseMatrix::TSOR_step (Vector &, const Vector &, TYPEMAT) const; +template void SparseMatrix::SSOR_step (Vector &, const Vector &, TYPEMAT) const; diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 94c277dc91..f16d3070cd 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -834,6 +834,22 @@ class SparseMatrix : public Subscriptor template void Tvmult (Vector& dst, const Vector& src) const; + /** + * Adding Matrix-vector multiplication. Add + * $M*src$ on $dst$ with $M$ being this matrix. + */ + template + void vmult_add (Vector& dst, const Vector& src) const; + + /** + * Adding Matrix-vector multiplication. Add + * $M^T*src$ to $dst$ with $M$ being this + * matrix. This function does the same as + * #vmult_add# but takes the transposed matrix. + */ + template + void Tvmult_add (Vector& dst, const Vector& src) const; + /** * Return the norm of the vector $v$ with * respect to the norm induced by this @@ -902,35 +918,73 @@ class SparseMatrix : public Subscriptor void precondition_Jacobi (Vector &dst, const Vector &src, const number om = 1.) const; - // + /** + * Apply SOR preconditioning to #src#. + */ template void precondition_SSOR (Vector &dst, const Vector &src, const number om = 1.) const; - // + + /** + * Apply SOR preconditioning matrix to #src#. + * The result of this method is + * $dst = (om D - L)^{-1} src$. + */ template void precondition_SOR (Vector &dst, const Vector &src, const number om = 1.) const; /** - * Perform an SSOR step in-place, i.e. - * without copying to a second vector. - * #omega# is the damping parameter. + * Perform SSOR preconditioning in-place. + * Apply the preconditioner matrix without + * copying to a second vector. + * #omega# is the relaxation parameter. */ template - void SSOR (Vector &dst, + void SSOR (Vector &v, const number omega = 1.) const; /** - * Perform an SOR step in-place, i.e. - * without copying to a second vector. + * Perform an SOR preconditioning in-place. + * The result is $v = (\omega D - L)^{-1} v$. * #omega# is the damping parameter. */ template - void SOR (Vector &dst, + void SOR (Vector &v, const number om = 1.) const; + /** + * Do one SOR step on #v#. + * Performs a direct SOR step with right hand + * side #b#. + */ + template + void SOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one adjoint SOR step on #v#. + * Performs a direct TSOR step with right hand + * side #b#. + */ + template + void TSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + + /** + * Do one adjoint SSOR step on #v#. + * Performs a direct SSOR step with right hand + * side #b# by performing TSOR after SOR. + */ + template + void SSOR_step (Vector &v, + const Vector &b, + const number om = 1.) const; + /** * Return a (constant) reference to the * underlying sparsity pattern of this diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 5818de14c4..3ebf0eca2f 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -226,6 +226,53 @@ SparseMatrix::Tvmult (Vector& dst, const Vector& +template +template +void +SparseMatrix::vmult_add (Vector& dst, const Vector& src) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size())); + Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size())); + + const unsigned int n_rows = m(); + const number *val_ptr = &val[cols->rowstart[0]]; + const int *colnum_ptr = &cols->colnums[cols->rowstart[0]]; + somenumber *dst_ptr = &dst(0); + for (unsigned int row=0; rowrowstart[row+1]]; + while (val_ptr != val_end_of_row) + s += *val_ptr++ * src(*colnum_ptr++); + *dst_ptr++ += s; + }; +}; + + +template +template +void +SparseMatrix::Tvmult_add (Vector& dst, const Vector& src) const +{ + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert(n() == dst.size(), ExcDimensionsDontMatch(n(),dst.size())); + Assert(m() == src.size(), ExcDimensionsDontMatch(m(),src.size())); + + for (unsigned int i=0;irowstart[i]; jrowstart[i+1] ;j++) + { + int p = cols->colnums[j]; + dst(p) += val[j] * src(i); + } + } +} + + + template template somenumber @@ -452,6 +499,65 @@ SparseMatrix::SOR (Vector& dst, const number om) const } +template +template +void +SparseMatrix::SOR_step (Vector &v, + const Vector &b, + const number om) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (m() == v.size(), ExcDimensionsDontMatch(m(),v.size())); + Assert (m() == b.size(), ExcDimensionsDontMatch(m(),b.size())); + + for (unsigned int row=0; rowrowstart[row]; jrowstart[row+1]; ++j) + { + s -= val[j] * v(cols->colnums[j]); + } + v(row) += s * om / val[cols->rowstart[row]]; + } +} + +template +template +void +SparseMatrix::TSOR_step (Vector &v, + const Vector &b, + const number om) const +{ + Assert (cols != 0, ExcMatrixNotInitialized()); + Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (m() == v.size(), ExcDimensionsDontMatch(m(),v.size())); + Assert (m() == b.size(), ExcDimensionsDontMatch(m(),b.size())); + + for (int row=m()-1; row>=0; --row) + { + somenumber s = b(row); + for (unsigned int j=cols->rowstart[row]; jrowstart[row+1]; ++j) + { + s -= val[j] * v(cols->colnums[j]); + } + v(row) += s * om / val[cols->rowstart[row]]; + } +} + +template +template +void +SparseMatrix::SSOR_step (Vector &v, + const Vector &b, + const number om) const +{ + SOR_step(v,b,om); + TSOR_step(v,b,om); +} + template template @@ -460,6 +566,8 @@ SparseMatrix::SSOR (Vector& dst, const number om) const { Assert (cols != 0, ExcMatrixNotInitialized()); Assert (val != 0, ExcMatrixNotInitialized()); + Assert (m() == n(), ExcMatrixNotSquare()); + Assert (m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size())); int p; const unsigned int n = dst.size(); diff --git a/deal.II/lac/source/sparse_matrix.double.cc b/deal.II/lac/source/sparse_matrix.double.cc index 3ad663c052..d129b74e99 100644 --- a/deal.II/lac/source/sparse_matrix.double.cc +++ b/deal.II/lac/source/sparse_matrix.double.cc @@ -18,71 +18,14 @@ template class SparseMatrix; #define TYPE2 float -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include #undef TYPE2 #define TYPE2 double -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); +#include -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; +#undef TYPE2 +#define TYPE2 long double -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include diff --git a/deal.II/lac/source/sparse_matrix.float.cc b/deal.II/lac/source/sparse_matrix.float.cc index 88c4460cf0..f7d2f03495 100644 --- a/deal.II/lac/source/sparse_matrix.float.cc +++ b/deal.II/lac/source/sparse_matrix.float.cc @@ -18,71 +18,14 @@ template class SparseMatrix; #define TYPE2 float -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include #undef TYPE2 #define TYPE2 double -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); +#include -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template double -SparseMatrix::matrix_norm (const Vector &) const; - -template double SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; +#undef TYPE2 +#define TYPE2 long double -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include diff --git a/deal.II/lac/source/sparse_matrix.long_double.cc b/deal.II/lac/source/sparse_matrix.long_double.cc index 7334d98cce..dbd0e62aa0 100644 --- a/deal.II/lac/source/sparse_matrix.long_double.cc +++ b/deal.II/lac/source/sparse_matrix.long_double.cc @@ -18,179 +18,14 @@ template class SparseMatrix; #define TYPE2 float -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include #undef TYPE2 #define TYPE2 double -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include #undef TYPE2 #define TYPE2 long double -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; - -#undef TYPEMAT -#define TYPEMAT double - -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; - -#undef TYPEMAT -#define TYPEMAT float - -template SparseMatrix & -SparseMatrix::copy_from (const SparseMatrix &); - -template void SparseMatrix::add_scaled (const TYPEMAT, - const SparseMatrix &); - -template void SparseMatrix::vmult (Vector &, - const Vector &) const; -template void SparseMatrix::Tvmult (Vector &, - const Vector &) const; - -template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; - -template TYPE2 SparseMatrix::residual (Vector &, - const Vector &, - const Vector &) const; - -template void SparseMatrix::precondition_SSOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_SOR (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::precondition_Jacobi (Vector &, - const Vector &, - TYPEMAT) const; - -template void SparseMatrix::SOR (Vector &, TYPEMAT) const; -template void SparseMatrix::SSOR (Vector &, TYPEMAT) const; +#include -- 2.39.5