2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
38 #include <types/simple.h>
41 #include "gmx_lapack.h"
42 #include "lapack_limits.h"
47 F77_FUNC(dlarrex,DLARREX)(const char *range,
77 double eps, tau, nrm, tmp, vvl, vvu, offd;
78 int iend, jblk, till, itmp;
79 double rtol, delta, sigma;
103 if (*range=='A' || *range=='a')
105 else if (*range=='V' || *range=='v')
107 else if (*range=='I' || *range=='i')
112 eps = GMX_DOUBLE_EPS;
116 for (i__ = 1; i__ <= i__1; ++i__) {
117 if (fabs(e[i__]) <= *tol) {
118 isplit[*nsplit] = i__;
122 isplit[*nsplit] = *n;
126 for (jblk = 1; jblk <= i__1; ++jblk) {
128 if (ibegin == iend) {
137 in = iend - ibegin + 1;
139 gl = d__[ibegin] - fabs(e[ibegin]);
140 gu = d__[ibegin] + fabs(e[ibegin]);
141 gersch[(ibegin << 1) - 1] = gl;
142 gersch[ibegin * 2] = gu;
143 gersch[(iend << 1) - 1] = d__[iend] - fabs(e[iend - 1]);
144 gersch[iend * 2] = d__[iend] + fabs(e[iend - 1]);
145 d__1 = gersch[(iend << 1) - 1];
146 gl = (d__1<gl) ? d__1 : gl;
147 d__1 = gersch[iend * 2];
148 gu = (d__1>gu) ? d__1 : gu;
150 for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
151 offd = fabs(e[i__ - 1]) + fabs(e[i__]);
152 gersch[(i__ << 1) - 1] = d__[i__] - offd;
153 d__1 = gersch[(i__ << 1) - 1];
154 gl = (d__1<gl) ? d__1 : gl;
155 gersch[i__ * 2] = d__[i__] + offd;
156 d__1 = gersch[i__ * 2];
157 gu = (d__1>gu) ? d__1 : gu;
159 d__1 = fabs(gl), d__2 = fabs(gu);
160 nrm = (d__1>d__2) ? d__1 : d__2;
164 for (i__ = ibegin; i__ <= i__2; ++i__) {
165 work[i__] = e[i__] * e[i__];
167 for (j = 1; j <= 2; ++j) {
169 tau = gl + width * .25;
171 tau = gu - width * .25;
173 tmp = d__[ibegin] - tau;
180 for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
181 tmp = d__[i__] - tau - work[i__ - 1] / tmp;
188 } else if (cnt == in) {
196 if (in - cnt > maxcnt) {
207 sigma -= delta * tau;
208 work[1] = d__[ibegin] - sigma;
211 for (i__ = 1; i__ <= i__2; ++i__) {
212 work[(in << 1) + i__] = 1. / work[i__];
213 tmp = e[j] * work[(in << 1) + i__];
214 work[i__ + 1] = d__[j + 1] - sigma - tmp * e[j];
215 work[in + i__] = tmp;
218 for (i__ = in; i__ >= 1; --i__) {
219 tmp = sgndef * work[i__];
220 if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_DOUBLE_MIN || ! (tmp > 0. || tmp < 1.)) {
226 F77_FUNC(dcopy,DCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
228 F77_FUNC(dcopy,DCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
230 for (i__ = 1; i__ <= i__2; ++i__) {
231 work[in * 3 + i__] = work[i__] * work[in + i__];
232 work[(in << 2) + i__] = work[in * 3 + i__] * work[in + i__];
236 work[1] = (gl + gu) / 2. - sigma;
238 work[(in << 1) + 1] = (gu - gl) / 2.;
241 work[in] = (gl + gu) / 2. - sigma;
243 work[in * 3] = (gu - gl) / 2.;
246 F77_FUNC(dlarrbx,DLARRBX)(&in, &d__[ibegin], &e[ibegin], &work[in * 3 + 1], &work[(in <<
247 2) + 1], &cnt, &cnt, &rtol, &rtol, &c__0, &work[1], &work[in
248 + 1], &work[(in << 1) + 1], &work[in * 5 + 1], &iwork[1], &
251 tau = work[1] - work[(in << 1) + 1];
253 tau = work[in] + work[in * 3];
264 for (i__ = 1; i__ <= i__2; ++i__) {
265 work[i__] = d__[j] + s;
266 work[(in << 1) + i__] = 1. / work[i__];
267 work[in + i__] = e[j] * d__[j] * work[(in << 1) + i__];
268 s = s * work[in + i__] * e[j] - tau;
271 work[in] = d__[iend] + s;
273 for (i__ = in; i__ >= 1; --i__) {
274 tmp = sgndef * work[i__];
275 if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_DOUBLE_MIN || ! (tmp > 0. || tmp < 1.)) {
282 F77_FUNC(dcopy,DCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
284 F77_FUNC(dcopy,DCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
286 tmp = (double) in * 4. * eps * (fabs(sigma) + fabs(tau));
288 for (i__ = ibegin; i__ <= i__2; ++i__) {
289 gersch[(i__ << 1) - 1] = gersch[(i__ << 1) - 1] - sigma - tmp;
290 gersch[i__ * 2] = gersch[i__ * 2] - sigma + tmp;
295 for (i__ = 1; i__ <= i__2; ++i__) {
296 work[(i__ << 1) - 1] = fabs(d__[j]);
297 work[i__ * 2] = e[j] * e[j] * work[(i__ << 1) - 1];
300 work[(in << 1) - 1] = fabs(d__[iend]);
302 F77_FUNC(dlasq2,DLASQ2)(&in, &work[1], info);
309 for (i__ = 1; i__ <= i__2; ++i__) {
311 w[*m] = work[in - i__ + 1];
317 for (i__ = 1; i__ <= i__2; ++i__) {
332 for (i__ = 1; i__ <= i__1; ++i__) {
337 for (j = ibegin; j <= i__2; ++j) {
338 if (vvl <= w[j] && w[j] <= vvu) {
342 indexw[*m] = j - ibegin + 1;
347 } else if (irange == 3) {
351 for (i__ = 1; i__ <= i__1; ++i__) {
352 w[i__] = w[*il + i__ - 1];
353 indexw[i__] = *il + i__ - 1;
358 for (i__ = 1; i__ <= i__1; ++i__) {
361 for (j = ibegin; j <= i__2; ++j) {
362 work[j] = w[j] + e[iend];
367 for (i__ = 1; i__ <= i__1; ++i__) {
369 iwork[*n + i__] = iblock[i__];
371 F77_FUNC(dlasrt2,DLASRT2)("I", n, &work[1], &iwork[1], &iinfo);
373 for (i__ = 1; i__ <= i__1; ++i__) {
374 itmp = iwork[*il + i__ - 1];
376 iblock[i__] = iwork[*n + itmp];
379 for (i__ = 1; i__ <= i__1; ++i__) {
380 iwork[*n + i__] = iwork[*il + i__ - 1];
383 F77_FUNC(ilasrt2,ILASRT2)("I", m, &iblock[1], &iwork[1], &iinfo);
386 cnt = iwork[*n + iwork[j]];
390 ibegin = isplit[itmp - 1] + 1;
393 for (i__ = 1; i__ <= i__1; ++i__) {
394 w[i__] = work[iwork[i__]];
395 if (iblock[i__] != itmp || i__ == *m) {
396 if (iblock[i__] == itmp) {
402 F77_FUNC(dlasrt,DLASRT)("I", &i__2, &w[j], &iinfo);
403 cnt = cnt - ibegin + 1;
405 for (k = j; k <= i__2; ++k) {
406 indexw[k] = cnt + k - j;
410 cnt = iwork[*n + iwork[j]];
411 ibegin = isplit[itmp - 1] + 1;
412 if (i__ == *m && till < *m) {
413 indexw[*m] = cnt - ibegin + 1;
416 i__2 = cnt, i__3 = iwork[*n + iwork[i__]];
417 cnt = (i__2<i__3) ? i__2 : i__3;