From 2ed864bda49bb303d6cd7b2401f2faa2d9d1175a Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Thu, 5 Nov 2026 16:27:49 +0100 Subject: [PATCH] numerical thresholds collected into a traits template --- conjgrad.h | 2 +- davidson.h | 4 ++-- gmres.h | 2 +- la_traits.h | 25 ++++++++++++++++++++++--- lanczos.h | 2 +- matexp.h | 4 ++-- nonclass.cc | 17 ++++++++--------- polynomial.h | 2 +- tensor.h | 2 +- 9 files changed, 39 insertions(+), 21 deletions(-) diff --git a/conjgrad.h b/conjgrad.h index 5775b20..b9b9031 100644 --- a/conjgrad.h +++ b/conjgrad.h @@ -34,7 +34,7 @@ namespace LA { template -extern bool conjgrad(const Matrix &bigmat, const NRVec &b, NRVec &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 &b, NRVec &x, const bool doguess=true, const double tol=LA_precision::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(); diff --git a/davidson.h b/davidson.h index 0283b61..41068a3 100644 --- a/davidson.h +++ b/davidson.h @@ -47,7 +47,7 @@ namespace LA { template void davidson(const Matrix &bigmat, NRVec &eivals, NRVec *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::normtype>::eigensolver_eps(), const bool incore=1, int maxit=100, const int maxkrylov = 500, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL) { @@ -235,7 +235,7 @@ eival_n = r[nroot]; for(i=0; i::normtype denom = LA_traits::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::normtype)0.1,std::abs(denom)):std::max((typename LA_traits::normtype)0.1,std::abs(denom)); vec1[i] /= denom; } } diff --git a/gmres.h b/gmres.h index 9a449cc..bd00077 100644 --- a/gmres.h +++ b/gmres.h @@ -73,7 +73,7 @@ for(int i=0;i<=k;i++) //x contains ev. initial guess and on return the solution template -bool gmres(const Matrix &bigmat, const NRVec &b, NRVec &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 &b, NRVec &x, const bool doguess=1, const double eps=LA_precision::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(); diff --git a/la_traits.h b/la_traits.h index 226b01a..207d6f1 100644 --- a/la_traits.h +++ b/la_traits.h @@ -457,9 +457,28 @@ template struct LA_traits : LA_traits_aux::scalar_type> {}; -//precision-dependent default thresholds (the double values are those used so far) -template struct LA_precision {static T svd_threshold() {return (T)1e-12;}}; -template<> struct LA_precision {static float svd_threshold() {return 1e-5f;}}; +//precision-dependent default thresholds; only float and double are supported (the double values are those used so far) +template struct LA_precision; +template<> struct LA_precision + { + 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 + { + 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 diff --git a/lanczos.h b/lanczos.h index f5e6529..6f974ed 100644 --- a/lanczos.h +++ b/lanczos.h @@ -39,7 +39,7 @@ namespace LA { template extern void lanczos(const Matrix &bigmat, NRVec &eivals, NRVec *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::normtype>::eigensolver_eps(), const bool incore=1, int maxit=100, const int maxkrylov = 1000, void (*initguess)(NRVec &)=NULL, const typename LA_traits::normtype *target=NULL) { diff --git a/matexp.h b/matexp.h index 85aa814..b73b2b1 100644 --- a/matexp.h +++ b/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::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::normtype>::taylor_precision(); //further decreasing brings nothing double s,t; s=mnorm*scale; int n=0; diff --git a/nonclass.cc b/nonclass.cc index ade946a..7e796a1 100644 --- a/nonclass.cc +++ b/nonclass.cc @@ -50,7 +50,6 @@ INSTANTIZE(unsigned int) INSTANTIZE(unsigned long) INSTANTIZE(unsigned long long) -#define EPSDET 1e-300 template void lawritemat(FILE *file,const T *a,int r,int c,const char *form0, @@ -168,7 +167,7 @@ static void linear_solve_do(NRMat &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::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 &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::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 &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::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 &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::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 > &a, std::complex::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 > &a, std::complex 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::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 > &A, std::complex 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 t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < LA_precision::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 > &A, std::complex if (det && r==0) { *det = 1.; //take into account some numerical instabilities in zgesv for singular matrices - for (int i=0; i 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 t=A[i][i]; if(!std::isfinite(t.real()) || !std::isfinite(t.imag()) || std::abs(t) < LA_precision::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 companion() const; //matrix which has this characteristic polynomial NRVec::complextype> roots() const; //implemented for complex and double only NRVec realroots(const typename LA_traits::normtype thr) const; - T newton(const T x0, const typename LA_traits::normtype thr=1e-14, const int maxit=1000) const; //solve root from the guess + T newton(const T x0, const typename LA_traits::normtype thr=LA_precision::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 diff --git a/tensor.h b/tensor.h index dc3bdd7..9b75be2 100644 --- a/tensor.h +++ b/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 &nl, int symmetry=0) const {return merge_indices(findindexlist(nl),symmetry);}; - NRVec > Tucker(typename LA_traits::normtype thr=1e-12, bool inverseorder=false); //HOSVD-Tucker decomposition, return core tensor in *this, flattened + NRVec > Tucker(typename LA_traits::normtype thr=LA_precision::normtype>::svd_threshold(), bool inverseorder=false); //HOSVD-Tucker decomposition, return core tensor in *this, flattened Tensor inverseTucker(const NRVec > &x, bool inverseorder=false) const; //rebuild the original tensor from Tucker Tensor linear_transform(const NRVec > &x) const; //linear transform by a different matrix per each index group, preserving group symmetries Tensor linear_transform(const int g, const Tensor &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)