2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
5 #include "gromacs/utility/real.h"
8 F77_FUNC(dlasq6,DLASQ6)(int *i0,
26 const double safemin = GMX_DOUBLE_MIN*(1.0+GMX_DOUBLE_EPS);
30 if (*n0 - *i0 - 1 <= 0) {
34 j4 = (*i0 << 2) + *pp - 3;
41 for (j4 = *i0*4; j4 <= i__1; j4 += 4) {
42 z__[j4 - 2] = d__ + z__[j4 - 1];
43 if (fabs(z__[j4 - 2])<GMX_DOUBLE_MIN) {
48 } else if (safemin * z__[j4 + 1] < z__[j4 - 2] && safemin * z__[j4
50 temp = z__[j4 + 1] / z__[j4 - 2];
51 z__[j4] = z__[j4 - 1] * temp;
54 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
55 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
60 d__1 = emin, d__2 = z__[j4];
61 emin = (d__1<d__2) ? d__1 : d__2;
65 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
66 z__[j4 - 3] = d__ + z__[j4];
67 if (fabs(z__[j4 - 3])<GMX_DOUBLE_MIN) {
72 } else if (safemin * z__[j4 + 2] < z__[j4 - 3] && safemin * z__[j4
74 temp = z__[j4 + 2] / z__[j4 - 3];
75 z__[j4 - 1] = z__[j4] * temp;
78 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
79 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
83 d__1 = emin, d__2 = z__[j4 - 1];
84 emin = (d__1<d__2) ? d__1 : d__2;
90 j4 = 4*(*n0 - 2) - *pp;
91 j4p2 = j4 + (*pp << 1) - 1;
92 z__[j4 - 2] = *dnm2 + z__[j4p2];
93 if (fabs(z__[j4 - 2])<GMX_DOUBLE_MIN) {
95 *dnm1 = z__[j4p2 + 2];
98 } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
100 temp = z__[j4p2 + 2] / z__[j4 - 2];
101 z__[j4] = z__[j4p2] * temp;
102 *dnm1 = *dnm2 * temp;
104 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
105 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
112 j4p2 = j4 + (*pp << 1) - 1;
113 z__[j4 - 2] = *dnm1 + z__[j4p2];
114 if (fabs(z__[j4 - 2])<GMX_DOUBLE_MIN) {
119 } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
121 temp = z__[j4p2 + 2] / z__[j4 - 2];
122 z__[j4] = z__[j4p2] * temp;
125 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
126 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
132 z__[(*n0 << 2) - *pp] = emin;