4 #include "types/simple.h"
6 #include "../gmx_blas.h"
7 #include "../gmx_lapack.h"
8 #include "lapack_limits.h"
13 F77_FUNC(slarrex,SLARREX)(const char *range,
43 float eps, tau, nrm, tmp, vvl, vvu, offd;
44 int iend, jblk, till, itmp;
45 float rtol, delta, sigma;
69 if (*range=='A' || *range=='a')
71 else if (*range=='V' || *range=='v')
73 else if (*range=='I' || *range=='i')
82 for (i__ = 1; i__ <= i__1; ++i__) {
83 if (fabs(e[i__]) <= *tol) {
84 isplit[*nsplit] = i__;
92 for (jblk = 1; jblk <= i__1; ++jblk) {
103 in = iend - ibegin + 1;
105 gl = d__[ibegin] - fabs(e[ibegin]);
106 gu = d__[ibegin] + fabs(e[ibegin]);
107 gersch[(ibegin << 1) - 1] = gl;
108 gersch[ibegin * 2] = gu;
109 gersch[(iend << 1) - 1] = d__[iend] - fabs(e[iend - 1]);
110 gersch[iend * 2] = d__[iend] + fabs(e[iend - 1]);
111 d__1 = gersch[(iend << 1) - 1];
112 gl = (d__1<gl) ? d__1 : gl;
113 d__1 = gersch[iend * 2];
114 gu = (d__1>gu) ? d__1 : gu;
116 for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
117 offd = fabs(e[i__ - 1]) + fabs(e[i__]);
118 gersch[(i__ << 1) - 1] = d__[i__] - offd;
119 d__1 = gersch[(i__ << 1) - 1];
120 gl = (d__1<gl) ? d__1 : gl;
121 gersch[i__ * 2] = d__[i__] + offd;
122 d__1 = gersch[i__ * 2];
123 gu = (d__1>gu) ? d__1 : gu;
125 d__1 = fabs(gl), d__2 = fabs(gu);
126 nrm = (d__1>d__2) ? d__1 : d__2;
130 for (i__ = ibegin; i__ <= i__2; ++i__) {
131 work[i__] = e[i__] * e[i__];
133 for (j = 1; j <= 2; ++j) {
135 tau = gl + width * .25;
137 tau = gu - width * .25;
139 tmp = d__[ibegin] - tau;
146 for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
147 tmp = d__[i__] - tau - work[i__ - 1] / tmp;
154 } else if (cnt == in) {
162 if (in - cnt > maxcnt) {
173 sigma -= delta * tau;
174 work[1] = d__[ibegin] - sigma;
177 for (i__ = 1; i__ <= i__2; ++i__) {
178 work[(in << 1) + i__] = 1. / work[i__];
179 tmp = e[j] * work[(in << 1) + i__];
180 work[i__ + 1] = d__[j + 1] - sigma - tmp * e[j];
181 work[in + i__] = tmp;
184 for (i__ = in; i__ >= 1; --i__) {
185 tmp = sgndef * work[i__];
186 if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_FLOAT_MIN || ! (tmp > 0. || tmp < 1.)) {
192 F77_FUNC(scopy,SCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
194 F77_FUNC(scopy,SCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
196 for (i__ = 1; i__ <= i__2; ++i__) {
197 work[in * 3 + i__] = work[i__] * work[in + i__];
198 work[(in << 2) + i__] = work[in * 3 + i__] * work[in + i__];
202 work[1] = (gl + gu) / 2. - sigma;
204 work[(in << 1) + 1] = (gu - gl) / 2.;
207 work[in] = (gl + gu) / 2. - sigma;
209 work[in * 3] = (gu - gl) / 2.;
212 F77_FUNC(slarrbx,SLARRBX)(&in, &d__[ibegin], &e[ibegin], &work[in * 3 + 1], &work[(in <<
213 2) + 1], &cnt, &cnt, &rtol, &rtol, &c__0, &work[1], &work[in
214 + 1], &work[(in << 1) + 1], &work[in * 5 + 1], &iwork[1], &
217 tau = work[1] - work[(in << 1) + 1];
219 tau = work[in] + work[in * 3];
230 for (i__ = 1; i__ <= i__2; ++i__) {
231 work[i__] = d__[j] + s;
232 work[(in << 1) + i__] = 1. / work[i__];
233 work[in + i__] = e[j] * d__[j] * work[(in << 1) + i__];
234 s = s * work[in + i__] * e[j] - tau;
237 work[in] = d__[iend] + s;
239 for (i__ = in; i__ >= 1; --i__) {
240 tmp = sgndef * work[i__];
241 if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_FLOAT_MIN || ! (tmp > 0. || tmp < 1.)) {
248 F77_FUNC(scopy,SCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
250 F77_FUNC(scopy,SCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
252 tmp = (float) in * 4. * eps * (fabs(sigma) + fabs(tau));
254 for (i__ = ibegin; i__ <= i__2; ++i__) {
255 gersch[(i__ << 1) - 1] = gersch[(i__ << 1) - 1] - sigma - tmp;
256 gersch[i__ * 2] = gersch[i__ * 2] - sigma + tmp;
261 for (i__ = 1; i__ <= i__2; ++i__) {
262 work[(i__ << 1) - 1] = fabs(d__[j]);
263 work[i__ * 2] = e[j] * e[j] * work[(i__ << 1) - 1];
266 work[(in << 1) - 1] = fabs(d__[iend]);
268 F77_FUNC(slasq2,SLASQ2)(&in, &work[1], info);
275 for (i__ = 1; i__ <= i__2; ++i__) {
277 w[*m] = work[in - i__ + 1];
283 for (i__ = 1; i__ <= i__2; ++i__) {
298 for (i__ = 1; i__ <= i__1; ++i__) {
303 for (j = ibegin; j <= i__2; ++j) {
304 if (vvl <= w[j] && w[j] <= vvu) {
308 indexw[*m] = j - ibegin + 1;
313 } else if (irange == 3) {
317 for (i__ = 1; i__ <= i__1; ++i__) {
318 w[i__] = w[*il + i__ - 1];
319 indexw[i__] = *il + i__ - 1;
324 for (i__ = 1; i__ <= i__1; ++i__) {
327 for (j = ibegin; j <= i__2; ++j) {
328 work[j] = w[j] + e[iend];
333 for (i__ = 1; i__ <= i__1; ++i__) {
335 iwork[*n + i__] = iblock[i__];
337 F77_FUNC(slasrt2,SLASRT2)("I", n, &work[1], &iwork[1], &iinfo);
339 for (i__ = 1; i__ <= i__1; ++i__) {
340 itmp = iwork[*il + i__ - 1];
342 iblock[i__] = iwork[*n + itmp];
345 for (i__ = 1; i__ <= i__1; ++i__) {
346 iwork[*n + i__] = iwork[*il + i__ - 1];
349 F77_FUNC(ilasrt2,ILASRT2)("I", m, &iblock[1], &iwork[1], &iinfo);
352 cnt = iwork[*n + iwork[j]];
356 ibegin = isplit[itmp - 1] + 1;
359 for (i__ = 1; i__ <= i__1; ++i__) {
360 w[i__] = work[iwork[i__]];
361 if (iblock[i__] != itmp || i__ == *m) {
362 if (iblock[i__] == itmp) {
368 F77_FUNC(slasrt,SLASRT)("I", &i__2, &w[j], &iinfo);
369 cnt = cnt - ibegin + 1;
371 for (k = j; k <= i__2; ++k) {
372 indexw[k] = cnt + k - j;
376 cnt = iwork[*n + iwork[j]];
377 ibegin = isplit[itmp - 1] + 1;
378 if (i__ == *m && till < *m) {
379 indexw[*m] = cnt - ibegin + 1;
382 i__2 = cnt, i__3 = iwork[*n + iwork[i__]];
383 cnt = (i__2<i__3) ? i__2 : i__3;