2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
5 #include "gromacs/utility/real.h"
8 #pragma warning(disable: 4723) /*division by zero - is used on purpose here*/
12 F77_FUNC(slasq2,SLASQ2)(int *n,
28 float dmin__, emin, emax;
30 float qmin, temp, qmax, zmax;
32 float desig, trace, sigma;
36 float oldemn, safmin, minval;
37 float posinf,neginf,negzro,newzro;
45 minval = GMX_FLOAT_MIN;
46 safmin = minval*(1.0+eps);
66 if (z__[2] < 0. || z__[3] < 0.) {
69 } else if (z__[3] > z__[1]) {
74 z__[5] = z__[1] + z__[2] + z__[3];
75 if (z__[2] > z__[3] * tol2) {
76 t = (z__[1] - z__[3] + z__[2]) * .5;
77 s = z__[3] * (z__[2] / t);
79 s = z__[3] * (z__[2] / (t * (sqrt(s / t + 1.) + 1.)));
81 s = z__[3] * (z__[2] / (t + sqrt(t) * sqrt(t + s)));
83 t = z__[1] + (s + z__[2]);
88 z__[6] = z__[2] + z__[1];
99 for (k = 1; k <= i__1; k += 2) {
103 } else if (z__[k + 1] < 0.) {
109 d__1 = qmax, d__2 = z__[k];
110 qmax = (d__1>d__2) ? d__1 : d__2;
111 d__1 = emin, d__2 = z__[k + 1];
112 emin = (d__1<d__2) ? d__1 : d__2;
113 d__1 = (qmax>zmax) ? qmax : zmax;
115 zmax = (d__1>d__2) ? d__1 : d__2;
117 if (z__[(*n << 1) - 1] < 0.) {
118 *info = -((*n << 1) + 199);
121 d__ += z__[(*n << 1) - 1];
122 d__1 = qmax, d__2 = z__[(*n << 1) - 1];
123 qmax = (d__1>d__2) ? d__1 : d__2;
124 zmax = (qmax>zmax) ? qmax : zmax;
126 if (fabs(e)<GMX_FLOAT_MIN) {
128 for (k = 2; k <= i__1; ++k) {
129 z__[k] = z__[(k << 1) - 1];
131 F77_FUNC(slasrt,SLASRT)("D", n, &z__[1], &iinfo);
132 z__[(*n << 1) - 1] = d__;
138 if (fabs(trace)<GMX_FLOAT_MIN) {
139 z__[(*n << 1) - 1] = 0.;
150 negzro = one/(neginf+one);
151 if(fabs(negzro)>GMX_FLOAT_MIN)
156 newzro = negzro + zero;
157 if(fabs(newzro-zero)>GMX_FLOAT_MIN)
159 posinf = one /newzro;
162 neginf = neginf*posinf;
165 posinf = posinf*posinf;
169 for (k = *n << 1; k >= 2; k += -2) {
171 z__[(k << 1) - 1] = z__[k];
172 z__[(k << 1) - 2] = 0.;
173 z__[(k << 1) - 3] = z__[k - 1];
179 if (z__[(i0 << 2) - 3] * 1.5 < z__[(n0 << 2) - 3]) {
181 i__1 = 2*(i0 + n0 - 1);
182 for (i4 = i0 << 2; i4 <= i__1; i4 += 4) {
184 z__[i4 - 3] = z__[ipn4 - i4 - 3];
185 z__[ipn4 - i4 - 3] = temp;
187 z__[i4 - 1] = z__[ipn4 - i4 - 5];
188 z__[ipn4 - i4 - 5] = temp;
194 for (k = 1; k <= 2; ++k) {
196 d__ = z__[(n0 << 2) + pp - 3];
197 i__1 = (i0 << 2) + pp;
198 for (i4 = 4*(n0 - 1) + pp; i4 >= i__1; i4 += -4) {
199 if (z__[i4 - 1] <= tol2 * d__) {
203 d__ = z__[i4 - 3] * (d__ / (d__ + z__[i4 - 1]));
207 emin = z__[(i0 << 2) + pp + 1];
208 d__ = z__[(i0 << 2) + pp - 3];
209 i__1 = 4*(n0 - 1) + pp;
210 for (i4 = (i0 << 2) + pp; i4 <= i__1; i4 += 4) {
211 z__[i4 - (pp << 1) - 2] = d__ + z__[i4 - 1];
212 if (z__[i4 - 1] <= tol2 * d__) {
214 z__[i4 - (pp << 1) - 2] = d__;
215 z__[i4 - (pp << 1)] = 0.;
217 } else if (safmin * z__[i4 + 1] < z__[i4 - (pp << 1) - 2] &&
218 safmin * z__[i4 - (pp << 1) - 2] < z__[i4 + 1]) {
219 temp = z__[i4 + 1] / z__[i4 - (pp << 1) - 2];
220 z__[i4 - (pp << 1)] = z__[i4 - 1] * temp;
223 z__[i4 - (pp << 1)] = z__[i4 + 1] * (z__[i4 - 1] / z__[i4 - (
225 d__ = z__[i4 + 1] * (d__ / z__[i4 - (pp << 1) - 2]);
227 d__1 = emin, d__2 = z__[i4 - (pp << 1)];
228 emin = (d__1<d__2) ? d__1 : d__2;
230 z__[(n0 << 2) - pp - 2] = d__;
233 qmax = z__[(i0 << 2) - pp - 2];
234 i__1 = (n0 << 2) - pp - 2;
235 for (i4 = (i0 << 2) - pp + 2; i4 <= i__1; i4 += 4) {
236 d__1 = qmax, d__2 = z__[i4];
237 qmax = (d__1>d__2) ? d__1 : d__2;
248 for (iwhila = 1; iwhila <= i__1; ++iwhila) {
257 sigma = -z__[(n0 << 2) - 1];
266 emin = fabs(z__[(n0 << 2) - 5]);
270 qmin = z__[(n0 << 2) - 3];
272 for (i4 = n0 << 2; i4 >= 8; i4 += -4) {
273 if (z__[i4 - 5] <= 0.) {
276 if (qmin >= emax * 4.) {
277 d__1 = qmin, d__2 = z__[i4 - 3];
278 qmin = (d__1<d__2) ? d__1 : d__2;
279 d__1 = emax, d__2 = z__[i4 - 5];
280 emax = (d__1>d__2) ? d__1 : d__2;
282 d__1 = qmax, d__2 = z__[i4 - 7] + z__[i4 - 5];
283 qmax = (d__1>d__2) ? d__1 : d__2;
284 d__1 = emin, d__2 = z__[i4 - 5];
285 emin = (d__1<d__2) ? d__1 : d__2;
294 dee = z__[(i0 << 2) - 3];
297 i__2 = (n0 << 2) - 3;
298 for (i4 = (i0 << 2) - 3; i4 <= i__2; i4 += 4) {
299 dee = z__[i4] * (dee / (dee + z__[i4 - 2]));
305 if (2*(kmin - i0) < n0 - kmin && deemin <= z__[(n0 << 2) - 3] *
309 i__2 = 2*(i0 + n0 - 1);
310 for (i4 = i0 << 2; i4 <= i__2; i4 += 4) {
312 z__[i4 - 3] = z__[ipn4 - i4 - 3];
313 z__[ipn4 - i4 - 3] = temp;
315 z__[i4 - 2] = z__[ipn4 - i4 - 2];
316 z__[ipn4 - i4 - 2] = temp;
318 z__[i4 - 1] = z__[ipn4 - i4 - 5];
319 z__[ipn4 - i4 - 5] = temp;
321 z__[i4] = z__[ipn4 - i4 - 4];
322 z__[ipn4 - i4 - 4] = temp;
328 d__1 = 0., d__2 = qmin - sqrt(qmin) * 2. * sqrt(emax);
329 dmin__ = -((d__1>d__2) ? d__1 : d__2);
331 nbig = (n0 - i0 + 1) * 30;
333 for (iwhilb = 1; iwhilb <= i__2; ++iwhilb) {
338 F77_FUNC(slasq3,SLASQ3)(&i0, &n0, &z__[1], &pp, &dmin__, &sigma, &desig, &qmax, &
339 nfail, &iter, &ndiv, &ieee);
343 if (pp == 0 && n0 - i0 >= 3) {
344 if (z__[n0 * 4] <= tol2 * qmax || z__[(n0 << 2) - 1] <= tol2 *
347 qmax = z__[(i0 << 2) - 3];
348 emin = z__[(i0 << 2) - 1];
349 oldemn = z__[i0 * 4];
351 for (i4 = i0 << 2; i4 <= i__3; i4 += 4) {
352 if (z__[i4] <= tol2 * z__[i4 - 3] || z__[i4 - 1] <=
354 z__[i4 - 1] = -sigma;
358 oldemn = z__[i4 + 4];
360 d__1 = qmax, d__2 = z__[i4 + 1];
361 qmax = (d__1>d__2) ? d__1 : d__2;
362 d__1 = emin, d__2 = z__[i4 - 1];
363 emin = (d__1<d__2) ? d__1 : d__2;
364 d__1 = oldemn, d__2 = z__[i4];
365 oldemn = (d__1<d__2) ? d__1 : d__2;
368 z__[(n0 << 2) - 1] = emin;
369 z__[n0 * 4] = oldemn;
389 for (k = 2; k <= i__1; ++k) {
390 z__[k] = z__[(k << 2) - 3];
393 F77_FUNC(slasrt,SLASRT)("D", n, &z__[1], &iinfo);
396 for (k = *n; k >= 1; --k) {
401 z__[(*n << 1) + 1] = trace;
402 z__[(*n << 1) + 2] = e;
403 z__[(*n << 1) + 3] = (float) iter;
405 z__[(*n << 1) + 4] = (float) ndiv / (float) (i__1 * i__1);
406 z__[(*n << 1) + 5] = nfail * 100. / (float) iter;