From 52b504c9137879ea1696b8fab17c6e38f6c9eea8 Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@dune-project.org>
Date: Fri, 12 Nov 2004 15:52:13 +0000
Subject: [PATCH] introcued ctype into ISTLPrecision but currently uses default
ctype, so not completely usefull. We first need the ISTL traits ... otherwise
we can't have usefull values for ctype=complex<...>.
[[Imported from SVN: r1068]]
---
istl/fmatrix.hh | 16 ++++++++--------
istl/precision.hh | 31 +++++++++++++++++--------------
2 files changed, 25 insertions(+), 22 deletions(-)
diff --git a/istl/fmatrix.hh b/istl/fmatrix.hh
index 59c3b1adc..ccf474604 100644
--- a/istl/fmatrix.hh
+++ b/istl/fmatrix.hh
@@ -245,8 +245,8 @@ namespace Dune {
// Gaussian elimination with maximum column pivot
double norm=A.infinity_norm_real(); // for relative thresholds
- double pivthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::pivoting_limit());
- double singthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::singular_limit());
+ double pivthres = std::max(ISTLPrecision<>::absolute_limit(),norm*ISTLPrecision<>::pivoting_limit());
+ double singthres = std::max(ISTLPrecision<>::absolute_limit(),norm*ISTLPrecision<>::singular_limit());
V& rhs = x; // use x to store rhs
rhs = b; // copy data
@@ -299,7 +299,7 @@ namespace Dune {
inline void fm_solve (const FieldMatrix<K,1,1>& A, V& x, const V& b)
{
#ifdef DUNE_ISTL_WITH_CHECKING
- if (fvmeta_absreal(A[0][0])<ISTLPrecision::absolute_limit())
+ if (fvmeta_absreal(A[0][0])<ISTLPrecision<>::absolute_limit())
DUNE_THROW(ISTLError,"matrix is singular");
#endif
x[0] = b[0]/A[0][0];
@@ -311,7 +311,7 @@ namespace Dune {
{
#ifdef DUNE_ISTL_WITH_CHECKING
K detinv = A[0][0]*A[1][1]-A[0][1]*A[1][0];
- if (fvmeta_absreal(detinv)<ISTLPrecision::absolute_limit())
+ if (fvmeta_absreal(detinv)<ISTLPrecision<>::absolute_limit())
DUNE_THROW(ISTLError,"matrix is singular");
detinv = 1/detinv;
#else
@@ -334,8 +334,8 @@ namespace Dune {
FieldMatrix<K,n,n>& U=A;
double norm=A.infinity_norm_real(); // for relative thresholds
- double pivthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::pivoting_limit());
- double singthres = std::max(ISTLPrecision::absolute_limit(),norm*ISTLPrecision::singular_limit());
+ double pivthres = std::max(ISTLPrecision<>::absolute_limit(),norm*ISTLPrecision<>::pivoting_limit());
+ double singthres = std::max(ISTLPrecision<>::absolute_limit(),norm*ISTLPrecision<>::singular_limit());
// LU decomposition of A in A
for (int i=0; i<n; i++) // loop over all rows
@@ -397,7 +397,7 @@ namespace Dune {
void fm_invert (FieldMatrix<K,1,1>& A)
{
#ifdef DUNE_ISTL_WITH_CHECKING
- if (fvmeta_absreal(A[0][0])<ISTLPrecision::absolute_limit())
+ if (fvmeta_absreal(A[0][0])<ISTLPrecision<>::absolute_limit())
DUNE_THROW(ISTLError,"matrix is singular");
#endif
A[0][0] = 1/A[0][0];
@@ -409,7 +409,7 @@ namespace Dune {
{
K detinv = A[0][0]*A[1][1]-A[0][1]*A[1][0];
#ifdef DUNE_ISTL_WITH_CHECKING
- if (fvmeta_absreal(detinv)<ISTLPrecision::absolute_limit())
+ if (fvmeta_absreal(detinv)<ISTLPrecision<>::absolute_limit())
DUNE_THROW(ISTLError,"matrix is singular");
#endif
detinv = 1/detinv;
diff --git a/istl/precision.hh b/istl/precision.hh
index 75da98365..6c936d859 100644
--- a/istl/precision.hh
+++ b/istl/precision.hh
@@ -12,55 +12,58 @@ namespace Dune {
@{
*/
- class ISTLPrecision { // uses standard malloc and free
+ template <class ctype = double>
+ class ISTLPrecision {
public:
//! return threshold to do pivoting
- static double pivoting_limit ()
+ static ctype pivoting_limit ()
{
return _pivoting;
}
//! set pivoting threshold
- static void set_pivoting_limit (double pivthres)
+ static void set_pivoting_limit (ctype pivthres)
{
_pivoting = pivthres;
}
//! return threshold to declare matrix singular
- static double singular_limit ()
+ static ctype singular_limit ()
{
return _singular;
}
//! set singular threshold
- static void set_singular_limit (double singthres)
+ static void set_singular_limit (ctype singthres)
{
_singular = singthres;
}
//! return threshold to declare matrix singular
- static double absolute_limit ()
+ static ctype absolute_limit ()
{
return _absolute;
}
//! set singular threshold
- static void set_absolute_limit (double absthres)
+ static void set_absolute_limit (ctype absthres)
{
_absolute = absthres;
}
private:
// just to demonstrate some state information
- static double _pivoting;
- static double _singular;
- static double _absolute;
+ static ctype _pivoting;
+ static ctype _singular;
+ static ctype _absolute;
};
- double ISTLPrecision::_pivoting = 1E-8;
- double ISTLPrecision::_singular = 1E-14;
- double ISTLPrecision::_absolute = 1E-80;
-
+ template <class ctype>
+ ctype ISTLPrecision<ctype>::_pivoting = 1E-8;
+ template <class ctype>
+ ctype ISTLPrecision<ctype>::_singular = 1E-14;
+ template <class ctype>
+ ctype ISTLPrecision<ctype>::_absolute = 1E-80;
/** @} end documentation */
--
GitLab