2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2004 David van der Spoel, Erik Lindahl, University of Groningen.
5 * Copyright (c) 2012,2013, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gromacs/legacyheaders/types/simple.h"
40 #include "gmx_arpack.h"
42 #include "gmx_lapack.h"
44 F77_FUNC(dstqrb, DSTQRB) (int * n,
60 int l1, ii, mm, lm1, mm1, nm1;
64 int lend, jtot, lendm1, lendp1, iscale;
66 int lendsv, nmaxit, icompz;
67 double ssfmax, ssfmin, safmin, minval, safmax, anorm;
97 minval = GMX_DOUBLE_MIN;
98 safmin = minval / GMX_DOUBLE_EPS;
100 ssfmax = sqrt(safmax) / 3.;
101 ssfmin = sqrt(safmin) / eps2;
106 for (j = 1; j <= i__1; ++j)
132 for (m = l1; m <= i__1; ++m)
139 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps)
160 anorm = F77_FUNC(dlanst, DLANST) ("i", &i__1, &d__[l], &e[l]);
170 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
173 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
176 else if (anorm < ssfmin)
180 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
183 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
187 if (fabs(d__[lend]) < fabs(d__[l]))
201 for (m = l; m <= i__1; ++m)
205 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin)
229 F77_FUNC(dlaev2, DLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
231 work[*n - 1 + l] = s;
234 z__[l + 1] = c__ * tst - s * z__[l];
235 z__[l] = s * tst + c__ * z__[l];
239 F77_FUNC(dlae2, DLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
258 g = (d__[l + 1] - p) / (e[l] * 2.);
259 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
260 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
268 for (i__ = mm1; i__ >= i__1; --i__)
272 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
277 g = d__[i__ + 1] - p;
278 r__ = (d__[i__] - g) * s + c__ * 2. * b;
280 d__[i__ + 1] = g + p;
286 work[*n - 1 + i__] = -s;
295 F77_FUNC(dlasr, DLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
322 for (m = l; m >= i__1; --m)
324 d__2 = fabs(e[m - 1]);
326 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin)
350 F77_FUNC(dlaev2, DLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
354 z__[l] = c__ * tst - s * z__[l - 1];
355 z__[l - 1] = s * tst + c__ * z__[l - 1];
359 F77_FUNC(dlae2, DLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
379 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
380 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
381 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
389 for (i__ = m; i__ <= i__1; ++i__)
393 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
399 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
407 work[*n - 1 + i__] = s;
416 F77_FUNC(dlasr, DLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
439 i__1 = lendsv - lsv + 1;
440 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
443 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
446 else if (iscale == 2)
448 i__1 = lendsv - lsv + 1;
449 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
452 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
461 for (i__ = 1; i__ <= i__1; ++i__)
474 F77_FUNC(dlasrt, DLASRT) ("i", n, &d__[1], info);
481 for (ii = 2; ii <= i__1; ++ii)
487 for (j = ii; j <= i__2; ++j)
513 F77_FUNC(dgetv0, DGETV0) (int * ido,
515 int gmx_unused * itry,
532 int v_dim1, v_offset, i__1;
540 v_offset = 1 + v_dim1;
556 F77_FUNC(dlarnv, DLARNV) (&idist, &iwork[1], n, &resid[1]);
563 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
582 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
588 else if (*bmat == 'I')
590 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
599 workd[*n * 3 + 4] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
600 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
602 else if (*bmat == 'I')
604 workd[*n * 3 + 4] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
606 *rnorm = workd[*n * 3 + 4];
616 F77_FUNC(dgemv, DGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
617 &workd[*n + 1], &c__1);
619 F77_FUNC(dgemv, DGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
620 c_b22, &resid[1], &c__1);
624 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
630 else if (*bmat == 'I')
632 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
639 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
640 *rnorm = sqrt(fabs(*rnorm));
642 else if (*bmat == 'I')
644 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
647 if (*rnorm > workd[*n * 3 + 4] * .717f)
656 workd[*n * 3 + 4] = *rnorm;
663 for (jj = 1; jj <= i__1; ++jj)
684 F77_FUNC(dsapps, DSAPPS) (int * n,
701 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
705 double r__, s, a1, a2, a3, a4;
716 v_offset = 1 + v_dim1;
719 h_offset = 1 + h_dim1;
722 q_offset = 1 + q_dim1;
725 epsmch = GMX_DOUBLE_EPS;
731 F77_FUNC(dlaset, DLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
739 for (jj = 1; jj <= i__1; ++jj)
747 for (i__ = istart; i__ <= i__2; ++i__)
749 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
750 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
752 h__[i__ + 1 + h_dim1] = 0.;
763 f = h__[istart + (h_dim1 << 1)] - shift[jj];
764 g = h__[istart + 1 + h_dim1];
765 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
767 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
769 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
771 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
773 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
775 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
776 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
777 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
780 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
781 for (j = 1; j <= i__2; ++j)
783 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
785 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
786 c__ * q[j + (istart + 1) * q_dim1];
787 q[j + istart * q_dim1] = a1;
792 for (i__ = istart + 1; i__ <= i__2; ++i__)
795 f = h__[i__ + h_dim1];
796 g = s * h__[i__ + 1 + h_dim1];
798 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
799 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
808 h__[i__ + h_dim1] = r__;
810 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
812 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
814 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
816 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
819 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
820 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
821 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
824 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
825 for (j = 1; j <= i__3; ++j)
827 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
829 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
830 c__ * q[j + (i__ + 1) * q_dim1];
831 q[j + i__ * q_dim1] = a1;
840 if (h__[iend + h_dim1] < 0.)
842 h__[iend + h_dim1] = -h__[iend + h_dim1];
843 F77_FUNC(dscal, DSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
852 for (i__ = itop; i__ <= i__2; ++i__)
854 if (h__[i__ + 1 + h_dim1] > 0.)
866 for (i__ = itop; i__ <= i__1; ++i__)
868 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
869 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
871 h__[i__ + 1 + h_dim1] = 0.;
876 if (h__[*kev + 1 + h_dim1] > 0.)
878 F77_FUNC(dgemv, DGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
879 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
883 for (i__ = 1; i__ <= i__1; ++i__)
885 i__2 = kplusp - i__ + 1;
886 F77_FUNC(dgemv, DGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
887 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
888 F77_FUNC(dcopy, DCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
893 F77_FUNC(dlacpy, DLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
895 if (h__[*kev + 1 + h_dim1] > 0.)
897 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
900 F77_FUNC(dscal, DSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
901 if (h__[*kev + 1 + h_dim1] > 0.)
903 F77_FUNC(daxpy, DAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
917 F77_FUNC(dsortr, DSORTR) (const char * which,
932 if (!strncmp(which, "SA", 2))
941 for (i__ = igap; i__ <= i__1; ++i__)
951 if (x1[j] < x1[j + igap])
954 x1[j] = x1[j + igap];
959 x2[j] = x2[j + igap];
976 else if (!strncmp(which, "SM", 2))
985 for (i__ = igap; i__ <= i__1; ++i__)
995 if (fabs(x1[j]) < fabs(x1[j + igap]))
998 x1[j] = x1[j + igap];
1003 x2[j] = x2[j + igap];
1004 x2[j + igap] = temp;
1020 else if (!strncmp(which, "LA", 2))
1029 for (i__ = igap; i__ <= i__1; ++i__)
1039 if (x1[j] > x1[j + igap])
1042 x1[j] = x1[j + igap];
1043 x1[j + igap] = temp;
1047 x2[j] = x2[j + igap];
1048 x2[j + igap] = temp;
1064 else if (!strncmp(which, "LM", 2))
1074 for (i__ = igap; i__ <= i__1; ++i__)
1084 if (fabs(x1[j]) > fabs(x1[j + igap]))
1087 x1[j] = x1[j + igap];
1088 x1[j + igap] = temp;
1092 x2[j] = x2[j + igap];
1093 x2[j + igap] = temp;
1118 F77_FUNC(dsesrt, DSESRT) (const char * which,
1126 int a_dim1, a_offset, i__1;
1133 a_offset = 1 + a_dim1 * 0;
1138 if (!strncmp(which, "SA", 2))
1147 for (i__ = igap; i__ <= i__1; ++i__)
1157 if (x[j] < x[j + igap])
1164 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1165 a_dim1 + 1], &c__1);
1181 else if (!strncmp(which, "SM", 2))
1190 for (i__ = igap; i__ <= i__1; ++i__)
1200 if (fabs(x[j]) < fabs(x[j + igap]))
1207 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1208 a_dim1 + 1], &c__1);
1224 else if (!strncmp(which, "LA", 2))
1233 for (i__ = igap; i__ <= i__1; ++i__)
1243 if (x[j] > x[j + igap])
1250 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1251 a_dim1 + 1], &c__1);
1267 else if (!strncmp(which, "LM", 2))
1276 for (i__ = igap; i__ <= i__1; ++i__)
1286 if (fabs(x[j]) > fabs(x[j + igap]))
1293 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1294 a_dim1 + 1], &c__1);
1319 F77_FUNC(dsgets, DSGETS) (int * ishift,
1335 if (!strncmp(which, "BE", 2))
1338 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1342 i__1 = (kevd2 < *np) ? kevd2 : *np;
1343 i__2 = (kevd2 > *np) ? kevd2 : *np;
1344 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[1], &c__1,
1345 &ritz[i__2 + 1], &c__1);
1346 i__1 = (kevd2 < *np) ? kevd2 : *np;
1347 i__2 = (kevd2 > *np) ? kevd2 : *np;
1348 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[1], &c__1,
1349 &bounds[i__2 + 1], &c__1);
1356 F77_FUNC(dsortr, DSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
1359 if (*ishift == 1 && *np > 0)
1362 F77_FUNC(dsortr, DSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
1363 F77_FUNC(dcopy, DCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
1373 F77_FUNC(dsconv, DSCONV) (int * n,
1389 eps23 = GMX_DOUBLE_EPS;
1390 eps23 = pow(eps23, c_b3);
1394 for (i__ = 1; i__ <= i__1; ++i__)
1398 d__3 = fabs(ritz[i__]);
1399 temp = (d__2 > d__3) ? d__2 : d__3;
1400 if (bounds[i__] <= *tol * temp)
1411 F77_FUNC(dseigt, DSEIGT) (double * rnorm,
1421 int h_dim1, h_offset, i__1;
1430 h_offset = 1 + h_dim1;
1433 F77_FUNC(dcopy, DCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1435 F77_FUNC(dcopy, DCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1436 F77_FUNC(dstqrb, DSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1443 for (k = 1; k <= i__1; ++k)
1445 bounds[k] = *rnorm * fabs(bounds[k]);
1458 F77_FUNC(dsaitr, DSAITR) (int * ido,
1482 int h_dim1, h_offset, v_dim1, v_offset, i__1;
1486 double safmin, minval;
1492 v_offset = 1 + v_dim1;
1495 h_offset = 1 + h_dim1;
1499 minval = GMX_DOUBLE_MIN;
1500 safmin = minval / GMX_DOUBLE_EPS;
1514 iwork[9] = iwork[8] + *n;
1515 iwork[10] = iwork[9] + *n;
1553 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
1554 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
1567 *info = iwork[12] - 1;
1574 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1575 if (*rnorm >= safmin)
1577 temp1 = 1. / *rnorm;
1578 F77_FUNC(dscal, DSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1579 F77_FUNC(dscal, DSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
1584 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
1585 v_dim1 + 1], n, &infol);
1586 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
1591 F77_FUNC(dcopy, DCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1592 ipntr[1] = iwork[10];
1593 ipntr[2] = iwork[9];
1594 ipntr[3] = iwork[8];
1603 F77_FUNC(dcopy, DCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1612 ipntr[1] = iwork[9];
1613 ipntr[2] = iwork[8];
1618 else if (*bmat == 'I')
1620 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1630 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
1632 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1634 else if (*bmat == 'G')
1636 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1638 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1640 else if (*bmat == 'I')
1642 workd[*n * 3 + 3] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1647 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
1648 &c__1, &c_b42, &workd[iwork[9]], &c__1);
1650 else if (*mode == 2)
1652 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
1653 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
1656 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1657 c__1, &c_b18, &resid[1], &c__1);
1659 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1660 if (iwork[12] == 1 || iwork[4] == 1)
1662 h__[iwork[12] + h_dim1] = 0.;
1666 h__[iwork[12] + h_dim1] = *rnorm;
1674 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1675 ipntr[1] = iwork[9];
1676 ipntr[2] = iwork[8];
1681 else if (*bmat == 'I')
1683 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1691 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1692 *rnorm = sqrt(fabs(*rnorm));
1694 else if (*bmat == 'I')
1696 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1699 if (*rnorm > workd[*n * 3 + 3] * .717f)
1706 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
1707 c__1, &c_b42, &workd[iwork[9]], &c__1);
1709 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1710 c__1, &c_b18, &resid[1], &c__1);
1712 if (iwork[12] == 1 || iwork[4] == 1)
1714 h__[iwork[12] + h_dim1] = 0.;
1716 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1721 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1722 ipntr[1] = iwork[9];
1723 ipntr[2] = iwork[8];
1728 else if (*bmat == 'I')
1730 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1737 workd[*n * 3 + 2] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1739 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
1741 else if (*bmat == 'I')
1743 workd[*n * 3 + 2] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1747 if (workd[*n * 3 + 2] > *rnorm * .717f)
1750 *rnorm = workd[*n * 3 + 2];
1756 *rnorm = workd[*n * 3 + 2];
1764 for (jj = 1; jj <= i__1; ++jj)
1776 if (h__[iwork[12] + h_dim1] < 0.)
1778 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1779 if (iwork[12] < *k + *np)
1781 F77_FUNC(dscal, DSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1785 F77_FUNC(dscal, DSCAL) (n, &c_b50, &resid[1], &c__1);
1790 if (iwork[12] > *k + *np)
1810 F77_FUNC(dsaup2, DSAUP2) (int * ido,
1819 int gmx_unused * iupd,
1840 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
1859 v_offset = 1 + v_dim1;
1862 h_offset = 1 + h_dim1;
1865 q_offset = 1 + q_dim1;
1869 eps23 = GMX_DOUBLE_EPS;
1870 eps23 = pow(eps23, c_b3);
1883 iwork[7] = iwork[9] + iwork[10];
1906 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
1907 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
1915 if (workd[*n * 3 + 1] == 0.)
1940 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
1941 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1967 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1],
1968 &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1986 F77_FUNC(dseigt, DSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
1987 bounds[1], &workl[1], &ierr);
1995 F77_FUNC(dcopy, DCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1996 F77_FUNC(dcopy, DCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
2000 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2002 F77_FUNC(dcopy, DCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
2003 F77_FUNC(dsconv, DSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
2007 for (j = 1; j <= i__1; ++j)
2009 if (bounds[j] == 0.)
2016 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
2019 if (!strncmp(which, "BE", 2))
2022 strncpy(wprime, "SA", 2);
2023 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2025 nevm2 = *nev - nevd2;
2028 i__1 = (nevd2 < *np) ? nevd2 : *np;
2029 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
2030 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
2031 &ritz[((i__2 > i__3) ? i__2 : i__3)],
2033 i__1 = (nevd2 < *np) ? nevd2 : *np;
2034 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
2035 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
2036 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
2044 if (!strncmp(which, "LM", 2))
2046 strncpy(wprime, "SM", 2);
2048 if (!strncmp(which, "SM", 2))
2050 strncpy(wprime, "LM", 2);
2052 if (!strncmp(which, "LA", 2))
2054 strncpy(wprime, "SA", 2);
2056 if (!strncmp(which, "SA", 2))
2058 strncpy(wprime, "LA", 2);
2061 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2066 for (j = 1; j <= i__1; ++j)
2069 d__3 = fabs(ritz[j]);
2070 temp = (d__2 > d__3) ? d__2 : d__3;
2074 strncpy(wprime, "LA", 2);
2075 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
2078 for (j = 1; j <= i__1; ++j)
2081 d__3 = fabs(ritz[j]);
2082 temp = (d__2 > d__3) ? d__2 : d__3;
2086 if (!strncmp(which, "BE", 2))
2089 strncpy(wprime, "LA", 2);
2090 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2095 F77_FUNC(dsortr, DSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2099 h__[h_dim1 + 1] = workd[*n * 3 + 1];
2102 if (iwork[6] > *mxiter && iwork[8] < *nev)
2106 if (*np == 0 && iwork[8] < iwork[9])
2115 else if (iwork[8] < *nev && *ishift == 1)
2118 i__1 = iwork[8], i__2 = *np / 2;
2119 *nev += (i__1 < i__2) ? i__1 : i__2;
2120 if (*nev == 1 && iwork[7] >= 6)
2122 *nev = iwork[7] / 2;
2124 else if (*nev == 1 && iwork[7] > 2)
2128 *np = iwork[7] - *nev;
2133 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2153 F77_FUNC(dcopy, DCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
2156 F77_FUNC(dsapps, DSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
2157 resid[1], &q[q_offset], ldq, &workd[1]);
2162 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
2169 else if (*bmat == 'I')
2171 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
2178 workd[*n * 3 + 1] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
2179 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
2181 else if (*bmat == 'I')
2183 workd[*n * 3 + 1] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
2205 F77_FUNC(dsaupd, DSAUPD) (int * ido,
2223 int v_dim1, v_offset, i__1, i__2;
2229 v_offset = 1 + v_dim1;
2241 iwork[5] = iparam[1];
2242 iwork[10] = iparam[3];
2243 iwork[12] = iparam[4];
2246 iwork[11] = iparam[7];
2257 else if (*ncv <= *nev || *ncv > *n)
2263 iwork[15] = *ncv - *nev;
2269 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
2270 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
2271 strncmp(which, "BE", 2))
2275 if (*bmat != 'I' && *bmat != 'G')
2281 if (*lworkl < i__1 * i__1 + (*ncv << 3))
2285 if (iwork[11] < 1 || iwork[11] > 5)
2289 else if (iwork[11] == 1 && *bmat == 'G')
2293 else if (iwork[5] < 0 || iwork[5] > 1)
2297 else if (*nev == 1 && !strncmp(which, "BE", 2))
2315 *tol = GMX_DOUBLE_EPS;
2318 iwork[15] = *ncv - *nev;
2321 i__1 = i__2 * i__2 + (*ncv << 3);
2322 for (j = 1; j <= i__1; ++j)
2330 iwork[16] = iwork[3] + (iwork[8] << 1);
2331 iwork[1] = iwork[16] + *ncv;
2332 iwork[4] = iwork[1] + *ncv;
2334 iwork[7] = iwork[4] + i__1 * i__1;
2335 iwork[14] = iwork[7] + *ncv * 3;
2337 ipntr[4] = iwork[14];
2338 ipntr[5] = iwork[3];
2339 ipntr[6] = iwork[16];
2340 ipntr[7] = iwork[1];
2341 ipntr[11] = iwork[7];
2344 F77_FUNC(dsaup2, DSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
2345 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
2346 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
2347 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
2352 iparam[8] = iwork[15];
2359 iparam[3] = iwork[10];
2360 iparam[5] = iwork[15];
2380 F77_FUNC(dseupd, DSEUPD) (int * rvec,
2381 const char * howmny,
2406 int v_dim1, v_offset, z_dim1, z_offset, i__1;
2407 double d__1, d__2, d__3;
2409 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2421 double thres1 = 0, thres2 = 0;
2425 int leftptr, rghtptr;
2431 z_offset = 1 + z_dim1;
2436 v_offset = 1 + v_dim1;
2464 if (*ncv <= *nev || *ncv > *n)
2468 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
2469 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
2470 strncmp(which, "BE", 2))
2474 if (*bmat != 'I' && *bmat != 'G')
2478 if (*howmny != 'A' && *howmny != 'P' &&
2479 *howmny != 'S' && *rvec)
2483 if (*rvec && *howmny == 'S')
2488 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
2493 if (mode == 1 || mode == 2)
2495 strncpy(type__, "REGULR", 6);
2499 strncpy(type__, "SHIFTI", 6);
2503 strncpy(type__, "BUCKLE", 6);
2507 strncpy(type__, "CAYLEY", 6);
2513 if (mode == 1 && *bmat == 'G')
2517 if (*nev == 1 && !strncmp(which, "BE", 2))
2536 iw = iq + ldh * *ncv;
2537 next = iw + (*ncv << 1);
2543 irz = ipntr[11] + *ncv;
2547 eps23 = GMX_DOUBLE_EPS;
2548 eps23 = pow(eps23, c_b21);
2555 else if (*bmat == 'G')
2557 bnorm2 = F77_FUNC(dnrm2, DNRM2) (n, &workd[1], &c__1);
2563 if (!strncmp(which, "LM", 2) || !strncmp(which, "SM", 2) ||
2564 !strncmp(which, "LA", 2) || !strncmp(which, "SA", 2))
2568 else if (!strncmp(which, "BE", 2))
2572 ism = (*nev > nconv) ? *nev : nconv;
2575 thres1 = workl[ism];
2576 thres2 = workl[ilg];
2584 for (j = 0; j <= i__1; ++j)
2587 if (!strncmp(which, "LM", 2))
2589 if (fabs(workl[irz + j]) >= fabs(thres1))
2592 d__3 = fabs(workl[irz + j]);
2593 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2594 if (workl[ibd + j] <= *tol * tempbnd)
2600 else if (!strncmp(which, "SM", 2))
2602 if (fabs(workl[irz + j]) <= fabs(thres1))
2605 d__3 = fabs(workl[irz + j]);
2606 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2607 if (workl[ibd + j] <= *tol * tempbnd)
2613 else if (!strncmp(which, "LA", 2))
2615 if (workl[irz + j] >= thres1)
2618 d__3 = fabs(workl[irz + j]);
2619 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2620 if (workl[ibd + j] <= *tol * tempbnd)
2626 else if (!strncmp(which, "SA", 2))
2628 if (workl[irz + j] <= thres1)
2631 d__3 = fabs(workl[irz + j]);
2632 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2633 if (workl[ibd + j] <= *tol * tempbnd)
2639 else if (!strncmp(which, "BE", 2))
2641 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
2644 d__3 = fabs(workl[irz + j]);
2645 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2646 if (workl[ibd + j] <= *tol * tempbnd)
2654 reord = select[j + 1] || reord;
2663 F77_FUNC(dcopy, DCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2664 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2666 F77_FUNC(dsteqr, DSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
2688 if (select[leftptr])
2694 else if (!select[rghtptr])
2703 temp = workl[ihd + leftptr - 1];
2704 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2705 workl[ihd + rghtptr - 1] = temp;
2706 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
2708 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
2709 iq + *ncv * (leftptr - 1)], &c__1);
2710 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
2717 if (leftptr < rghtptr)
2726 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2732 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2733 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2736 if (!strncmp(type__, "REGULR", 6))
2741 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2745 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2752 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2753 if (!strncmp(type__, "SHIFTI", 6))
2756 for (k = 1; k <= i__1; ++k)
2758 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2761 else if (!strncmp(type__, "BUCKLE", 6))
2764 for (k = 1; k <= i__1; ++k)
2766 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
2770 else if (!strncmp(type__, "CAYLEY", 6))
2773 for (k = 1; k <= i__1; ++k)
2775 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
2776 workl[ihd + k - 1] - 1.);
2780 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2781 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2784 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2788 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2789 d__1 = bnorm2 / rnorm;
2790 F77_FUNC(dscal, DSCAL) (ncv, &d__1, &workl[ihb], &c__1);
2791 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2796 if (*rvec && *howmny == 'A')
2799 F77_FUNC(dgeqr2, DGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2802 F77_FUNC(dorm2r, DORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
2803 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
2804 F77_FUNC(dlacpy, DLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2807 for (j = 1; j <= i__1; ++j)
2809 workl[ihb + j - 1] = 0.;
2811 workl[ihb + *ncv - 1] = 1.;
2812 F77_FUNC(dorm2r, DORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
2813 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
2816 else if (*rvec && *howmny == 'S')
2821 if (!strncmp(type__, "REGULR", 6) && *rvec)
2825 for (j = 1; j <= i__1; ++j)
2827 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
2831 else if (strncmp(type__, "REGULR", 6) && *rvec)
2834 F77_FUNC(dscal, DSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
2835 if (!strncmp(type__, "SHIFTI", 6))
2839 for (k = 1; k <= i__1; ++k)
2841 d__2 = workl[iw + k - 1];
2842 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
2846 else if (!strncmp(type__, "BUCKLE", 6))
2850 for (k = 1; k <= i__1; ++k)
2852 d__2 = workl[iw + k - 1] - 1.;
2853 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
2857 else if (!strncmp(type__, "CAYLEY", 6))
2861 for (k = 1; k <= i__1; ++k)
2863 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2871 if (*rvec && (!strncmp(type__, "SHIFTI", 6) || !strncmp(type__, "CAYLEY", 6)))
2875 for (k = 0; k <= i__1; ++k)
2877 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2881 else if (*rvec && !strncmp(type__, "BUCKLE", 6))
2885 for (k = 0; k <= i__1; ++k)
2887 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
2893 if (strncmp(type__, "REGULR", 6))
2895 F77_FUNC(dger, DGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
2909 /* Selected single precision arpack routines */
2913 F77_FUNC(sstqrb, SSTQRB) (int * n,
2927 int i__, j, k, l, m;
2929 int l1, ii, mm, lm1, mm1, nm1;
2930 float rt1, rt2, eps;
2933 int lend, jtot, lendm1, lendp1, iscale;
2935 int lendsv, nmaxit, icompz;
2936 float ssfmax, ssfmin, safmin, minval, safmax, anorm;
2962 eps = GMX_FLOAT_EPS;
2966 minval = GMX_FLOAT_MIN;
2967 safmin = minval / GMX_FLOAT_EPS;
2968 safmax = 1. / safmin;
2969 ssfmax = sqrt(safmax) / 3.;
2970 ssfmin = sqrt(safmin) / eps2;
2975 for (j = 1; j <= i__1; ++j)
3001 for (m = l1; m <= i__1; ++m)
3008 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps)
3028 i__1 = lend - l + 1;
3029 anorm = F77_FUNC(slanst, SLANST) ("i", &i__1, &d__[l], &e[l]);
3038 i__1 = lend - l + 1;
3039 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
3042 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
3045 else if (anorm < ssfmin)
3048 i__1 = lend - l + 1;
3049 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
3052 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
3056 if (fabs(d__[lend]) < fabs(d__[l]))
3070 for (m = l; m <= i__1; ++m)
3074 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin)
3098 F77_FUNC(slaev2, SLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
3100 work[*n - 1 + l] = s;
3103 z__[l + 1] = c__ * tst - s * z__[l];
3104 z__[l] = s * tst + c__ * z__[l];
3108 F77_FUNC(slae2, SLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
3127 g = (d__[l + 1] - p) / (e[l] * 2.);
3128 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3129 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
3137 for (i__ = mm1; i__ >= i__1; --i__)
3141 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3146 g = d__[i__ + 1] - p;
3147 r__ = (d__[i__] - g) * s + c__ * 2. * b;
3149 d__[i__ + 1] = g + p;
3155 work[*n - 1 + i__] = -s;
3164 F77_FUNC(slasr, SLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
3191 for (m = l; m >= i__1; --m)
3193 d__2 = fabs(e[m - 1]);
3195 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin)
3219 F77_FUNC(slaev2, SLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
3223 z__[l] = c__ * tst - s * z__[l - 1];
3224 z__[l - 1] = s * tst + c__ * z__[l - 1];
3228 F77_FUNC(slae2, SLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
3248 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
3249 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3250 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
3258 for (i__ = m; i__ <= i__1; ++i__)
3262 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3268 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
3276 work[*n - 1 + i__] = s;
3285 F77_FUNC(slasr, SLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
3308 i__1 = lendsv - lsv + 1;
3309 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
3311 i__1 = lendsv - lsv;
3312 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
3315 else if (iscale == 2)
3317 i__1 = lendsv - lsv + 1;
3318 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
3320 i__1 = lendsv - lsv;
3321 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
3330 for (i__ = 1; i__ <= i__1; ++i__)
3343 F77_FUNC(slasrt, SLASRT) ("i", n, &d__[1], info);
3350 for (ii = 2; ii <= i__1; ++ii)
3356 for (j = ii; j <= i__2; ++j)
3382 F77_FUNC(sgetv0, SGETV0) (int * ido,
3384 int gmx_unused * itry,
3401 int v_dim1, v_offset, i__1;
3409 v_offset = 1 + v_dim1;
3425 F77_FUNC(slarnv, SLARNV) (&idist, &iwork[1], n, &resid[1]);
3432 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3451 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
3457 else if (*bmat == 'I')
3459 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3468 workd[*n * 3 + 4] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3469 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
3471 else if (*bmat == 'I')
3473 workd[*n * 3 + 4] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3475 *rnorm = workd[*n * 3 + 4];
3485 F77_FUNC(sgemv, SGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
3486 &workd[*n + 1], &c__1);
3488 F77_FUNC(sgemv, SGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
3489 c_b22, &resid[1], &c__1);
3493 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3499 else if (*bmat == 'I')
3501 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3508 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3509 *rnorm = sqrt(fabs(*rnorm));
3511 else if (*bmat == 'I')
3513 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3516 if (*rnorm > workd[*n * 3 + 4] * .717f)
3525 workd[*n * 3 + 4] = *rnorm;
3532 for (jj = 1; jj <= i__1; ++jj)
3553 F77_FUNC(ssapps, SSAPPS) (int * n,
3570 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
3574 float r__, s, a1, a2, a3, a4;
3585 v_offset = 1 + v_dim1;
3588 h_offset = 1 + h_dim1;
3591 q_offset = 1 + q_dim1;
3594 epsmch = GMX_FLOAT_EPS;
3598 kplusp = *kev + *np;
3600 F77_FUNC(slaset, SLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3608 for (jj = 1; jj <= i__1; ++jj)
3616 for (i__ = istart; i__ <= i__2; ++i__)
3618 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
3619 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3621 h__[i__ + 1 + h_dim1] = 0.;
3632 f = h__[istart + (h_dim1 << 1)] - shift[jj];
3633 g = h__[istart + 1 + h_dim1];
3634 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3636 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
3638 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
3640 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
3642 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
3644 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
3645 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3646 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
3649 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
3650 for (j = 1; j <= i__2; ++j)
3652 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
3654 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
3655 c__ * q[j + (istart + 1) * q_dim1];
3656 q[j + istart * q_dim1] = a1;
3661 for (i__ = istart + 1; i__ <= i__2; ++i__)
3664 f = h__[i__ + h_dim1];
3665 g = s * h__[i__ + 1 + h_dim1];
3667 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3668 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3677 h__[i__ + h_dim1] = r__;
3679 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
3681 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
3683 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
3685 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
3688 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
3689 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3690 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
3693 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
3694 for (j = 1; j <= i__3; ++j)
3696 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
3698 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
3699 c__ * q[j + (i__ + 1) * q_dim1];
3700 q[j + i__ * q_dim1] = a1;
3709 if (h__[iend + h_dim1] < 0.)
3711 h__[iend + h_dim1] = -h__[iend + h_dim1];
3712 F77_FUNC(sscal, SSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3721 for (i__ = itop; i__ <= i__2; ++i__)
3723 if (h__[i__ + 1 + h_dim1] > 0.)
3735 for (i__ = itop; i__ <= i__1; ++i__)
3737 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
3738 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3740 h__[i__ + 1 + h_dim1] = 0.;
3745 if (h__[*kev + 1 + h_dim1] > 0.)
3747 F77_FUNC(sgemv, SGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
3748 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
3752 for (i__ = 1; i__ <= i__1; ++i__)
3754 i__2 = kplusp - i__ + 1;
3755 F77_FUNC(sgemv, SGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
3756 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
3757 F77_FUNC(scopy, SCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
3762 F77_FUNC(slacpy, SLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3764 if (h__[*kev + 1 + h_dim1] > 0.)
3766 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3769 F77_FUNC(sscal, SSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3770 if (h__[*kev + 1 + h_dim1] > 0.)
3772 F77_FUNC(saxpy, SAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
3786 F77_FUNC(ssortr, SSORTR) (const char * which,
3801 if (!strncmp(which, "SA", 2))
3810 for (i__ = igap; i__ <= i__1; ++i__)
3820 if (x1[j] < x1[j + igap])
3823 x1[j] = x1[j + igap];
3824 x1[j + igap] = temp;
3828 x2[j] = x2[j + igap];
3829 x2[j + igap] = temp;
3845 else if (!strncmp(which, "SM", 2))
3854 for (i__ = igap; i__ <= i__1; ++i__)
3864 if (fabs(x1[j]) < fabs(x1[j + igap]))
3867 x1[j] = x1[j + igap];
3868 x1[j + igap] = temp;
3872 x2[j] = x2[j + igap];
3873 x2[j + igap] = temp;
3889 else if (!strncmp(which, "LA", 2))
3898 for (i__ = igap; i__ <= i__1; ++i__)
3908 if (x1[j] > x1[j + igap])
3911 x1[j] = x1[j + igap];
3912 x1[j + igap] = temp;
3916 x2[j] = x2[j + igap];
3917 x2[j + igap] = temp;
3933 else if (!strncmp(which, "LM", 2))
3943 for (i__ = igap; i__ <= i__1; ++i__)
3953 if (fabs(x1[j]) > fabs(x1[j + igap]))
3956 x1[j] = x1[j + igap];
3957 x1[j + igap] = temp;
3961 x2[j] = x2[j + igap];
3962 x2[j + igap] = temp;
3987 F77_FUNC(ssesrt, SSESRT) (const char * which,
3995 int a_dim1, a_offset, i__1;
4002 a_offset = 1 + a_dim1 * 0;
4007 if (!strncmp(which, "SA", 2))
4016 for (i__ = igap; i__ <= i__1; ++i__)
4026 if (x[j] < x[j + igap])
4033 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4034 a_dim1 + 1], &c__1);
4050 else if (!strncmp(which, "SM", 2))
4059 for (i__ = igap; i__ <= i__1; ++i__)
4069 if (fabs(x[j]) < fabs(x[j + igap]))
4076 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4077 a_dim1 + 1], &c__1);
4093 else if (!strncmp(which, "LA", 2))
4102 for (i__ = igap; i__ <= i__1; ++i__)
4112 if (x[j] > x[j + igap])
4119 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4120 a_dim1 + 1], &c__1);
4136 else if (!strncmp(which, "LM", 2))
4145 for (i__ = igap; i__ <= i__1; ++i__)
4155 if (fabs(x[j]) > fabs(x[j + igap]))
4162 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4163 a_dim1 + 1], &c__1);
4188 F77_FUNC(ssgets, SSGETS) (int * ishift,
4204 if (!strncmp(which, "BE", 2))
4207 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
4211 i__1 = (kevd2 < *np) ? kevd2 : *np;
4212 i__2 = (kevd2 > *np) ? kevd2 : *np;
4213 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[1], &c__1,
4214 &ritz[i__2 + 1], &c__1);
4215 i__1 = (kevd2 < *np) ? kevd2 : *np;
4216 i__2 = (kevd2 > *np) ? kevd2 : *np;
4217 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[1], &c__1,
4218 &bounds[i__2 + 1], &c__1);
4225 F77_FUNC(ssortr, SSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
4228 if (*ishift == 1 && *np > 0)
4231 F77_FUNC(ssortr, SSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
4232 F77_FUNC(scopy, SCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
4242 F77_FUNC(ssconv, SSCONV) (int * n,
4258 eps23 = GMX_FLOAT_EPS;
4259 eps23 = pow(eps23, c_b3);
4263 for (i__ = 1; i__ <= i__1; ++i__)
4267 d__3 = fabs(ritz[i__]);
4268 temp = (d__2 > d__3) ? d__2 : d__3;
4269 if (bounds[i__] <= *tol * temp)
4280 F77_FUNC(sseigt, SSEIGT) (float * rnorm,
4290 int h_dim1, h_offset, i__1;
4299 h_offset = 1 + h_dim1;
4302 F77_FUNC(scopy, SCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
4304 F77_FUNC(scopy, SCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
4305 F77_FUNC(sstqrb, SSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
4312 for (k = 1; k <= i__1; ++k)
4314 bounds[k] = *rnorm * fabs(bounds[k]);
4327 F77_FUNC(ssaitr, SSAITR) (int * ido,
4351 int h_dim1, h_offset, v_dim1, v_offset, i__1;
4355 float safmin, minval;
4361 v_offset = 1 + v_dim1;
4364 h_offset = 1 + h_dim1;
4368 minval = GMX_FLOAT_MIN;
4369 safmin = minval / GMX_FLOAT_EPS;
4383 iwork[9] = iwork[8] + *n;
4384 iwork[10] = iwork[9] + *n;
4422 F77_FUNC(sgetv0, sgetv0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
4423 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
4436 *info = iwork[12] - 1;
4443 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
4444 if (*rnorm >= safmin)
4446 temp1 = 1. / *rnorm;
4447 F77_FUNC(sscal, SSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
4448 F77_FUNC(sscal, SSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
4453 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
4454 v_dim1 + 1], n, &infol);
4455 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
4460 F77_FUNC(scopy, SCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
4461 ipntr[1] = iwork[10];
4462 ipntr[2] = iwork[9];
4463 ipntr[3] = iwork[8];
4472 F77_FUNC(scopy, SCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
4481 ipntr[1] = iwork[9];
4482 ipntr[2] = iwork[8];
4487 else if (*bmat == 'I')
4489 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4499 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
4501 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
4503 else if (*bmat == 'G')
4505 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4507 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
4509 else if (*bmat == 'I')
4511 workd[*n * 3 + 3] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4516 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
4517 &c__1, &c_b42, &workd[iwork[9]], &c__1);
4519 else if (*mode == 2)
4521 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
4522 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
4525 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4526 c__1, &c_b18, &resid[1], &c__1);
4528 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
4529 if (iwork[12] == 1 || iwork[4] == 1)
4531 h__[iwork[12] + h_dim1] = 0.;
4535 h__[iwork[12] + h_dim1] = *rnorm;
4543 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4544 ipntr[1] = iwork[9];
4545 ipntr[2] = iwork[8];
4550 else if (*bmat == 'I')
4552 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4560 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4561 *rnorm = sqrt(fabs(*rnorm));
4563 else if (*bmat == 'I')
4565 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4568 if (*rnorm > workd[*n * 3 + 3] * .717f)
4575 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
4576 c__1, &c_b42, &workd[iwork[9]], &c__1);
4578 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4579 c__1, &c_b18, &resid[1], &c__1);
4581 if (iwork[12] == 1 || iwork[4] == 1)
4583 h__[iwork[12] + h_dim1] = 0.;
4585 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
4590 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4591 ipntr[1] = iwork[9];
4592 ipntr[2] = iwork[8];
4597 else if (*bmat == 'I')
4599 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4606 workd[*n * 3 + 2] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4608 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
4610 else if (*bmat == 'I')
4612 workd[*n * 3 + 2] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4616 if (workd[*n * 3 + 2] > *rnorm * .717f)
4619 *rnorm = workd[*n * 3 + 2];
4625 *rnorm = workd[*n * 3 + 2];
4633 for (jj = 1; jj <= i__1; ++jj)
4645 if (h__[iwork[12] + h_dim1] < 0.)
4647 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4648 if (iwork[12] < *k + *np)
4650 F77_FUNC(sscal, SSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4654 F77_FUNC(sscal, SSCAL) (n, &c_b50, &resid[1], &c__1);
4659 if (iwork[12] > *k + *np)
4679 F77_FUNC(ssaup2, SSAUP2) (int * ido,
4688 int gmx_unused * iupd,
4709 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
4728 v_offset = 1 + v_dim1;
4731 h_offset = 1 + h_dim1;
4734 q_offset = 1 + q_dim1;
4738 eps23 = GMX_FLOAT_EPS;
4739 eps23 = pow(eps23, c_b3);
4752 iwork[7] = iwork[9] + iwork[10];
4775 F77_FUNC(sgetv0, SGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
4776 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
4784 if (workd[*n * 3 + 1] == 0.)
4809 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
4810 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
4836 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[
4837 v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[
4855 F77_FUNC(sseigt, SSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
4856 bounds[1], &workl[1], &ierr);
4864 F77_FUNC(scopy, SCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4865 F77_FUNC(scopy, SCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4869 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4871 F77_FUNC(scopy, SCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4872 F77_FUNC(ssconv, SSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4877 for (j = 1; j <= i__1; ++j)
4879 if (bounds[j] == 0.)
4886 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
4889 if (!strncmp(which, "BE", 2))
4892 strncpy(wprime, "SA", 2);
4893 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4895 nevm2 = *nev - nevd2;
4898 i__1 = (nevd2 < *np) ? nevd2 : *np;
4899 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4900 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
4901 &ritz[((i__2 > i__3) ? i__2 : i__3)],
4903 i__1 = (nevd2 < *np) ? nevd2 : *np;
4904 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4905 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
4906 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
4914 if (!strncmp(which, "LM", 2))
4916 strncpy(wprime, "SM", 2);
4918 if (!strncmp(which, "SM", 2))
4920 strncpy(wprime, "LM", 2);
4922 if (!strncmp(which, "LA", 2))
4924 strncpy(wprime, "SA", 2);
4926 if (!strncmp(which, "SA", 2))
4928 strncpy(wprime, "LA", 2);
4931 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4936 for (j = 1; j <= i__1; ++j)
4939 d__3 = fabs(ritz[j]);
4940 temp = (d__2 > d__3) ? d__2 : d__3;
4944 strncpy(wprime, "LA", 2);
4945 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4948 for (j = 1; j <= i__1; ++j)
4951 d__3 = fabs(ritz[j]);
4952 temp = (d__2 > d__3) ? d__2 : d__3;
4956 if (!strncmp(which, "BE", 2))
4959 strncpy(wprime, "LA", 2);
4960 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4965 F77_FUNC(ssortr, SSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4969 h__[h_dim1 + 1] = workd[*n * 3 + 1];
4972 if (iwork[6] > *mxiter && iwork[8] < *nev)
4976 if (*np == 0 && iwork[8] < iwork[9])
4985 else if (iwork[8] < *nev && *ishift == 1)
4988 i__1 = iwork[8], i__2 = *np / 2;
4989 *nev += (i__1 < i__2) ? i__1 : i__2;
4990 if (*nev == 1 && iwork[7] >= 6)
4992 *nev = iwork[7] / 2;
4994 else if (*nev == 1 && iwork[7] > 2)
4998 *np = iwork[7] - *nev;
5003 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
5023 F77_FUNC(scopy, SCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
5026 F77_FUNC(ssapps, SSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
5027 resid[1], &q[q_offset], ldq, &workd[1]);
5032 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
5039 else if (*bmat == 'I')
5041 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
5048 workd[*n * 3 + 1] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
5049 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
5051 else if (*bmat == 'I')
5053 workd[*n * 3 + 1] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
5075 F77_FUNC(ssaupd, SSAUPD) (int * ido,
5093 int v_dim1, v_offset, i__1, i__2;
5099 v_offset = 1 + v_dim1;
5111 iwork[5] = iparam[1];
5112 iwork[10] = iparam[3];
5113 iwork[12] = iparam[4];
5116 iwork[11] = iparam[7];
5127 else if (*ncv <= *nev || *ncv > *n)
5133 iwork[15] = *ncv - *nev;
5139 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
5140 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
5141 strncmp(which, "BE", 2))
5145 if (*bmat != 'I' && *bmat != 'G')
5151 if (*lworkl < i__1 * i__1 + (*ncv << 3))
5155 if (iwork[11] < 1 || iwork[11] > 5)
5159 else if (iwork[11] == 1 && *bmat == 'G')
5163 else if (iwork[5] < 0 || iwork[5] > 1)
5167 else if (*nev == 1 && !strncmp(which, "BE", 2))
5185 *tol = GMX_FLOAT_EPS;
5188 iwork[15] = *ncv - *nev;
5191 i__1 = i__2 * i__2 + (*ncv << 3);
5192 for (j = 1; j <= i__1; ++j)
5200 iwork[16] = iwork[3] + (iwork[8] << 1);
5201 iwork[1] = iwork[16] + *ncv;
5202 iwork[4] = iwork[1] + *ncv;
5204 iwork[7] = iwork[4] + i__1 * i__1;
5205 iwork[14] = iwork[7] + *ncv * 3;
5207 ipntr[4] = iwork[14];
5208 ipntr[5] = iwork[3];
5209 ipntr[6] = iwork[16];
5210 ipntr[7] = iwork[1];
5211 ipntr[11] = iwork[7];
5214 F77_FUNC(ssaup2, SSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
5215 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
5216 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
5217 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
5222 iparam[8] = iwork[15];
5229 iparam[3] = iwork[10];
5230 iparam[5] = iwork[15];
5250 F77_FUNC(sseupd, SSEUPD) (int * rvec,
5251 const char * howmny,
5276 int v_dim1, v_offset, z_dim1, z_offset, i__1;
5277 float d__1, d__2, d__3;
5279 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
5291 float thres1 = 0, thres2 = 0;
5295 int leftptr, rghtptr;
5301 z_offset = 1 + z_dim1;
5306 v_offset = 1 + v_dim1;
5334 if (*ncv <= *nev || *ncv > *n)
5338 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
5339 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
5340 strncmp(which, "BE", 2))
5344 if (*bmat != 'I' && *bmat != 'G')
5348 if (*howmny != 'A' && *howmny != 'P' &&
5349 *howmny != 'S' && *rvec)
5353 if (*rvec && *howmny == 'S')
5358 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
5363 if (mode == 1 || mode == 2)
5365 strncpy(type__, "REGULR", 6);
5369 strncpy(type__, "SHIFTI", 6);
5373 strncpy(type__, "BUCKLE", 6);
5377 strncpy(type__, "CAYLEY", 6);
5383 if (mode == 1 && *bmat == 'G')
5387 if (*nev == 1 && !strncmp(which, "BE", 2))
5406 iw = iq + ldh * *ncv;
5407 next = iw + (*ncv << 1);
5413 irz = ipntr[11] + *ncv;
5417 eps23 = GMX_FLOAT_EPS;
5418 eps23 = pow(eps23, c_b21);
5425 else if (*bmat == 'G')
5427 bnorm2 = F77_FUNC(snrm2, SNRM2) (n, &workd[1], &c__1);
5433 if (!strncmp(which, "LM", 2) || !strncmp(which, "SM", 2) ||
5434 !strncmp(which, "LA", 2) || !strncmp(which, "SA", 2))
5438 else if (!strncmp(which, "BE", 2))
5442 ism = (*nev > nconv) ? *nev : nconv;
5445 thres1 = workl[ism];
5446 thres2 = workl[ilg];
5454 for (j = 0; j <= i__1; ++j)
5457 if (!strncmp(which, "LM", 2))
5459 if (fabs(workl[irz + j]) >= fabs(thres1))
5462 d__3 = fabs(workl[irz + j]);
5463 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5464 if (workl[ibd + j] <= *tol * tempbnd)
5470 else if (!strncmp(which, "SM", 2))
5472 if (fabs(workl[irz + j]) <= fabs(thres1))
5475 d__3 = fabs(workl[irz + j]);
5476 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5477 if (workl[ibd + j] <= *tol * tempbnd)
5483 else if (!strncmp(which, "LA", 2))
5485 if (workl[irz + j] >= thres1)
5488 d__3 = fabs(workl[irz + j]);
5489 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5490 if (workl[ibd + j] <= *tol * tempbnd)
5496 else if (!strncmp(which, "SA", 2))
5498 if (workl[irz + j] <= thres1)
5501 d__3 = fabs(workl[irz + j]);
5502 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5503 if (workl[ibd + j] <= *tol * tempbnd)
5509 else if (!strncmp(which, "BE", 2))
5511 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
5514 d__3 = fabs(workl[irz + j]);
5515 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5516 if (workl[ibd + j] <= *tol * tempbnd)
5524 reord = select[j + 1] || reord;
5533 F77_FUNC(scopy, SCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
5534 F77_FUNC(scopy, SCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
5536 F77_FUNC(ssteqr, SSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
5558 if (select[leftptr])
5564 else if (!select[rghtptr])
5573 temp = workl[ihd + leftptr - 1];
5574 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
5575 workl[ihd + rghtptr - 1] = temp;
5576 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
5578 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
5579 iq + *ncv * (leftptr - 1)], &c__1);
5580 F77_FUNC(scopy, SCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
5587 if (leftptr < rghtptr)
5596 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5602 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
5603 F77_FUNC(scopy, SCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
5606 if (!strncmp(type__, "REGULR", 6))
5611 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5615 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5622 F77_FUNC(scopy, SCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
5623 if (!strncmp(type__, "SHIFTI", 6))
5626 for (k = 1; k <= i__1; ++k)
5628 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
5631 else if (!strncmp(type__, "BUCKLE", 6))
5634 for (k = 1; k <= i__1; ++k)
5636 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
5640 else if (!strncmp(type__, "CAYLEY", 6))
5643 for (k = 1; k <= i__1; ++k)
5645 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
5646 workl[ihd + k - 1] - 1.);
5650 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5651 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
5654 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5658 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5659 d__1 = bnorm2 / rnorm;
5660 F77_FUNC(sscal, SSCAL) (ncv, &d__1, &workl[ihb], &c__1);
5661 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
5666 if (*rvec && *howmny == 'A')
5669 F77_FUNC(sgeqr2, SGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
5672 F77_FUNC(sorm2r, SORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
5673 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
5674 F77_FUNC(slacpy, SLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
5677 for (j = 1; j <= i__1; ++j)
5679 workl[ihb + j - 1] = 0.;
5681 workl[ihb + *ncv - 1] = 1.;
5682 F77_FUNC(sorm2r, SORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
5683 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
5686 else if (*rvec && *howmny == 'S')
5691 if (!strncmp(type__, "REGULR", 6) && *rvec)
5695 for (j = 1; j <= i__1; ++j)
5697 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
5701 else if (strncmp(type__, "REGULR", 6) && *rvec)
5704 F77_FUNC(sscal, SSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
5705 if (!strncmp(type__, "SHIFTI", 6))
5709 for (k = 1; k <= i__1; ++k)
5711 d__2 = workl[iw + k - 1];
5712 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
5716 else if (!strncmp(type__, "BUCKLE", 6))
5720 for (k = 1; k <= i__1; ++k)
5722 d__2 = workl[iw + k - 1] - 1.;
5723 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
5727 else if (!strncmp(type__, "CAYLEY", 6))
5731 for (k = 1; k <= i__1; ++k)
5733 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
5741 if (*rvec && (!strncmp(type__, "SHIFTI", 6) || !strncmp(type__, "CAYLEY", 6)))
5745 for (k = 0; k <= i__1; ++k)
5747 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
5751 else if (*rvec && !strncmp(type__, "BUCKLE", 6))
5755 for (k = 0; k <= i__1; ++k)
5757 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
5763 if (strncmp(type__, "REGULR", 6))
5765 F77_FUNC(sger, SGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[