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

@@ -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;

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 */