diff --git a/ChangeLog b/ChangeLog index f5e9d329..22ab44ca 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,5 @@ Version 5.2.0 (2023-01-31): - * Optimize SR-ERI integral screen + * Optimize SR-ERI integral screening Version 5.1.9 (2023-01-09): * Fix approx_log bug for non-x86 architecture diff --git a/src/rys_roots.c b/src/rys_roots.c index a4c87885..c6fb9e39 100644 --- a/src/rys_roots.c +++ b/src/rys_roots.c @@ -24,11 +24,11 @@ #define PIE4 0.78539816339744827900 #define THRESHOLD_ZERO (DBL_EPSILON * 8) -static void rys_root1(double x, double *roots, double *weights); -static void rys_root2(double x, double *roots, double *weights); -static void rys_root3(double x, double *roots, double *weights); -static void rys_root4(double x, double *roots, double *weights); -static void rys_root5(double x, double *roots, double *weights); +static int rys_root1(double x, double *roots, double *weights); +static int rys_root2(double x, double *roots, double *weights); +static int rys_root3(double x, double *roots, double *weights); +static int rys_root4(double x, double *roots, double *weights); +static int rys_root5(double x, double *roots, double *weights); typedef int QuadratureFunction(int n, double x, double lower, double *roots, double *weights); #ifndef HAVE_QUADMATH_H #define CINTqrys_schmidt CINTlrys_schmidt @@ -48,12 +48,6 @@ static int segment_solve(int n, double x, double lower, double *u, double *w, if (error) { error = CINTqrys_schmidt(n, x, lower, u, w); } - if (error) { - fprintf(stderr, "libcint rys_roots failed. nroots=%d\n", n); -#ifndef KEEP_GOING - exit(1); -#endif - } return error; } @@ -69,41 +63,49 @@ void CINTrys_roots(int nroots, double x, double *u, double *w) return; } + int err; switch (nroots) { case 1: - rys_root1(x, u, w); + err = rys_root1(x, u, w); break; case 2: - rys_root2(x, u, w); + err = rys_root2(x, u, w); break; case 3: - rys_root3(x, u, w); + err = rys_root3(x, u, w); break; case 4: - rys_root4(x, u, w); + err = rys_root4(x, u, w); break; case 5: - rys_root5(x, u, w); + err = rys_root5(x, u, w); break; case 6: case 7: - segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTrys_schmidt); + err = segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTrys_schmidt); break; case 8: - segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTlrys_schmidt); + err = segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTlrys_schmidt); break; case 9: - segment_solve(nroots, x, 0., u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, 0., u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); break; case 10: case 11: - segment_solve(nroots, x, 0., u, w, 18, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, 0., u, w, 18, CINTlrys_jacobi, CINTlrys_laguerre); break; case 12: - segment_solve(nroots, x, 0., u, w, 22, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, 0., u, w, 22, CINTlrys_jacobi, CINTlrys_laguerre); break; default: - segment_solve(nroots, x, 0., u, w, 50, CINTqrys_jacobi, CINTqrys_laguerre); + err = segment_solve(nroots, x, 0., u, w, 50, CINTqrys_jacobi, CINTqrys_laguerre); break; } + if (err) { + fprintf(stderr, "rys_roots fails: nroots=%d x=%g\n", + nroots, x); +#ifndef KEEP_GOING + exit(err); +#endif + } } #ifdef WITH_RANGE_COULOMB @@ -124,89 +126,79 @@ static int segment_solve1(int n, double x, double lower, double *u, double *w, } else if (lower < lower_bp2) { error = fn3(n, x, lower, u, w); } else { - fprintf(stderr, "libcint SR-rys_roots does not support nroots=%d x=%g lower=%g\n", - n, x, lower); -#ifndef KEEP_GOING - exit(1); -#endif return 1; } if (error) { error = CINTqrys_schmidt(n, x, lower, u, w); - if (error) { - fprintf(stderr, "libcint SR-rys_roots failed. nroots=%d\n", n); -#ifndef KEEP_GOING - exit(1); -#endif - } } return error; } void CINTsr_rys_roots(int nroots, double x, double lower, double *u, double *w) { + int err; switch (nroots) { case 1: - CINTrys_schmidt(nroots, x, lower, u, w); + err = CINTrys_schmidt(nroots, x, lower, u, w); break; case 2: if (lower < 0.99) { - CINTrys_schmidt(nroots, x, lower, u, w); + err = CINTrys_schmidt(nroots, x, lower, u, w); } else { - CINTqrys_jacobi(nroots, x, lower, u, w); + err = CINTqrys_jacobi(nroots, x, lower, u, w); } break; case 3: if (lower < 0.93) { - CINTrys_schmidt(nroots, x, lower, u, w); + err = CINTrys_schmidt(nroots, x, lower, u, w); } else if (lower < 0.97) { - segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); } else { - CINTqrys_jacobi(nroots, x, lower, u, w); + err = CINTqrys_jacobi(nroots, x, lower, u, w); } break; case 4: - if (lower < 0.85) { - CINTrys_schmidt(nroots, x, lower, u, w); + if (lower < 0.8) { + err = CINTrys_schmidt(nroots, x, lower, u, w); } else if (lower < 0.9) { - segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); } else { - CINTqrys_jacobi(nroots, x, lower, u, w); + err = CINTqrys_jacobi(nroots, x, lower, u, w); } break; case 5: - if (lower < 0.45) { - CINTrys_schmidt(nroots, x, lower, u, w); + if (lower < 0.4) { + err = segment_solve(nroots, x, lower, u, w, 50, CINTrys_schmidt, CINTlrys_laguerre); } else if (lower < 0.8) { - segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); } else { - CINTqrys_jacobi(nroots, x, lower, u, w); + err = CINTqrys_jacobi(nroots, x, lower, u, w); } break; case 6: - if (lower < 0.35) { - CINTrys_schmidt(nroots, x, lower, u, w); + if (lower < 0.25) { + err = segment_solve(nroots, x, lower, u, w, 60, CINTrys_schmidt, CINTlrys_laguerre); } else if (lower < 0.8) { - segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); + err = segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre); } else { - CINTqrys_jacobi(nroots, x, lower, u, w); + err = CINTqrys_jacobi(nroots, x, lower, u, w); } break; case 7: - segment_solve1(nroots, x, lower, u, w, 0.55, 1., 60, CINTlrys_jacobi, CINTlrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.5, 1., 60, CINTlrys_jacobi, CINTlrys_laguerre, CINTqrys_jacobi); break; case 8: case 9: case 10: //CINTqrys_jacobi(nroots, x, lower, u, w); - segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 11: case 12: - segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 13: case 14: - segment_solve1(nroots, x, lower, u, w, 0.25, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.25, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 15: case 16: - segment_solve1(nroots, x, lower, u, w, 0.25, 0.75, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.25, 0.75, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 17: segment_solve1(nroots, x, lower, u, w, 0.25, 0.65, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); @@ -215,35 +207,42 @@ void CINTsr_rys_roots(int nroots, double x, double lower, double *u, double *w) segment_solve1(nroots, x, lower, u, w, 0.15, 0.65, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 19: - segment_solve1(nroots, x, lower, u, w, 0.15, 0.55, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.15, 0.55, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 20: case 21: - segment_solve1(nroots, x, lower, u, w, 0.25, 0.45, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.25, 0.45, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; case 22: case 23: case 24: - segment_solve1(nroots, x, lower, u, w, 0.25, 0.35, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); + err = segment_solve1(nroots, x, lower, u, w, 0.25, 0.35, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi); break; default: fprintf(stderr, "libcint SR-rys_roots does not support nroots=%d\n", nroots); #ifndef KEEP_GOING exit(1); +#endif + } + if (err) { + fprintf(stderr, "sr_rys_roots fails: nroots=%d x=%g lower=%g\n", + nroots, x, lower); +#ifndef KEEP_GOING + exit(err); #endif } } #endif -static void rys_root1(double X, double *roots, double *weights) +static int rys_root1(double X, double *roots, double *weights) { double Y, F1; if (X > 33.) { weights[0] = sqrt(PIE4/X); roots[0] = 0.5E+00/(X-0.5E+00); - return; + return 0; } else if (X < 3.e-7) { weights[0] = 1.0E+00 -X/3.0E+00; roots[0] = 0.5E+00 -X/5.0E+00; - return; + return 0; } double E = exp(-X); @@ -292,9 +291,10 @@ static void rys_root1(double X, double *roots, double *weights) double WW1 = 2. * X * F1 + E; weights[0] = WW1; roots[0] = F1 / (WW1 - F1); + return 0; } -static void rys_root2(double X, double *roots, double *weights) +static int rys_root2(double X, double *roots, double *weights) { double R12, R22, W22; @@ -450,9 +450,10 @@ static void rys_root2(double X, double *roots, double *weights) roots[1] = RT2; weights[0] = WW1; weights[1] = WW2; + return 0; } -static void rys_root3(double X, double *roots, double *weights) +static int rys_root3(double X, double *roots, double *weights) { double R13, R23, W23, R33, W33; @@ -720,9 +721,10 @@ static void rys_root3(double X, double *roots, double *weights) weights[0] = WW1; weights[1] = WW2; weights[2] = WW3; + return 0; } -static void rys_root4(double X, double *roots, double *weights) +static int rys_root4(double X, double *roots, double *weights) { double R14, R24, W24, R34, W34, R44, W44; double RT1, RT2, RT3, RT4, WW1, WW2, WW3, WW4; @@ -1085,9 +1087,10 @@ static void rys_root4(double X, double *roots, double *weights) weights[1] = WW2; weights[2] = WW3; weights[3] = WW4; + return 0; } -static void rys_root5(double X, double *roots, double *weights) +static int rys_root5(double X, double *roots, double *weights) { double R15,R25,W25,R35,W35,R45,W45,R55,W55; double RT1, RT2, RT3, RT4, RT5, WW1, WW2, WW3, WW4, WW5; @@ -1588,6 +1591,7 @@ static void rys_root5(double X, double *roots, double *weights) weights[2] = WW3; weights[3] = WW4; weights[4] = WW5; + return 0; } #define POLYNOMIAL_VALUE1(p, a, order, x) \ @@ -1617,7 +1621,8 @@ static int R_dnode(double *a, double *roots, int order) continue; } if (p0 * p1init > 0) { - fprintf(stderr, "ROOT NUMBER %d WAS NOT FOUND FOR POLYNOMIAL OF ORDER %d\n", m, order); + fprintf(stderr, "ROOT NUMBER %d WAS NOT FOUND FOR POLYNOMIAL OF ORDER %d\n", + m, order); return 1; } if (x0 <= x1init) { @@ -1761,15 +1766,9 @@ static int _rdk_rys_roots(int nroots, double *fmt_ints, } int error = R_dsmit(cs, fmt_ints, nroots1); -#ifndef KEEP_GOING if (error) { - exit(error); - } -#else - if (error == 1) { return 1; } -#endif dum = sqrt(cs[2*nroots1+1] * cs[2*nroots1+1] - 4 * cs[2*nroots1+0] * cs[2*nroots1+2]); rt[0] = .5 * (-cs[2*nroots1+1] - dum) / cs[2*nroots1+2]; @@ -1783,11 +1782,7 @@ static int _rdk_rys_roots(int nroots, double *fmt_ints, a = cs + order * nroots1; error = R_dnode(a, rt, order); if (error) { -#ifndef KEEP_GOING - exit(error); -#else return error; -#endif } } @@ -1936,15 +1931,9 @@ int CINTlrys_schmidt(int nroots, double x, double lower, double *roots, double * rt[0] = fmt_ints[1] / fmt_ints[0]; } else { error = R_lsmit(qcs, fmt_ints, nroots1); -#ifndef KEEP_GOING if (error) { - exit(error); - } -#else - if (error == 1) { - return 1; + return error; } -#endif for (k = 1; k < nroots1; k++) { for (i = 0; i <= k; i++) { cs[k * nroots1 + i] = qcs[k * nroots1 + i]; @@ -1963,11 +1952,7 @@ int CINTlrys_schmidt(int nroots, double x, double lower, double *roots, double * a = cs + order * nroots1; error = R_dnode(a, rt, order); if (error) { -#ifndef KEEP_GOING - exit(error); -#else return error; -#endif } } @@ -2073,15 +2058,9 @@ int CINTqrys_schmidt(int nroots, double x, double lower, double *roots, double * rt[0] = fmt_ints[1] / fmt_ints[0]; } else { error = R_qsmit(qcs, fmt_ints, nroots1); -#ifndef KEEP_GOING if (error) { - exit(error); - } -#else - if (error == 1) { - return 1; + return error; } -#endif for (k = 1; k < nroots1; k++) { for (i = 0; i <= k; i++) { cs[k * nroots1 + i] = qcs[k * nroots1 + i]; @@ -2100,11 +2079,7 @@ int CINTqrys_schmidt(int nroots, double x, double lower, double *roots, double * a = cs + order * nroots1; error = R_dnode(a, rt, order); if (error) { -#ifndef KEEP_GOING - exit(error); -#else return error; -#endif } }