From 23336544be13c489286df9bca9e58814936fc0ed Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Sun, 18 Nov 2007 22:20:58 +0000 Subject: [PATCH] First batch of changes to make the FullMatrix deal with complex numbers git-svn-id: https://svn.dealii.org/trunk@15517 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/full_matrix.h | 55 ++++++-- .../lac/include/lac/full_matrix.templates.h | 118 ++++++++++-------- deal.II/lac/source/full_matrix.double.cc | 107 ++++++++++++++++ 3 files changed, 217 insertions(+), 63 deletions(-) diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index e147ac1cb9..ca5e181a63 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -62,6 +63,29 @@ class FullMatrix : public Table<2,number> * the STL container classes. */ 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; + class const_iterator; @@ -430,8 +454,17 @@ class FullMatrix : public Table<2,number> * 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 number2 matrix_norm_square (const Vector &v) const; @@ -444,6 +477,12 @@ class FullMatrix : public Table<2,number> * the cellwise scalar product of * two functions in the finite * element context. + * + * 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 number2 matrix_scalar_product (const Vector &u, @@ -457,7 +496,7 @@ class FullMatrix : public Table<2,number> * |M_{ij}|$ (maximum of * the sums over columns). */ - number l1_norm () const; + real_type l1_norm () const; /** * Return the $l_\infty$-norm of @@ -466,7 +505,7 @@ class FullMatrix : public Table<2,number> * |M_{ij}|$ (maximum of the sums * over rows). */ - number linfty_norm () const; + real_type linfty_norm () const; /** * Compute the Frobenius norm of @@ -480,7 +519,7 @@ class FullMatrix : public Table<2,number> * l2-norm of * the vector space. */ - number frobenius_norm () const; + real_type frobenius_norm () const; /** * Compute the relative norm of @@ -495,7 +534,7 @@ class FullMatrix : public Table<2,number> * symmetric within a certain * accuracy, or not. */ - number relative_symmetry_norm2 () const; + real_type relative_symmetry_norm2 () const; /** * Computes the determinant of a @@ -531,7 +570,7 @@ class FullMatrix : public Table<2,number> /** * Print the matrix and allow * formatting of entries. - * + * * The parameters allow for a * flexible setting of the output * format: @@ -1141,7 +1180,7 @@ template FullMatrix & FullMatrix::operator = (const number d) { - Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); + Assert (d==number(0), ExcScalarAssignmentOnlyForZeroValue()); if (this->n_elements() != 0) std::fill_n (this->val, this->n_elements(), number()); diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 42df4a50bb..0d3407c153 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -314,9 +314,9 @@ FullMatrix::vmult (Vector& dst, { for (unsigned int i=0; i::vmult (Vector& dst, { for (unsigned int i=0; i::Tvmult (Vector &dst, { for (unsigned int i=0; iel(j,i); + s += number(src(j)) * this->el(j,i); dst(i) = s; }; } @@ -365,12 +365,12 @@ void FullMatrix::Tvmult (Vector &dst, { for (unsigned int i=0; iel(j,i); + s += number(src(j)) * this->el(j,i); dst(i) += s; - }; - }; + } + } } @@ -388,14 +388,14 @@ number FullMatrix::residual (Vector& dst, Assert (&src != &dst, ExcSourceEqualsDestination()); - number2 s, res = 0.; + number res = 0.; const unsigned int size_m = m(), size_n = n(); for (unsigned int i=0; iel(i,j); + s -= number(src(j)) * this->el(i,j); dst(i) = s; res += s*s; } @@ -416,11 +416,11 @@ void FullMatrix::forward (Vector &dst, unsigned int i,j; unsigned int nu = ( (m()el(i,j); + number s = src(i); + for (j=0; jel(i,j); dst(i) = s/this->el(i,i); } } @@ -436,12 +436,12 @@ void FullMatrix::backward (Vector &dst, unsigned int j; unsigned int nu = (m()=0; --i) { - s = src(i); - for (j=i+1; jel(i,j); - dst(i) = s/this->el(i,i); + number2 s = src(i); + for (j=i+1; jel(i,j)); + dst(i) = s/number2(this->el(i,i)); } } @@ -946,7 +946,8 @@ void FullMatrix::Tmmult (FullMatrix &dst, template template -number2 FullMatrix::matrix_norm_square (const Vector &v) const +number2 +FullMatrix::matrix_norm_square (const Vector &v) const { Assert (!this->empty(), ExcEmptyMatrix()); @@ -960,14 +961,14 @@ number2 FullMatrix::matrix_norm_square (const Vector &v) const for (unsigned int row=0; row::matrix_norm_square (const Vector &v) const template template -number2 FullMatrix::matrix_scalar_product (const Vector &u, +number2 +FullMatrix::matrix_scalar_product (const Vector &u, const Vector &v) const { Assert (!this->empty(), ExcEmptyMatrix()); @@ -991,14 +993,14 @@ number2 FullMatrix::matrix_scalar_product (const Vector &u, for (unsigned int row=0; row::symmetrize () for (unsigned int i=0; iel(i,j) + this->el(j,i)) / 2; + const number t = (this->el(i,j) + this->el(j,i)) / 2.; this->el(i,j) = this->el(j,i) = t; }; } template -number FullMatrix::l1_norm () const +typename FullMatrix::real_type +FullMatrix::l1_norm () const { Assert (!this->empty(), ExcEmptyMatrix()); - number sum=0, max=0; + real_type sum=0, max=0; const unsigned int n_rows = m(), n_cols = n(); for (unsigned int col=0; colel(row,col)); + sum += std::abs(this->el(row,col)); if (sum > max) max = sum; } @@ -1043,18 +1046,19 @@ number FullMatrix::l1_norm () const template -number FullMatrix::linfty_norm () const +typename FullMatrix::real_type +FullMatrix::linfty_norm () const { Assert (!this->empty(), ExcEmptyMatrix()); - number sum=0, max=0; + real_type sum=0, max=0; const unsigned int n_rows = m(), n_cols = n(); for (unsigned int row=0; rowel(row,col)); + sum += std::abs(this->el(row,col)); if (sum > max) max = sum; } @@ -1950,33 +1954,37 @@ FullMatrix::trace () const template -number +typename FullMatrix::real_type FullMatrix::frobenius_norm () const { Assert (!this->empty(), ExcEmptyMatrix()); - number s = 0.; + real_type s = 0.; for (unsigned int i=0; in_rows()*this->n_cols(); ++i) - s += this->data()[i]*this->data()[i]; + s += numbers::NumberTraits::abs_square(this->data()[i]); return std::sqrt(s); } template -number +typename FullMatrix::real_type FullMatrix::relative_symmetry_norm2 () const { Assert (!this->empty(), ExcEmptyMatrix()); - number s = 0.; - number a = 0.; + real_type s = 0.; + real_type a = 0.; for (unsigned int i=0; in_rows(); ++i) for (unsigned int j=0; jn_cols(); ++j) { - a += ((*this)(i,j)-(*this)(j,i))*((*this)(i,j)-(*this)(j,i)); - s += (*this)(i,j)*(*this)(i,j); + const number x_ij = (*this)(i,j); + const number x_ji = (*this)(j,i); + + a += numbers::NumberTraits::abs_square(x_ij-x_ji); + s += numbers::NumberTraits::abs_square(x_ij); } + if (s!=0.) return std::sqrt(a)/std::sqrt(s); return 0; @@ -2084,7 +2092,7 @@ FullMatrix::invert (const FullMatrix &M) const number t63 = t43*t20-t43*t22-t46*t33+t46*t35+ t49*t50-t49*t52-t54*t25+t54*t27+ t57*t38-t57*t40-t60*t50+t60*t52; - const number t65 = 1/(t42+t63); + const number t65 = 1./(t42+t63); const number t71 = M.el(0,2)*M.el(2,1); const number t73 = M.el(0,3)*M.el(2,1); const number t75 = M.el(0,2)*M.el(3,1); @@ -2173,7 +2181,7 @@ FullMatrix::precondition_Jacobi (Vector &dst, const somenumber *src_ptr = src.begin(); for (unsigned int i=0; iel(i,i); + *dst_ptr = somenumber(om) * *src_ptr / somenumber(this->el(i,i)); } @@ -2213,7 +2221,7 @@ FullMatrix::print_formatted ( for (unsigned int i=0; iel(i,j)) > threshold) + if (std::abs(this->el(i,j)) > threshold) out << std::setw(width) << this->el(i,j) * denominator << ' '; else @@ -2249,7 +2257,7 @@ FullMatrix::gauss_jordan () // regular double diagonal_sum = 0; for (unsigned int i=0; iel(i,i)); + diagonal_sum += std::abs(this->el(i,i)); const double typical_diagonal_element = diagonal_sum/N; // initialize the array that holds @@ -2265,13 +2273,13 @@ FullMatrix::gauss_jordan () // part of the line on and // right of the diagonal for // the largest element - number max = std::fabs(this->el(j,j)); + real_type max = std::abs(this->el(j,j)); unsigned int r = j; for (unsigned int i=j+1; iel(i,j)) > max) + if (std::abs(this->el(i,j)) > max) { - max = std::fabs(this->el(i,j)); + max = std::abs(this->el(i,j)); r = i; } } diff --git a/deal.II/lac/source/full_matrix.double.cc b/deal.II/lac/source/full_matrix.double.cc index 29c5e79095..0361ce0a50 100644 --- a/deal.II/lac/source/full_matrix.double.cc +++ b/deal.II/lac/source/full_matrix.double.cc @@ -117,4 +117,111 @@ FullMatrix::residual(Vector&, const Vector&) const; + +// #undef TYPEMAT + +// #define TYPEMAT std::complex + +// template class FullMatrix; + +// template void FullMatrix::print( +// LogStream&, const unsigned int, const unsigned int) const; +// template void FullMatrix::print( +// std::ostream&, const unsigned int, const unsigned int) const; + +// template FullMatrix& FullMatrix::operator =( +// const FullMatrix >&); + +// #undef TYPEMAT2 +// #define TYPEMAT2 std::complex + +// template void FullMatrix::fill ( +// const FullMatrix&, unsigned, unsigned, unsigned, unsigned); +// template void FullMatrix::add (const TYPEMAT, const FullMatrix&); +// template void FullMatrix::add (const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&); +// template void FullMatrix::add (const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&); +// template void FullMatrix::add ( +// const FullMatrix&, std::complex, unsigned, unsigned, unsigned, unsigned); +// template void FullMatrix::Tadd (const TYPEMAT, const FullMatrix&); +// template void FullMatrix::Tadd ( +// const FullMatrix&, std::complex, unsigned, unsigned, unsigned, unsigned); +// template void FullMatrix::equ (const TYPEMAT, const FullMatrix&); +// template void FullMatrix::equ (const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&); +// template void FullMatrix::equ (const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&, +// const TYPEMAT, const FullMatrix&); +// template void FullMatrix::mmult (FullMatrix&, const FullMatrix&, const bool) const; +// template void FullMatrix::Tmmult (FullMatrix&, const FullMatrix&, const bool) const; +// template void FullMatrix::add_diag (const TYPEMAT, const FullMatrix&); +// template void FullMatrix::invert (const FullMatrix&); + + +// #undef TYPEVEC +// #undef TYPERES +// #define TYPEVEC std::complex +// #define TYPERES std::complex + +// template void FullMatrix::fill_permutation ( +// const FullMatrix&, +// const std::vector&, +// const std::vector&); +// template void FullMatrix::vmult( +// Vector&, const Vector&, bool) const; +// template void FullMatrix::Tvmult( +// Vector&, const Vector&, bool) const; +// template std::complex FullMatrix::residual( +// Vector&, const Vector&, const Vector&) const; +// template TYPEVEC FullMatrix::matrix_norm_square ( +// const Vector &) const; +// template TYPEVEC FullMatrix::matrix_scalar_product( +// const Vector&, const Vector&) const; +// template void FullMatrix::forward( +// Vector&, const Vector&) const; +// template void FullMatrix::backward( +// Vector&, const Vector&) const; + +// template +// void FullMatrix::precondition_Jacobi ( +// Vector &, const Vector &, const TYPEMAT) const; + +// #undef TYPEVEC +// #define TYPEVEC std::complex + +// template void FullMatrix::fill_permutation ( +// const FullMatrix&, +// const std::vector&, +// const std::vector&); +// template void FullMatrix::vmult( +// Vector&, const Vector&, bool) const; +// template void FullMatrix::Tvmult( +// Vector&, const Vector&, bool) const; +// template std::complex FullMatrix::residual( +// Vector&, const Vector&, const Vector&) const; +// template TYPEVEC FullMatrix::matrix_norm_square ( +// const Vector &) const; +// template TYPEVEC FullMatrix::matrix_scalar_product( +// const Vector&, const Vector&) const; +// template void FullMatrix::forward( +// Vector&, const Vector&) const; +// template void FullMatrix::backward( +// Vector&, const Vector&) const; +// template +// void FullMatrix::precondition_Jacobi ( +// Vector &, const Vector &, const TYPEMAT) const; + + +// #undef TYPERES +// #define TYPERES std::complex + +// template +// std::complex +// FullMatrix::residual(Vector&, +// const Vector&, +// const Vector&) const; + + DEAL_II_NAMESPACE_CLOSE -- 2.39.5