From 847d7ec66f1e3b61f84805e0c9f5b50266f59af4 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 21 Apr 2009 13:44:29 +0000 Subject: [PATCH] After half a year, finally check in some more of the conversion that allows to use the SparseMatrix also for complex-valued template arguments. git-svn-id: https://svn.dealii.org/trunk@18675 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix.h | 43 ++++++- .../lac/include/lac/sparse_matrix.templates.h | 46 ++++---- deal.II/lac/source/sparse_matrix.inst.in | 109 ++++++++++++++++++ 3 files changed, 172 insertions(+), 26 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index f3948e2ee0..b6cf3d02ba 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -441,7 +441,7 @@ namespace internals * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). * - * @author Essentially everyone who has ever worked on deal.II, 1994-2004 + * @author Essentially everyone who has ever worked on deal.II, 1994-2007 */ template class SparseMatrix : public virtual Subscriptor @@ -453,6 +453,27 @@ class SparseMatrix : public virtual Subscriptor */ typedef number value_type; + /** + * Declare a type that has holds + * real-valued numbers with the + * same precision as the template + * argument to this class. If the + * template argument of this + * class is a real data type, + * then real_type equals the + * template argument. If the + * template argument is a + * std::complex type then + * real_type equals the type + * underlying the complex + * numbers. + * + * This typedef is used to + * represent the return type of + * norms. + */ + typedef typename numbers::NumberTraits::real_type real_type; + /** * Typedef of an STL conforming iterator * class walking over all the nonzero @@ -1376,8 +1397,17 @@ class SparseMatrix : public virtual Subscriptor * nodal values of the finite * element function. * - * Obviously, the matrix needs to - * be quadratic for this operation. + * Obviously, the matrix needs to be + * quadratic for this operation, and for + * the result to actually be a norm it + * also needs to be either real symmetric + * or complex hermitian. + * + * The underlying template types of both + * this matrix and the given vector + * should either both be real or + * complex-valued, but not mixed, for + * this function to make sense. */ template somenumber matrix_norm_square (const Vector &v) const; @@ -1389,6 +1419,7 @@ class SparseMatrix : public virtual Subscriptor template somenumber matrix_scalar_product (const Vector &u, const Vector &v) const; + /** * Compute the residual of an * equation Mx=b, where @@ -1426,7 +1457,7 @@ class SparseMatrix : public virtual Subscriptor * $|Mv|_1\leq |M|_1 |v|_1$. * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) */ - number l1_norm () const; + real_type l1_norm () const; /** * Return the linfty-norm of the @@ -1440,7 +1471,7 @@ class SparseMatrix : public virtual Subscriptor * $|Mv|_infty \leq |M|_infty |v|_infty$. * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) */ - number linfty_norm () const; + real_type linfty_norm () const; /** * Return the frobenius norm of the @@ -1448,7 +1479,7 @@ class SparseMatrix : public virtual Subscriptor * sum of squares of all entries in the * matrix. */ - number frobenius_norm () const; + real_type frobenius_norm () const; //@} /** * @name Preconditioning methods diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index dc3b73c97c..3fe59f41ff 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -41,6 +41,7 @@ #include #include + DEAL_II_NAMESPACE_OPEN @@ -778,11 +779,11 @@ SparseMatrix::matrix_norm_square (const Vector& v) const while (val_ptr != val_end_of_row) s += *val_ptr++ * v(*colnum_ptr++); - sum += s* v(row); - }; + sum += v(row) * numbers::NumberTraits::conjugate(s); + } return sum; - }; + } } @@ -810,8 +811,9 @@ threaded_matrix_norm_square (const Vector &v, while (val_ptr != val_end_of_row) s += *val_ptr++ * v(*colnum_ptr++); - sum += s* v(row); - }; + sum += v(row) * numbers::NumberTraits::conjugate(s); + } + *partial_sum = sum; } @@ -891,11 +893,11 @@ SparseMatrix::matrix_scalar_product (const Vector& u, while (val_ptr != val_end_of_row) s += *val_ptr++ * v(*colnum_ptr++); - sum += s* u(row); - }; + sum += u(row) * numbers::NumberTraits::conjugate(s); + } return sum; - }; + } } @@ -924,8 +926,9 @@ threaded_matrix_scalar_product (const Vector &u, while (val_ptr != val_end_of_row) s += *val_ptr++ * v(*colnum_ptr++); - sum += s* u(row); - }; + sum += u(row) * numbers::NumberTraits::conjugate(s); + } + *partial_sum = sum; } @@ -934,16 +937,17 @@ threaded_matrix_scalar_product (const Vector &u, template -number SparseMatrix::l1_norm () const +typename SparseMatrix::real_type +SparseMatrix::l1_norm () const { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Vector column_sums(n()); + Vector column_sums(n()); const unsigned int n_rows = m(); for (unsigned int row=0; rowrowstart[row]; jrowstart[row+1] ; ++j) - column_sums(cols->colnums[j]) += std::fabs(val[j]); + column_sums(cols->colnums[j]) += numbers::NumberTraits::abs(val[j]); return column_sums.linfty_norm(); } @@ -951,21 +955,22 @@ number SparseMatrix::l1_norm () const template -number SparseMatrix::linfty_norm () const +typename SparseMatrix::real_type +SparseMatrix::linfty_norm () const { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); const number *val_ptr = &val[cols->rowstart[0]]; - number sum, max=0; + real_type max=0; const unsigned int n_rows = m(); for (unsigned int row=0; rowrowstart[row+1]]; while (val_ptr != val_end_of_row) - sum += std::fabs(*val_ptr++); + sum += numbers::NumberTraits::abs(*val_ptr++); if (sum > max) max = sum; } @@ -975,16 +980,17 @@ number SparseMatrix::linfty_norm () const template -number SparseMatrix::frobenius_norm () const +typename SparseMatrix::real_type +SparseMatrix::frobenius_norm () const { // simply add up all entries in the // sparsity pattern, without taking any // reference to rows or columns - number norm_sqr = 0; + real_type norm_sqr = 0; const unsigned int n_rows = m(); for (const number *ptr = &val[0]; ptr != &val[cols->rowstart[n_rows]]; ++ptr) - norm_sqr += *ptr * *ptr; + norm_sqr += numbers::NumberTraits::abs_square(*ptr); return std::sqrt (norm_sqr); } diff --git a/deal.II/lac/source/sparse_matrix.inst.in b/deal.II/lac/source/sparse_matrix.inst.in index 0edd36c17d..a8dafdb7cb 100644 --- a/deal.II/lac/source/sparse_matrix.inst.in +++ b/deal.II/lac/source/sparse_matrix.inst.in @@ -12,6 +12,8 @@ //---------------------------- sparse_matrix.inst.in --------------------------- +// real instantiations + for (S : REAL_SCALARS) { template class SparseMatrix; @@ -128,3 +130,110 @@ for (S1, S2, S3 : REAL_SCALARS; template void SparseMatrix:: Tvmult_add (V1 &, const V2 &) const; } + + + +// complex instantiations + +// for (S : COMPLEX_SCALARS) +// { +// template class SparseMatrix; +// } + + + +// for (S1, S2 : COMPLEX_SCALARS) +// { +// template SparseMatrix & +// SparseMatrix::copy_from (const SparseMatrix &); + +// template +// void SparseMatrix::copy_from (const FullMatrix &); + +// template void SparseMatrix::add (const S1, +// const SparseMatrix &); +// } + + +// for (S1, S2 : COMPLEX_SCALARS) +// { +// template S2 +// SparseMatrix:: +// matrix_norm_square (const Vector &) const; + +// template S2 +// SparseMatrix:: +// matrix_scalar_product (const Vector &, +// const Vector &) const; + +// template S2 SparseMatrix:: +// residual (Vector &, +// const Vector &, +// const Vector &) const; + +// template void SparseMatrix:: +// precondition_SSOR (Vector &, +// const Vector &, +// const S1) const; + +// template void SparseMatrix:: +// precondition_SOR (Vector &, +// const Vector &, +// const S1) const; + +// template void SparseMatrix:: +// precondition_TSOR (Vector &, +// const Vector &, +// const S1) const; + +// template void SparseMatrix:: +// precondition_Jacobi (Vector &, +// const Vector &, +// const S1) const; + +// template void SparseMatrix:: +// SOR (Vector &, +// const S1) const; +// template void SparseMatrix:: +// TSOR (Vector &, +// const S1) const; +// template void SparseMatrix:: +// SSOR (Vector &, +// const S1) const; +// template void SparseMatrix:: +// PSOR (Vector &, +// const std::vector&, +// const std::vector&, +// const S1) const; +// template void SparseMatrix:: +// TPSOR (Vector &, +// const std::vector&, +// const std::vector&, +// const S1) const; +// template void SparseMatrix:: +// SOR_step (Vector &, +// const Vector &, +// const S1) const; +// template void SparseMatrix:: +// TSOR_step (Vector &, +// const Vector &, +// const S1) const; +// template void SparseMatrix:: +// SSOR_step (Vector &, +// const Vector &, +// const S1) const; +// } + + +// for (S1, S2, S3 : COMPLEX_SCALARS; +// V1, V2 : DEAL_II_VEC_TEMPLATES) +// { +// template void SparseMatrix:: +// vmult (V1 &, const V2 &) const; +// template void SparseMatrix:: +// Tvmult (V1 &, const V2 &) const; +// template void SparseMatrix:: +// vmult_add (V1 &, const V2 &) const; +// template void SparseMatrix:: +// Tvmult_add (V1 &, const V2 &) const; +// } -- 2.39.5