2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
5 #include "gromacs/utility/real.h"
8 F77_FUNC(ssterf,SSTERF)(int *n,
20 float bb, rt1, rt2, eps, rte;
24 float gamma, alpha, sigma, anorm;
35 const float safmin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
55 ssfmax = sqrt(safmax) / 3.;
56 ssfmin = sqrt(safmin) / eps2;
66 F77_FUNC(slasrt,SLASRT)("I", n, &d__[1], info);
73 for (m = l1; m <= i__1; ++m) {
74 if (fabs(e[m]) <= sqrt(fabs(d__[m])) *
75 sqrt(fabs(d__[m + 1])) * eps) {
93 anorm = F77_FUNC(slanst,SLANST)("I", &i__1, &d__[l], &e[l]);
98 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
101 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
103 } else if (anorm < ssfmin) {
106 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
109 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
114 for (i__ = l; i__ <= i__1; ++i__) {
116 e[i__] = d__1 * d__1;
119 if (fabs(d__[lend]) < fabs(d__[l])) {
129 for (m = l; m <= i__1; ++m) {
130 if (fabs(e[m]) <= eps2 * fabs(d__[m] * d__[m + 1])) {
147 F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
158 if (jtot == nmaxit) {
164 sigma = (d__[l + 1] - p) / (rte * 2.);
165 r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
166 sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
170 gamma = d__[m] - sigma;
174 for (i__ = m - 1; i__ >= i__1; --i__) {
178 e[i__ + 1] = s * r__;
185 gamma = c__ * (alpha - sigma) - s * oldgam;
186 d__[i__ + 1] = oldgam + (alpha - gamma);
187 if (fabs(c__)>GMX_FLOAT_MIN) {
188 p = gamma * gamma / c__;
195 d__[l] = sigma + gamma;
211 for (m = l; m >= i__1; --m) {
212 if (fabs(e[m - 1]) <= eps2 * fabs(d__[m] * d__[m - 1])) {
228 rte = sqrt(e[l - 1]);
229 F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
240 if (jtot == nmaxit) {
245 rte = sqrt(e[l - 1]);
246 sigma = (d__[l - 1] - p) / (rte * 2.);
247 r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
248 sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
252 gamma = d__[m] - sigma;
256 for (i__ = m; i__ <= i__1; ++i__) {
260 e[i__ - 1] = s * r__;
266 alpha = d__[i__ + 1];
267 gamma = c__ * (alpha - sigma) - s * oldgam;
268 d__[i__] = oldgam + (alpha - gamma);
269 if (fabs(c__)>GMX_FLOAT_MIN) {
270 p = gamma * gamma / c__;
277 d__[l] = sigma + gamma;
293 i__1 = lendsv - lsv + 1;
294 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
298 i__1 = lendsv - lsv + 1;
299 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
307 for (i__ = 1; i__ <= i__1; ++i__) {
308 if (fabs(e[i__])>GMX_FLOAT_MIN) {