#include <base/config.h>
+#include <base/numbers.h>
#include <base/table.h>
#include <lac/exceptions.h>
#include <lac/identity_matrix.h>
* 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<number>::real_type real_type;
+
class const_iterator;
* 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<typename number2>
number2 matrix_norm_square (const Vector<number2> &v) const;
* 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<typename number2>
number2 matrix_scalar_product (const Vector<number2> &u,
* |M_{ij}|$ (maximum of
* the sums over columns).
*/
- number l1_norm () const;
+ real_type l1_norm () const;
/**
* Return the $l_\infty$-norm of
* |M_{ij}|$ (maximum of the sums
* over rows).
*/
- number linfty_norm () const;
+ real_type linfty_norm () const;
/**
* Compute the Frobenius norm of
* <i>l<sub>2</sub></i>-norm of
* the vector space.
*/
- number frobenius_norm () const;
+ real_type frobenius_norm () const;
/**
* Compute the relative norm of
* symmetric within a certain
* accuracy, or not.
*/
- number relative_symmetry_norm2 () const;
+ real_type relative_symmetry_norm2 () const;
/**
* Computes the determinant of a
/**
* Print the matrix and allow
* formatting of entries.
- *
+ *
* The parameters allow for a
* flexible setting of the output
* format:
FullMatrix<number> &
FullMatrix<number>::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());
{
for (unsigned int i=0; i<size_m; ++i)
{
- number2 s = 0.;
+ number s = 0.;
for (unsigned int j=0; j<size_n; ++j)
- s += src(j) * *(e++);
+ s += number(src(j)) * *(e++);
dst(i) = s;
};
}
{
for (unsigned int i=0; i<size_m; ++i)
{
- number2 s = 0.;
+ number s = 0.;
for (unsigned int j=0; j<size_n; ++j)
- s += src(j) * *(e++);
+ s += number(src(j)) * *(e++);
dst(i) += s;
- };
- };
- };
+ }
+ }
+ }
}
{
for (unsigned int i=0; i<size_n; ++i)
{
- number2 s = 0.;
+ number s = 0.;
for (unsigned int j=0; j<size_m; ++j)
- s += src(j) * this->el(j,i);
+ s += number(src(j)) * this->el(j,i);
dst(i) = s;
};
}
{
for (unsigned int i=0; i<size_n; ++i)
{
- number2 s = 0.;
+ number s = 0.;
for (unsigned int j=0; j<size_m; ++j)
- s += src(j) * this->el(j,i);
+ s += number(src(j)) * this->el(j,i);
dst(i) += s;
- };
- };
+ }
+ }
}
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; i<size_n; ++i)
{
- s = right(i);
+ number s = right(i);
for (unsigned int j=0; j<size_m; ++j)
- s -= src(j) * this->el(i,j);
+ s -= number(src(j)) * this->el(i,j);
dst(i) = s;
res += s*s;
}
unsigned int i,j;
unsigned int nu = ( (m()<n()) ? m() : n());
- number2 s;
for (i=0; i<nu; ++i)
{
- s = src(i);
- for (j=0; j<i; ++j) s -= dst(j) * this->el(i,j);
+ number s = src(i);
+ for (j=0; j<i; ++j)
+ s -= number(dst(j)) * this->el(i,j);
dst(i) = s/this->el(i,i);
}
}
unsigned int j;
unsigned int nu = (m()<n() ? m() : n());
- number2 s;
for (int i=nu-1; i>=0; --i)
{
- s = src(i);
- for (j=i+1; j<nu; ++j) s -= dst(j) * this->el(i,j);
- dst(i) = s/this->el(i,i);
+ number2 s = src(i);
+ for (j=i+1; j<nu; ++j)
+ s -= dst(j) * number2(this->el(i,j));
+ dst(i) = s/number2(this->el(i,i));
}
}
template <typename number>
template <typename number2>
-number2 FullMatrix<number>::matrix_norm_square (const Vector<number2> &v) const
+number2
+FullMatrix<number>::matrix_norm_square (const Vector<number2> &v) const
{
Assert (!this->empty(), ExcEmptyMatrix());
for (unsigned int row=0; row<n_rows; ++row)
{
- number2 s = 0.;
+ number s = 0.;
const number * const val_end_of_row = val_ptr+n_rows;
v_ptr = v.begin();
while (val_ptr != val_end_of_row)
- s += *val_ptr++ * *v_ptr++;
+ s += number(*val_ptr++) * number(*v_ptr++);
- sum += s* v(row);
- };
+ sum += s* number(v(row));
+ }
return sum;
}
template <typename number>
template <typename number2>
-number2 FullMatrix<number>::matrix_scalar_product (const Vector<number2> &u,
+number2
+FullMatrix<number>::matrix_scalar_product (const Vector<number2> &u,
const Vector<number2> &v) const
{
Assert (!this->empty(), ExcEmptyMatrix());
for (unsigned int row=0; row<n_rows; ++row)
{
- number2 s = 0.;
+ number s = 0.;
const number * const val_end_of_row = val_ptr+n_cols;
v_ptr = v.begin();
while (val_ptr != val_end_of_row)
- s += *val_ptr++ * *v_ptr++;
+ s += number(*val_ptr++) * number(*v_ptr++);
- sum += s* u(row);
- };
+ sum += s * number(u(row));
+ }
return sum;
}
for (unsigned int i=0; i<N; ++i)
for (unsigned int j=i+1; j<N; ++j)
{
- const number t = (this->el(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 <typename number>
-number FullMatrix<number>::l1_norm () const
+typename FullMatrix<number>::real_type
+FullMatrix<number>::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; col<n_cols; ++col)
{
sum=0;
for (unsigned int row=0; row<n_rows; ++row)
- sum += std::fabs(this->el(row,col));
+ sum += std::abs(this->el(row,col));
if (sum > max)
max = sum;
}
template <typename number>
-number FullMatrix<number>::linfty_norm () const
+typename FullMatrix<number>::real_type
+FullMatrix<number>::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; row<n_rows; ++row)
{
sum=0;
for (unsigned int col=0; col<n_cols; ++col)
- sum += std::fabs(this->el(row,col));
+ sum += std::abs(this->el(row,col));
if (sum > max)
max = sum;
}
template <typename number>
-number
+typename FullMatrix<number>::real_type
FullMatrix<number>::frobenius_norm () const
{
Assert (!this->empty(), ExcEmptyMatrix());
- number s = 0.;
+ real_type s = 0.;
for (unsigned int i=0; i<this->n_rows()*this->n_cols(); ++i)
- s += this->data()[i]*this->data()[i];
+ s += numbers::NumberTraits<number>::abs_square(this->data()[i]);
return std::sqrt(s);
}
template <typename number>
-number
+typename FullMatrix<number>::real_type
FullMatrix<number>::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; i<this->n_rows(); ++i)
for (unsigned int j=0; j<this->n_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<number>::abs_square(x_ij-x_ji);
+ s += numbers::NumberTraits<number>::abs_square(x_ij);
}
+
if (s!=0.)
return std::sqrt(a)/std::sqrt(s);
return 0;
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);
const somenumber *src_ptr = src.begin();
for (unsigned int i=0; i<n; ++i, ++dst_ptr, ++src_ptr)
- *dst_ptr = om * *src_ptr / this->el(i,i);
+ *dst_ptr = somenumber(om) * *src_ptr / somenumber(this->el(i,i));
}
for (unsigned int i=0; i<m(); ++i)
{
for (unsigned int j=0; j<n(); ++j)
- if (std::fabs(this->el(i,j)) > threshold)
+ if (std::abs(this->el(i,j)) > threshold)
out << std::setw(width)
<< this->el(i,j) * denominator << ' ';
else
// regular
double diagonal_sum = 0;
for (unsigned int i=0; i<N; ++i)
- diagonal_sum += std::fabs(this->el(i,i));
+ diagonal_sum += std::abs(this->el(i,i));
const double typical_diagonal_element = diagonal_sum/N;
// initialize the array that holds
// 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; i<N; ++i)
{
- if (std::fabs(this->el(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;
}
}
const Vector<TYPERES>&) const;
+
+// #undef TYPEMAT
+
+// #define TYPEMAT std::complex<double>
+
+// template class FullMatrix<TYPEMAT>;
+
+// template void FullMatrix<TYPEMAT>::print(
+// LogStream&, const unsigned int, const unsigned int) const;
+// template void FullMatrix<TYPEMAT>::print(
+// std::ostream&, const unsigned int, const unsigned int) const;
+
+// template FullMatrix<TYPEMAT>& FullMatrix<TYPEMAT>::operator =(
+// const FullMatrix<std::complex<float> >&);
+
+// #undef TYPEMAT2
+// #define TYPEMAT2 std::complex<double>
+
+// template void FullMatrix<TYPEMAT>::fill<TYPEMAT2> (
+// const FullMatrix<TYPEMAT2>&, unsigned, unsigned, unsigned, unsigned);
+// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (
+// const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
+// template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (
+// const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
+// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+// const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::mmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
+// template void FullMatrix<TYPEMAT>::Tmmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
+// template void FullMatrix<TYPEMAT>::add_diag<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+// template void FullMatrix<TYPEMAT>::invert<TYPEMAT2> (const FullMatrix<TYPEMAT2>&);
+
+
+// #undef TYPEVEC
+// #undef TYPERES
+// #define TYPEVEC std::complex<double>
+// #define TYPERES std::complex<double>
+
+// template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
+// const FullMatrix<TYPEVEC>&,
+// const std::vector<unsigned int>&,
+// const std::vector<unsigned int>&);
+// template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+// template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+// template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
+// template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
+// const Vector<TYPEVEC> &) const;
+// template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
+// const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+// template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+// template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+
+// template
+// void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
+// Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
+
+// #undef TYPEVEC
+// #define TYPEVEC std::complex<float>
+
+// template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
+// const FullMatrix<TYPEVEC>&,
+// const std::vector<unsigned int>&,
+// const std::vector<unsigned int>&);
+// template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+// template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+// template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
+// template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
+// const Vector<TYPEVEC> &) const;
+// template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
+// const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+// template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+// template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
+// Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+// template
+// void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
+// Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
+
+
+// #undef TYPERES
+// #define TYPERES std::complex<float>
+
+// template
+// std::complex<double>
+// FullMatrix<TYPEMAT>::residual<TYPEVEC,TYPERES>(Vector<TYPEVEC>&,
+// const Vector<TYPEVEC>&,
+// const Vector<TYPERES>&) const;
+
+
DEAL_II_NAMESPACE_CLOSE