numerical thresholds collected into a traits template

conjgrad.h

@@ -34,7 +34,7 @@ namespace LA {
 
 
 template<typename T, typename Matrix>
-extern bool conjgrad(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=true, const double tol=1e-8, const int itmax=1000, const bool verbose=true, bool issquare=1,const bool precondition=1)
+extern bool conjgrad(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=true, const double tol=LA_precision<typename LA_traits<T>::normtype>::conjgrad_tol(), const int itmax=1000, const bool verbose=true, bool issquare=1,const bool precondition=1)
 {
 int m=bigmat.nrows();
 int n=bigmat.ncols();

davidson.h

@@ -47,7 +47,7 @@ namespace LA {
 
 template <typename T, typename Matrix>
 void davidson(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
-		int nroots=1,  const bool verbose=0, const double eps=1e-6,
+		int nroots=1,  const bool verbose=0, const double eps=LA_precision<typename LA_traits<T>::normtype>::eigensolver_eps(),
 	 	const bool incore=1, int maxit=100, const int maxkrylov = 500,
 		void (*initguess)(NRVec<T> &)=NULL, const typename LA_traits<T>::normtype *target=NULL)
 {
@@ -235,7 +235,7 @@ eival_n = r[nroot];
 for(i=0; i<n; ++i)
 	{
 	typename LA_traits<T>::normtype denom = LA_traits<T>::realpart(diagonal[i]) - eival_n;
-        denom = denom<0?-std::max(0.1,std::abs(denom)):std::max(0.1,std::abs(denom));
+        denom = denom<0?-std::max((typename LA_traits<T>::normtype)0.1,std::abs(denom)):std::max((typename LA_traits<T>::normtype)0.1,std::abs(denom));
         vec1[i] /= denom;
 	}
 }

gmres.h

@@ -73,7 +73,7 @@ for(int i=0;i<=k;i++)
 
 //x contains ev. initial guess and on return the solution
 template<typename T, typename Matrix>
-bool gmres(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=1, const double eps=1e-7, const int MAXIT=50, const bool verbose=1, bool square=1,const bool precondition=1, int neustart=0, const int incore=1)
+bool gmres(const Matrix &bigmat, const NRVec<T> &b, NRVec<T> &x, const bool doguess=1, const double eps=LA_precision<typename LA_traits<T>::normtype>::gmres_eps(), const int MAXIT=50, const bool verbose=1, bool square=1,const bool precondition=1, int neustart=0, const int incore=1)
 {
 int zeilen=bigmat.nrows();
 int spalten=bigmat.ncols();

la_traits.h

@@ -457,9 +457,28 @@ template<typename C>
 struct LA_traits : LA_traits_aux<C, typename isscalar<C>::scalar_type> {};
 
 
-//precision-dependent default thresholds (the double values are those used so far)
+//precision-dependent default thresholds; only float and double are supported (the double values are those used so far)
-template<typename T> struct LA_precision {static T svd_threshold() {return (T)1e-12;}};
+template<typename T> struct LA_precision;
-template<> struct LA_precision<float> {static float svd_threshold() {return 1e-5f;}};
+template<> struct LA_precision<double>
+	{
+	static double svd_threshold() {return 1e-12;} //singular values treated as zero (svd_solve, Tucker)
+	static double taylor_precision() {return 1e-14;} //truncation of the Taylor series in matexp
+	static double conjgrad_tol() {return 1e-8;} //default convergence threshold of conjgrad
+	static double gmres_eps() {return 1e-7;} //default relative residual of gmres
+	static double eigensolver_eps() {return 1e-6;} //default convergence threshold of davidson and lanczos
+	static double newton_threshold() {return 1e-14;} //convergence of Polynomial::newton
+	static double det_tiny() {return 1e-300;} //determinant underflow guard in linear_solve
+	};
+template<> struct LA_precision<float>
+	{
+	static float svd_threshold() {return 1e-5f;}
+	static float taylor_precision() {return 1e-7f;}
+	static float conjgrad_tol() {return 1e-5f;}
+	static float gmres_eps() {return 1e-5f;}
+	static float eigensolver_eps() {return 1e-4f;}
+	static float newton_threshold() {return 1e-6f;}
+	static float det_tiny() {return 1e-37f;}
+	};
 
 }//namespace
 

lanczos.h

@@ -39,7 +39,7 @@ namespace LA {
 
 template <typename T, typename Matrix>
 extern void lanczos(const Matrix &bigmat, NRVec<T> &eivals, NRVec<T> *eivecs, const char *eivecsfile, 
-		int nroots=1,  const bool verbose=0, const double eps=1e-6,
+		int nroots=1,  const bool verbose=0, const double eps=LA_precision<typename LA_traits<T>::normtype>::eigensolver_eps(),
 	 	const bool incore=1, int maxit=100, const int maxkrylov = 1000,
 		void (*initguess)(NRVec<T> &)=NULL, const typename LA_traits<T>::normtype *target=NULL)
 {

matexp.h

@@ -274,7 +274,7 @@ double scale=std::exp(-std::log(2.)*power);
 
 
 //find how long taylor expansion will be necessary
-const double precision=1e-14; //further decreasing brings nothing
+const double precision=LA_precision<typename LA_traits<C>::normtype>::taylor_precision(); //further decreasing brings nothing
 double s,t;
 s=mnorm*scale;
 int n=0;
@@ -311,7 +311,7 @@ if(maxpower>=0 && power>maxpower) power=maxpower;
 double scale=std::exp(-std::log(2.)*power);
 
 //find how long taylor expansion will be necessary
-const double precision=1e-14; //further decreasing brings nothing
+const double precision=LA_precision<typename LA_traits<C>::normtype>::taylor_precision(); //further decreasing brings nothing
 double s,t;
 s=mnorm*scale;
 int n=0;

nonclass.cc

@@ -50,7 +50,6 @@ INSTANTIZE(unsigned int)
 INSTANTIZE(unsigned long)
 INSTANTIZE(unsigned long long)
 
-#define EPSDET 1e-300
 
 template <typename T>
 void lawritemat(FILE *file,const T *a,int r,int c,const char *form0,
@@ -168,7 +167,7 @@ static void linear_solve_do(NRMat<double> &A, double *B, const int nrhs, const i
 	if (det && r==0) {
 		*det = 1.;
 		//take into account some numerical instabilities in dgesv for singular matrices
-		for (int i=0; i<n; ++i)  {double t=A[i][i]; if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
+		for (int i=0; i<n; ++i)  {double t=A[i][i]; if(!finite(t) || std::abs(t) < LA_precision<double>::det_tiny() ) {*det=0.; break;} else *det *=t;}
 		//find out whether ipiv are numbered from 0 or from 1
 		int shift=1;
 	        for (int i=0; i<n; ++i)
@@ -224,7 +223,7 @@ static void linear_solve_do(NRMat<float> &A, float *B, const int nrhs, const int
 	if (det && r==0) {
 		*det = 1.;
 		//take into account some numerical instabilities in dgesv for singular matrices
-		for (int i=0; i<n; ++i)  {float t=A[i][i]; if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
+		for (int i=0; i<n; ++i)  {float t=A[i][i]; if(!finite(t) || std::abs(t) < LA_precision<float>::det_tiny() ) {*det=0.; break;} else *det *=t;}
 		//find out whether ipiv are numbered from 0 or from 1
 		int shift=1;
 	        for (int i=0; i<n; ++i)
@@ -319,7 +318,7 @@ static void linear_solve_do(NRSMat<double> &a, double *b, const int nrhs, const
 			double t;
 			if(ipiv[i]<0) {t = a(i,i)*a(i+1,i+1)-a(i,i+1)*a(i,i+1); ++i;} //2x2 block of D (dsptrf with UPLO='U')
 			else t=a(i,i);
-			if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;
+			if(!finite(t) || std::abs(t) < LA_precision<double>::det_tiny() ) {*det=0.; break;} else *det *= t;
 			}
 		//ipiv is needed only for the block structure of D, the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
 	}
@@ -355,7 +354,7 @@ static void linear_solve_do(NRSMat<float> &a, float *b, const int nrhs, const in
 			float t;
 			if(ipiv[i]<0) {t = a(i,i)*a(i+1,i+1)-a(i,i+1)*a(i,i+1); ++i;} //2x2 block of D (dsptrf with UPLO='U')
 			else t=a(i,i);
-			if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;
+			if(!finite(t) || std::abs(t) < LA_precision<float>::det_tiny() ) {*det=0.; break;} else *det *= t;
 			}
 		//ipiv is needed only for the block structure of D, the permutation matrix occurs twice in the decomposition and signs thus cancel (man dspsv)
 	}
@@ -435,7 +434,7 @@ static void linear_solve_do(NRSMat<std::complex<double> > &a, std::complex<doubl
 			double t; //determinant of a hermitian matrix is real, det(conj(A)) = det(A)
 			if(ipiv[i]<0) {t = a(i,i).real()*a(i+1,i+1).real()-std::norm(a(i,i+1)); ++i;} //2x2 block of D (zhptrf with UPLO='U')
 			else t=a(i,i).real();
-			if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;
+			if(!finite(t) || std::abs(t) < LA_precision<double>::det_tiny() ) {*det=0.; break;} else *det *= t;
 			}
 		//ipiv is needed only for the block structure of D, the permutation matrix occurs twice in the decomposition and signs thus cancel (man zhpsv)
 	}
@@ -470,7 +469,7 @@ static void linear_solve_do(NRSMat<std::complex<float> > &a, std::complex<float>
 			float t; //determinant of a hermitian matrix is real, det(conj(A)) = det(A)
 			if(ipiv[i]<0) {t = a(i,i).real()*a(i+1,i+1).real()-std::norm(a(i,i+1)); ++i;} //2x2 block of D (zhptrf with UPLO='U')
 			else t=a(i,i).real();
-			if(!finite(t) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *= t;
+			if(!finite(t) || std::abs(t) < LA_precision<float>::det_tiny() ) {*det=0.; break;} else *det *= t;
 			}
 		//ipiv is needed only for the block structure of D, the permutation matrix occurs twice in the decomposition and signs thus cancel (man zhpsv)
 	}
@@ -551,7 +550,7 @@ static void linear_solve_do(NRMat<std::complex<double> > &A, std::complex<double
 	if (det && r==0) {
 		*det = 1.;
 		//take into account some numerical instabilities in zgesv for singular matrices
-		for (int i=0; i<n; ++i)  {std::complex<double> t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
+		for (int i=0; i<n; ++i)  {std::complex<double> t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < LA_precision<double>::det_tiny() ) {*det=0.; break;} else *det *=t;}
 		//find out whether ipiv are numbered from 0 or from 1
 		int shift=1;
 	        for (int i=0; i<n; ++i)
@@ -603,7 +602,7 @@ static void linear_solve_do(NRMat<std::complex<float> > &A, std::complex<float>
 	if (det && r==0) {
 		*det = 1.;
 		//take into account some numerical instabilities in zgesv for singular matrices
-		for (int i=0; i<n; ++i)  {std::complex<float> t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < EPSDET ) {*det=0.; break;} else *det *=t;}
+		for (int i=0; i<n; ++i)  {std::complex<float> t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < LA_precision<float>::det_tiny() ) {*det=0.; break;} else *det *=t;}
 		//find out whether ipiv are numbered from 0 or from 1
 		int shift=1;
 	        for (int i=0; i<n; ++i)

polynomial.h

@@ -156,7 +156,7 @@ public:
 	NRMat<T> companion() const; //matrix which has this characteristic polynomial
 	NRVec<typename LA_traits<T>::complextype> roots() const; //implemented for complex<double> and double only
 	NRVec<T> realroots(const typename LA_traits<T>::normtype thr) const;
-	T newton(const T x0, const typename LA_traits<T>::normtype thr=1e-14, const int maxit=1000) const; //solve root from the guess
+	T newton(const T x0, const typename LA_traits<T>::normtype thr=LA_precision<typename LA_traits<T>::normtype>::newton_threshold(), const int maxit=1000) const; //solve root from the guess
 	Polynomial pow(const int i) const; //integer power
 	Polynomial powx(const int i) const; //substitute x^i for x in the polynomial
 	void binomial(const int n); //(1+x)^n

tensor.h

@@ -446,7 +446,7 @@ public:
 	Tensor merge_indices(const INDEXLIST &il, int symmetry=0) const;	//opposite to flatten (merging with optional symmetrization/antisymmetrization and compression)
 	Tensor merge_indices(const NRVec<INDEXNAME> &nl, int symmetry=0) const {return merge_indices(findindexlist(nl),symmetry);};
 
-	NRVec<NRMat<T> > Tucker(typename LA_traits<T>::normtype thr=1e-12, bool inverseorder=false); //HOSVD-Tucker decomposition, return core tensor in *this, flattened 
+	NRVec<NRMat<T> > Tucker(typename LA_traits<T>::normtype thr=LA_precision<typename LA_traits<T>::normtype>::svd_threshold(), bool inverseorder=false); //HOSVD-Tucker decomposition, return core tensor in *this, flattened 
 	Tensor inverseTucker(const NRVec<NRMat<T> > &x, bool inverseorder=false) const; //rebuild the original tensor from Tucker
 	Tensor linear_transform(const NRVec<NRMat<T> > &x) const; //linear transform by a different matrix per each index group, preserving group symmetries
 	Tensor linear_transform(const int g, const Tensor<T> &x) const; //linear transform a single group, preserve its position and symmetry, x must be rank-2 flat tensor (this is useful for raising/lowering indices)