1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: m; c-basic-offset: 4 -*-
4 * This file is part of Gromacs Copyright (c) 1991-2004
5 * David van der Spoel, Erik Lindahl, University of Groningen.
7 * This file contains a subset of ARPACK functions to perform
8 * diagonalization and SVD for sparse matrices in Gromacs.
10 * The code has been translated to C to avoid being dependent on
11 * a Fotran compiler, and it has been made threadsafe by using
12 * additional workspace arrays to store data during reverse communication.
14 * You might prefer the original ARPACK library for general use, but
15 * in case you want to this version can be redistributed freely, just
16 * as the original library. However, please make clear that it is the
17 * hacked version from Gromacs so any bugs are blamed on us and not
18 * the original authors. You should also be aware that the double
19 * precision work array workd needs to be of size (3*N+4) here
20 * (4 more than the general library), and there is an extra argument
21 * iwork, which should be an integer work array of length 80.
23 * ARPACK was written by
25 * Danny Sorensen Phuong Vu
26 * Rihard Lehoucq CRPC / Rice University
27 * Dept. of Computational & Houston, Texas
35 #include "gromacs/legacyheaders/types/simple.h"
36 #include "gmx_arpack.h"
38 #include "gmx_lapack.h"
40 F77_FUNC(dstqrb, DSTQRB) (int * n,
56 int l1, ii, mm, lm1, mm1, nm1;
60 int lend, jtot, lendm1, lendp1, iscale;
62 int lendsv, nmaxit, icompz;
63 double ssfmax, ssfmin, safmin, minval, safmax, anorm;
93 minval = GMX_DOUBLE_MIN;
94 safmin = minval / GMX_DOUBLE_EPS;
96 ssfmax = sqrt(safmax) / 3.;
97 ssfmin = sqrt(safmin) / eps2;
102 for (j = 1; j <= i__1; ++j)
128 for (m = l1; m <= i__1; ++m)
135 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps)
156 anorm = F77_FUNC(dlanst, DLANST) ("i", &i__1, &d__[l], &e[l]);
166 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
169 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
172 else if (anorm < ssfmin)
176 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
179 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
183 if (fabs(d__[lend]) < fabs(d__[l]))
197 for (m = l; m <= i__1; ++m)
201 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin)
225 F77_FUNC(dlaev2, DLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
227 work[*n - 1 + l] = s;
230 z__[l + 1] = c__ * tst - s * z__[l];
231 z__[l] = s * tst + c__ * z__[l];
235 F77_FUNC(dlae2, DLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
254 g = (d__[l + 1] - p) / (e[l] * 2.);
255 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
256 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
264 for (i__ = mm1; i__ >= i__1; --i__)
268 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
273 g = d__[i__ + 1] - p;
274 r__ = (d__[i__] - g) * s + c__ * 2. * b;
276 d__[i__ + 1] = g + p;
282 work[*n - 1 + i__] = -s;
291 F77_FUNC(dlasr, DLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
318 for (m = l; m >= i__1; --m)
320 d__2 = fabs(e[m - 1]);
322 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin)
346 F77_FUNC(dlaev2, DLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
350 z__[l] = c__ * tst - s * z__[l - 1];
351 z__[l - 1] = s * tst + c__ * z__[l - 1];
355 F77_FUNC(dlae2, DLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
375 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
376 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
377 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
385 for (i__ = m; i__ <= i__1; ++i__)
389 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
395 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
403 work[*n - 1 + i__] = s;
412 F77_FUNC(dlasr, DLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
435 i__1 = lendsv - lsv + 1;
436 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
439 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
442 else if (iscale == 2)
444 i__1 = lendsv - lsv + 1;
445 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
448 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
457 for (i__ = 1; i__ <= i__1; ++i__)
470 F77_FUNC(dlasrt, DLASRT) ("i", n, &d__[1], info);
477 for (ii = 2; ii <= i__1; ++ii)
483 for (j = ii; j <= i__2; ++j)
509 F77_FUNC(dgetv0, DGETV0) (int * ido,
528 int v_dim1, v_offset, i__1;
536 v_offset = 1 + v_dim1;
552 F77_FUNC(dlarnv, DLARNV) (&idist, &iwork[1], n, &resid[1]);
559 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
578 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
584 else if (*bmat == 'I')
586 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
595 workd[*n * 3 + 4] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
596 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
598 else if (*bmat == 'I')
600 workd[*n * 3 + 4] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
602 *rnorm = workd[*n * 3 + 4];
612 F77_FUNC(dgemv, DGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
613 &workd[*n + 1], &c__1);
615 F77_FUNC(dgemv, DGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
616 c_b22, &resid[1], &c__1);
620 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
626 else if (*bmat == 'I')
628 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
635 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
636 *rnorm = sqrt(fabs(*rnorm));
638 else if (*bmat == 'I')
640 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
643 if (*rnorm > workd[*n * 3 + 4] * .717f)
652 workd[*n * 3 + 4] = *rnorm;
659 for (jj = 1; jj <= i__1; ++jj)
680 F77_FUNC(dsapps, DSAPPS) (int * n,
697 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
701 double r__, s, a1, a2, a3, a4;
712 v_offset = 1 + v_dim1;
715 h_offset = 1 + h_dim1;
718 q_offset = 1 + q_dim1;
721 epsmch = GMX_DOUBLE_EPS;
727 F77_FUNC(dlaset, DLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
735 for (jj = 1; jj <= i__1; ++jj)
743 for (i__ = istart; i__ <= i__2; ++i__)
745 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
746 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
748 h__[i__ + 1 + h_dim1] = 0.;
759 f = h__[istart + (h_dim1 << 1)] - shift[jj];
760 g = h__[istart + 1 + h_dim1];
761 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
763 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
765 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
767 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
769 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
771 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
772 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
773 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
776 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
777 for (j = 1; j <= i__2; ++j)
779 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
781 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
782 c__ * q[j + (istart + 1) * q_dim1];
783 q[j + istart * q_dim1] = a1;
788 for (i__ = istart + 1; i__ <= i__2; ++i__)
791 f = h__[i__ + h_dim1];
792 g = s * h__[i__ + 1 + h_dim1];
794 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
795 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
804 h__[i__ + h_dim1] = r__;
806 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
808 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
810 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
812 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
815 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
816 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
817 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
820 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
821 for (j = 1; j <= i__3; ++j)
823 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
825 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
826 c__ * q[j + (i__ + 1) * q_dim1];
827 q[j + i__ * q_dim1] = a1;
836 if (h__[iend + h_dim1] < 0.)
838 h__[iend + h_dim1] = -h__[iend + h_dim1];
839 F77_FUNC(dscal, DSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
848 for (i__ = itop; i__ <= i__2; ++i__)
850 if (h__[i__ + 1 + h_dim1] > 0.)
862 for (i__ = itop; i__ <= i__1; ++i__)
864 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
865 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
867 h__[i__ + 1 + h_dim1] = 0.;
872 if (h__[*kev + 1 + h_dim1] > 0.)
874 F77_FUNC(dgemv, DGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
875 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
879 for (i__ = 1; i__ <= i__1; ++i__)
881 i__2 = kplusp - i__ + 1;
882 F77_FUNC(dgemv, DGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
883 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
884 F77_FUNC(dcopy, DCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
889 F77_FUNC(dlacpy, DLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
891 if (h__[*kev + 1 + h_dim1] > 0.)
893 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
896 F77_FUNC(dscal, DSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
897 if (h__[*kev + 1 + h_dim1] > 0.)
899 F77_FUNC(daxpy, DAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
913 F77_FUNC(dsortr, DSORTR) (const char * which,
928 if (!strncmp(which, "SA", 2))
937 for (i__ = igap; i__ <= i__1; ++i__)
947 if (x1[j] < x1[j + igap])
950 x1[j] = x1[j + igap];
955 x2[j] = x2[j + igap];
972 else if (!strncmp(which, "SM", 2))
981 for (i__ = igap; i__ <= i__1; ++i__)
991 if (fabs(x1[j]) < fabs(x1[j + igap]))
994 x1[j] = x1[j + igap];
999 x2[j] = x2[j + igap];
1000 x2[j + igap] = temp;
1016 else if (!strncmp(which, "LA", 2))
1025 for (i__ = igap; i__ <= i__1; ++i__)
1035 if (x1[j] > x1[j + igap])
1038 x1[j] = x1[j + igap];
1039 x1[j + igap] = temp;
1043 x2[j] = x2[j + igap];
1044 x2[j + igap] = temp;
1060 else if (!strncmp(which, "LM", 2))
1070 for (i__ = igap; i__ <= i__1; ++i__)
1080 if (fabs(x1[j]) > fabs(x1[j + igap]))
1083 x1[j] = x1[j + igap];
1084 x1[j + igap] = temp;
1088 x2[j] = x2[j + igap];
1089 x2[j + igap] = temp;
1114 F77_FUNC(dsesrt, DSESRT) (const char * which,
1122 int a_dim1, a_offset, i__1;
1129 a_offset = 1 + a_dim1 * 0;
1134 if (!strncmp(which, "SA", 2))
1143 for (i__ = igap; i__ <= i__1; ++i__)
1153 if (x[j] < x[j + igap])
1160 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1161 a_dim1 + 1], &c__1);
1177 else if (!strncmp(which, "SM", 2))
1186 for (i__ = igap; i__ <= i__1; ++i__)
1196 if (fabs(x[j]) < fabs(x[j + igap]))
1203 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1204 a_dim1 + 1], &c__1);
1220 else if (!strncmp(which, "LA", 2))
1229 for (i__ = igap; i__ <= i__1; ++i__)
1239 if (x[j] > x[j + igap])
1246 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1247 a_dim1 + 1], &c__1);
1263 else if (!strncmp(which, "LM", 2))
1272 for (i__ = igap; i__ <= i__1; ++i__)
1282 if (fabs(x[j]) > fabs(x[j + igap]))
1289 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1290 a_dim1 + 1], &c__1);
1315 F77_FUNC(dsgets, DSGETS) (int * ishift,
1331 if (!strncmp(which, "BE", 2))
1334 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1338 i__1 = (kevd2 < *np) ? kevd2 : *np;
1339 i__2 = (kevd2 > *np) ? kevd2 : *np;
1340 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[1], &c__1,
1341 &ritz[i__2 + 1], &c__1);
1342 i__1 = (kevd2 < *np) ? kevd2 : *np;
1343 i__2 = (kevd2 > *np) ? kevd2 : *np;
1344 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[1], &c__1,
1345 &bounds[i__2 + 1], &c__1);
1352 F77_FUNC(dsortr, DSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
1355 if (*ishift == 1 && *np > 0)
1358 F77_FUNC(dsortr, DSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
1359 F77_FUNC(dcopy, DCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
1369 F77_FUNC(dsconv, DSCONV) (int * n,
1385 eps23 = GMX_DOUBLE_EPS;
1386 eps23 = pow(eps23, c_b3);
1390 for (i__ = 1; i__ <= i__1; ++i__)
1394 d__3 = fabs(ritz[i__]);
1395 temp = (d__2 > d__3) ? d__2 : d__3;
1396 if (bounds[i__] <= *tol * temp)
1407 F77_FUNC(dseigt, DSEIGT) (double * rnorm,
1417 int h_dim1, h_offset, i__1;
1426 h_offset = 1 + h_dim1;
1429 F77_FUNC(dcopy, DCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1431 F77_FUNC(dcopy, DCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1432 F77_FUNC(dstqrb, DSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1439 for (k = 1; k <= i__1; ++k)
1441 bounds[k] = *rnorm * fabs(bounds[k]);
1454 F77_FUNC(dsaitr, DSAITR) (int * ido,
1478 int h_dim1, h_offset, v_dim1, v_offset, i__1;
1482 double safmin, minval;
1488 v_offset = 1 + v_dim1;
1491 h_offset = 1 + h_dim1;
1495 minval = GMX_DOUBLE_MIN;
1496 safmin = minval / GMX_DOUBLE_EPS;
1510 iwork[9] = iwork[8] + *n;
1511 iwork[10] = iwork[9] + *n;
1549 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
1550 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
1563 *info = iwork[12] - 1;
1570 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1571 if (*rnorm >= safmin)
1573 temp1 = 1. / *rnorm;
1574 F77_FUNC(dscal, DSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1575 F77_FUNC(dscal, DSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
1580 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
1581 v_dim1 + 1], n, &infol);
1582 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
1587 F77_FUNC(dcopy, DCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1588 ipntr[1] = iwork[10];
1589 ipntr[2] = iwork[9];
1590 ipntr[3] = iwork[8];
1599 F77_FUNC(dcopy, DCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1608 ipntr[1] = iwork[9];
1609 ipntr[2] = iwork[8];
1614 else if (*bmat == 'I')
1616 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1626 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
1628 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1630 else if (*bmat == 'G')
1632 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1634 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1636 else if (*bmat == 'I')
1638 workd[*n * 3 + 3] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1643 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
1644 &c__1, &c_b42, &workd[iwork[9]], &c__1);
1646 else if (*mode == 2)
1648 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
1649 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
1652 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1653 c__1, &c_b18, &resid[1], &c__1);
1655 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1656 if (iwork[12] == 1 || iwork[4] == 1)
1658 h__[iwork[12] + h_dim1] = 0.;
1662 h__[iwork[12] + h_dim1] = *rnorm;
1670 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1671 ipntr[1] = iwork[9];
1672 ipntr[2] = iwork[8];
1677 else if (*bmat == 'I')
1679 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1687 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1688 *rnorm = sqrt(fabs(*rnorm));
1690 else if (*bmat == 'I')
1692 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1695 if (*rnorm > workd[*n * 3 + 3] * .717f)
1702 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
1703 c__1, &c_b42, &workd[iwork[9]], &c__1);
1705 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1706 c__1, &c_b18, &resid[1], &c__1);
1708 if (iwork[12] == 1 || iwork[4] == 1)
1710 h__[iwork[12] + h_dim1] = 0.;
1712 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1717 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1718 ipntr[1] = iwork[9];
1719 ipntr[2] = iwork[8];
1724 else if (*bmat == 'I')
1726 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1733 workd[*n * 3 + 2] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1735 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
1737 else if (*bmat == 'I')
1739 workd[*n * 3 + 2] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1743 if (workd[*n * 3 + 2] > *rnorm * .717f)
1746 *rnorm = workd[*n * 3 + 2];
1752 *rnorm = workd[*n * 3 + 2];
1760 for (jj = 1; jj <= i__1; ++jj)
1772 if (h__[iwork[12] + h_dim1] < 0.)
1774 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1775 if (iwork[12] < *k + *np)
1777 F77_FUNC(dscal, DSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1781 F77_FUNC(dscal, DSCAL) (n, &c_b50, &resid[1], &c__1);
1786 if (iwork[12] > *k + *np)
1806 F77_FUNC(dsaup2, DSAUP2) (int * ido,
1836 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
1855 v_offset = 1 + v_dim1;
1858 h_offset = 1 + h_dim1;
1861 q_offset = 1 + q_dim1;
1865 eps23 = GMX_DOUBLE_EPS;
1866 eps23 = pow(eps23, c_b3);
1879 iwork[7] = iwork[9] + iwork[10];
1902 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
1903 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
1911 if (workd[*n * 3 + 1] == 0.)
1936 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
1937 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1963 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1],
1964 &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1982 F77_FUNC(dseigt, DSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
1983 bounds[1], &workl[1], &ierr);
1991 F77_FUNC(dcopy, DCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1992 F77_FUNC(dcopy, DCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
1996 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
1998 F77_FUNC(dcopy, DCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
1999 F77_FUNC(dsconv, DSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
2003 for (j = 1; j <= i__1; ++j)
2005 if (bounds[j] == 0.)
2012 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
2015 if (!strncmp(which, "BE", 2))
2018 strncpy(wprime, "SA", 2);
2019 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2021 nevm2 = *nev - nevd2;
2024 i__1 = (nevd2 < *np) ? nevd2 : *np;
2025 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
2026 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
2027 &ritz[((i__2 > i__3) ? i__2 : i__3)],
2029 i__1 = (nevd2 < *np) ? nevd2 : *np;
2030 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
2031 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
2032 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
2040 if (!strncmp(which, "LM", 2))
2042 strncpy(wprime, "SM", 2);
2044 if (!strncmp(which, "SM", 2))
2046 strncpy(wprime, "LM", 2);
2048 if (!strncmp(which, "LA", 2))
2050 strncpy(wprime, "SA", 2);
2052 if (!strncmp(which, "SA", 2))
2054 strncpy(wprime, "LA", 2);
2057 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2062 for (j = 1; j <= i__1; ++j)
2065 d__3 = fabs(ritz[j]);
2066 temp = (d__2 > d__3) ? d__2 : d__3;
2070 strncpy(wprime, "LA", 2);
2071 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
2074 for (j = 1; j <= i__1; ++j)
2077 d__3 = fabs(ritz[j]);
2078 temp = (d__2 > d__3) ? d__2 : d__3;
2082 if (!strncmp(which, "BE", 2))
2085 strncpy(wprime, "LA", 2);
2086 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2091 F77_FUNC(dsortr, DSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2095 h__[h_dim1 + 1] = workd[*n * 3 + 1];
2098 if (iwork[6] > *mxiter && iwork[8] < *nev)
2102 if (*np == 0 && iwork[8] < iwork[9])
2111 else if (iwork[8] < *nev && *ishift == 1)
2114 i__1 = iwork[8], i__2 = *np / 2;
2115 *nev += (i__1 < i__2) ? i__1 : i__2;
2116 if (*nev == 1 && iwork[7] >= 6)
2118 *nev = iwork[7] / 2;
2120 else if (*nev == 1 && iwork[7] > 2)
2124 *np = iwork[7] - *nev;
2129 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2149 F77_FUNC(dcopy, DCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
2152 F77_FUNC(dsapps, DSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
2153 resid[1], &q[q_offset], ldq, &workd[1]);
2158 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
2165 else if (*bmat == 'I')
2167 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
2174 workd[*n * 3 + 1] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
2175 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
2177 else if (*bmat == 'I')
2179 workd[*n * 3 + 1] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
2201 F77_FUNC(dsaupd, DSAUPD) (int * ido,
2219 int v_dim1, v_offset, i__1, i__2;
2225 v_offset = 1 + v_dim1;
2237 iwork[5] = iparam[1];
2238 iwork[10] = iparam[3];
2239 iwork[12] = iparam[4];
2242 iwork[11] = iparam[7];
2253 else if (*ncv <= *nev || *ncv > *n)
2259 iwork[15] = *ncv - *nev;
2265 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
2266 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
2267 strncmp(which, "BE", 2))
2271 if (*bmat != 'I' && *bmat != 'G')
2277 if (*lworkl < i__1 * i__1 + (*ncv << 3))
2281 if (iwork[11] < 1 || iwork[11] > 5)
2285 else if (iwork[11] == 1 && *bmat == 'G')
2289 else if (iwork[5] < 0 || iwork[5] > 1)
2293 else if (*nev == 1 && !strncmp(which, "BE", 2))
2311 *tol = GMX_DOUBLE_EPS;
2314 iwork[15] = *ncv - *nev;
2317 i__1 = i__2 * i__2 + (*ncv << 3);
2318 for (j = 1; j <= i__1; ++j)
2326 iwork[16] = iwork[3] + (iwork[8] << 1);
2327 iwork[1] = iwork[16] + *ncv;
2328 iwork[4] = iwork[1] + *ncv;
2330 iwork[7] = iwork[4] + i__1 * i__1;
2331 iwork[14] = iwork[7] + *ncv * 3;
2333 ipntr[4] = iwork[14];
2334 ipntr[5] = iwork[3];
2335 ipntr[6] = iwork[16];
2336 ipntr[7] = iwork[1];
2337 ipntr[11] = iwork[7];
2340 F77_FUNC(dsaup2, DSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
2341 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
2342 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
2343 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
2348 iparam[8] = iwork[15];
2355 iparam[3] = iwork[10];
2356 iparam[5] = iwork[15];
2376 F77_FUNC(dseupd, DSEUPD) (int * rvec,
2377 const char * howmny,
2402 int v_dim1, v_offset, z_dim1, z_offset, i__1;
2403 double d__1, d__2, d__3;
2405 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2417 double thres1 = 0, thres2 = 0;
2421 int leftptr, rghtptr;
2427 z_offset = 1 + z_dim1;
2432 v_offset = 1 + v_dim1;
2460 if (*ncv <= *nev || *ncv > *n)
2464 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
2465 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
2466 strncmp(which, "BE", 2))
2470 if (*bmat != 'I' && *bmat != 'G')
2474 if (*howmny != 'A' && *howmny != 'P' &&
2475 *howmny != 'S' && *rvec)
2479 if (*rvec && *howmny == 'S')
2484 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
2489 if (mode == 1 || mode == 2)
2491 strncpy(type__, "REGULR", 6);
2495 strncpy(type__, "SHIFTI", 6);
2499 strncpy(type__, "BUCKLE", 6);
2503 strncpy(type__, "CAYLEY", 6);
2509 if (mode == 1 && *bmat == 'G')
2513 if (*nev == 1 && !strncmp(which, "BE", 2))
2532 iw = iq + ldh * *ncv;
2533 next = iw + (*ncv << 1);
2539 irz = ipntr[11] + *ncv;
2543 eps23 = GMX_DOUBLE_EPS;
2544 eps23 = pow(eps23, c_b21);
2551 else if (*bmat == 'G')
2553 bnorm2 = F77_FUNC(dnrm2, DNRM2) (n, &workd[1], &c__1);
2559 if (!strncmp(which, "LM", 2) || !strncmp(which, "SM", 2) ||
2560 !strncmp(which, "LA", 2) || !strncmp(which, "SA", 2))
2564 else if (!strncmp(which, "BE", 2))
2568 ism = (*nev > nconv) ? *nev : nconv;
2571 thres1 = workl[ism];
2572 thres2 = workl[ilg];
2580 for (j = 0; j <= i__1; ++j)
2583 if (!strncmp(which, "LM", 2))
2585 if (fabs(workl[irz + j]) >= fabs(thres1))
2588 d__3 = fabs(workl[irz + j]);
2589 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2590 if (workl[ibd + j] <= *tol * tempbnd)
2596 else if (!strncmp(which, "SM", 2))
2598 if (fabs(workl[irz + j]) <= fabs(thres1))
2601 d__3 = fabs(workl[irz + j]);
2602 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2603 if (workl[ibd + j] <= *tol * tempbnd)
2609 else if (!strncmp(which, "LA", 2))
2611 if (workl[irz + j] >= thres1)
2614 d__3 = fabs(workl[irz + j]);
2615 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2616 if (workl[ibd + j] <= *tol * tempbnd)
2622 else if (!strncmp(which, "SA", 2))
2624 if (workl[irz + j] <= thres1)
2627 d__3 = fabs(workl[irz + j]);
2628 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2629 if (workl[ibd + j] <= *tol * tempbnd)
2635 else if (!strncmp(which, "BE", 2))
2637 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
2640 d__3 = fabs(workl[irz + j]);
2641 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2642 if (workl[ibd + j] <= *tol * tempbnd)
2650 reord = select[j + 1] || reord;
2659 F77_FUNC(dcopy, DCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2660 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2662 F77_FUNC(dsteqr, DSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
2684 if (select[leftptr])
2690 else if (!select[rghtptr])
2699 temp = workl[ihd + leftptr - 1];
2700 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2701 workl[ihd + rghtptr - 1] = temp;
2702 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
2704 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
2705 iq + *ncv * (leftptr - 1)], &c__1);
2706 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
2713 if (leftptr < rghtptr)
2722 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2728 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2729 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2732 if (!strncmp(type__, "REGULR", 6))
2737 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2741 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2748 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2749 if (!strncmp(type__, "SHIFTI", 6))
2752 for (k = 1; k <= i__1; ++k)
2754 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2757 else if (!strncmp(type__, "BUCKLE", 6))
2760 for (k = 1; k <= i__1; ++k)
2762 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
2766 else if (!strncmp(type__, "CAYLEY", 6))
2769 for (k = 1; k <= i__1; ++k)
2771 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
2772 workl[ihd + k - 1] - 1.);
2776 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2777 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2780 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2784 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2785 d__1 = bnorm2 / rnorm;
2786 F77_FUNC(dscal, DSCAL) (ncv, &d__1, &workl[ihb], &c__1);
2787 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2792 if (*rvec && *howmny == 'A')
2795 F77_FUNC(dgeqr2, DGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2798 F77_FUNC(dorm2r, DORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
2799 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
2800 F77_FUNC(dlacpy, DLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2803 for (j = 1; j <= i__1; ++j)
2805 workl[ihb + j - 1] = 0.;
2807 workl[ihb + *ncv - 1] = 1.;
2808 F77_FUNC(dorm2r, DORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
2809 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
2812 else if (*rvec && *howmny == 'S')
2817 if (!strncmp(type__, "REGULR", 6) && *rvec)
2821 for (j = 1; j <= i__1; ++j)
2823 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
2827 else if (strncmp(type__, "REGULR", 6) && *rvec)
2830 F77_FUNC(dscal, DSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
2831 if (!strncmp(type__, "SHIFTI", 6))
2835 for (k = 1; k <= i__1; ++k)
2837 d__2 = workl[iw + k - 1];
2838 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
2842 else if (!strncmp(type__, "BUCKLE", 6))
2846 for (k = 1; k <= i__1; ++k)
2848 d__2 = workl[iw + k - 1] - 1.;
2849 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
2853 else if (!strncmp(type__, "CAYLEY", 6))
2857 for (k = 1; k <= i__1; ++k)
2859 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2867 if (*rvec && (!strncmp(type__, "SHIFTI", 6) || !strncmp(type__, "CAYLEY", 6)))
2871 for (k = 0; k <= i__1; ++k)
2873 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2877 else if (*rvec && !strncmp(type__, "BUCKLE", 6))
2881 for (k = 0; k <= i__1; ++k)
2883 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
2889 if (strncmp(type__, "REGULR", 6))
2891 F77_FUNC(dger, DGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
2905 /* Selected single precision arpack routines */
2909 F77_FUNC(sstqrb, SSTQRB) (int * n,
2923 int i__, j, k, l, m;
2925 int l1, ii, mm, lm1, mm1, nm1;
2926 float rt1, rt2, eps;
2929 int lend, jtot, lendm1, lendp1, iscale;
2931 int lendsv, nmaxit, icompz;
2932 float ssfmax, ssfmin, safmin, minval, safmax, anorm;
2958 eps = GMX_FLOAT_EPS;
2962 minval = GMX_FLOAT_MIN;
2963 safmin = minval / GMX_FLOAT_EPS;
2964 safmax = 1. / safmin;
2965 ssfmax = sqrt(safmax) / 3.;
2966 ssfmin = sqrt(safmin) / eps2;
2971 for (j = 1; j <= i__1; ++j)
2997 for (m = l1; m <= i__1; ++m)
3004 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps)
3024 i__1 = lend - l + 1;
3025 anorm = F77_FUNC(slanst, SLANST) ("i", &i__1, &d__[l], &e[l]);
3034 i__1 = lend - l + 1;
3035 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
3038 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
3041 else if (anorm < ssfmin)
3044 i__1 = lend - l + 1;
3045 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
3048 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
3052 if (fabs(d__[lend]) < fabs(d__[l]))
3066 for (m = l; m <= i__1; ++m)
3070 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin)
3094 F77_FUNC(slaev2, SLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
3096 work[*n - 1 + l] = s;
3099 z__[l + 1] = c__ * tst - s * z__[l];
3100 z__[l] = s * tst + c__ * z__[l];
3104 F77_FUNC(slae2, SLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
3123 g = (d__[l + 1] - p) / (e[l] * 2.);
3124 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3125 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
3133 for (i__ = mm1; i__ >= i__1; --i__)
3137 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3142 g = d__[i__ + 1] - p;
3143 r__ = (d__[i__] - g) * s + c__ * 2. * b;
3145 d__[i__ + 1] = g + p;
3151 work[*n - 1 + i__] = -s;
3160 F77_FUNC(slasr, SLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
3187 for (m = l; m >= i__1; --m)
3189 d__2 = fabs(e[m - 1]);
3191 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin)
3215 F77_FUNC(slaev2, SLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
3219 z__[l] = c__ * tst - s * z__[l - 1];
3220 z__[l - 1] = s * tst + c__ * z__[l - 1];
3224 F77_FUNC(slae2, SLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
3244 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
3245 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3246 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
3254 for (i__ = m; i__ <= i__1; ++i__)
3258 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3264 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
3272 work[*n - 1 + i__] = s;
3281 F77_FUNC(slasr, SLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
3304 i__1 = lendsv - lsv + 1;
3305 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
3307 i__1 = lendsv - lsv;
3308 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
3311 else if (iscale == 2)
3313 i__1 = lendsv - lsv + 1;
3314 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
3316 i__1 = lendsv - lsv;
3317 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
3326 for (i__ = 1; i__ <= i__1; ++i__)
3339 F77_FUNC(slasrt, SLASRT) ("i", n, &d__[1], info);
3346 for (ii = 2; ii <= i__1; ++ii)
3352 for (j = ii; j <= i__2; ++j)
3378 F77_FUNC(sgetv0, SGETV0) (int * ido,
3397 int v_dim1, v_offset, i__1;
3405 v_offset = 1 + v_dim1;
3421 F77_FUNC(slarnv, SLARNV) (&idist, &iwork[1], n, &resid[1]);
3428 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3447 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
3453 else if (*bmat == 'I')
3455 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3464 workd[*n * 3 + 4] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3465 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
3467 else if (*bmat == 'I')
3469 workd[*n * 3 + 4] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3471 *rnorm = workd[*n * 3 + 4];
3481 F77_FUNC(sgemv, SGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
3482 &workd[*n + 1], &c__1);
3484 F77_FUNC(sgemv, SGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
3485 c_b22, &resid[1], &c__1);
3489 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3495 else if (*bmat == 'I')
3497 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3504 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3505 *rnorm = sqrt(fabs(*rnorm));
3507 else if (*bmat == 'I')
3509 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3512 if (*rnorm > workd[*n * 3 + 4] * .717f)
3521 workd[*n * 3 + 4] = *rnorm;
3528 for (jj = 1; jj <= i__1; ++jj)
3549 F77_FUNC(ssapps, SSAPPS) (int * n,
3566 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
3570 float r__, s, a1, a2, a3, a4;
3581 v_offset = 1 + v_dim1;
3584 h_offset = 1 + h_dim1;
3587 q_offset = 1 + q_dim1;
3590 epsmch = GMX_FLOAT_EPS;
3594 kplusp = *kev + *np;
3596 F77_FUNC(slaset, SLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3604 for (jj = 1; jj <= i__1; ++jj)
3612 for (i__ = istart; i__ <= i__2; ++i__)
3614 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
3615 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3617 h__[i__ + 1 + h_dim1] = 0.;
3628 f = h__[istart + (h_dim1 << 1)] - shift[jj];
3629 g = h__[istart + 1 + h_dim1];
3630 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3632 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
3634 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
3636 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
3638 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
3640 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
3641 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3642 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
3645 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
3646 for (j = 1; j <= i__2; ++j)
3648 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
3650 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
3651 c__ * q[j + (istart + 1) * q_dim1];
3652 q[j + istart * q_dim1] = a1;
3657 for (i__ = istart + 1; i__ <= i__2; ++i__)
3660 f = h__[i__ + h_dim1];
3661 g = s * h__[i__ + 1 + h_dim1];
3663 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3664 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3673 h__[i__ + h_dim1] = r__;
3675 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
3677 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
3679 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
3681 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
3684 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
3685 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3686 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
3689 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
3690 for (j = 1; j <= i__3; ++j)
3692 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
3694 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
3695 c__ * q[j + (i__ + 1) * q_dim1];
3696 q[j + i__ * q_dim1] = a1;
3705 if (h__[iend + h_dim1] < 0.)
3707 h__[iend + h_dim1] = -h__[iend + h_dim1];
3708 F77_FUNC(sscal, SSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3717 for (i__ = itop; i__ <= i__2; ++i__)
3719 if (h__[i__ + 1 + h_dim1] > 0.)
3731 for (i__ = itop; i__ <= i__1; ++i__)
3733 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
3734 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3736 h__[i__ + 1 + h_dim1] = 0.;
3741 if (h__[*kev + 1 + h_dim1] > 0.)
3743 F77_FUNC(sgemv, SGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
3744 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
3748 for (i__ = 1; i__ <= i__1; ++i__)
3750 i__2 = kplusp - i__ + 1;
3751 F77_FUNC(sgemv, SGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
3752 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
3753 F77_FUNC(scopy, SCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
3758 F77_FUNC(slacpy, SLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3760 if (h__[*kev + 1 + h_dim1] > 0.)
3762 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3765 F77_FUNC(sscal, SSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3766 if (h__[*kev + 1 + h_dim1] > 0.)
3768 F77_FUNC(saxpy, SAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
3782 F77_FUNC(ssortr, SSORTR) (const char * which,
3797 if (!strncmp(which, "SA", 2))
3806 for (i__ = igap; i__ <= i__1; ++i__)
3816 if (x1[j] < x1[j + igap])
3819 x1[j] = x1[j + igap];
3820 x1[j + igap] = temp;
3824 x2[j] = x2[j + igap];
3825 x2[j + igap] = temp;
3841 else if (!strncmp(which, "SM", 2))
3850 for (i__ = igap; i__ <= i__1; ++i__)
3860 if (fabs(x1[j]) < fabs(x1[j + igap]))
3863 x1[j] = x1[j + igap];
3864 x1[j + igap] = temp;
3868 x2[j] = x2[j + igap];
3869 x2[j + igap] = temp;
3885 else if (!strncmp(which, "LA", 2))
3894 for (i__ = igap; i__ <= i__1; ++i__)
3904 if (x1[j] > x1[j + igap])
3907 x1[j] = x1[j + igap];
3908 x1[j + igap] = temp;
3912 x2[j] = x2[j + igap];
3913 x2[j + igap] = temp;
3929 else if (!strncmp(which, "LM", 2))
3939 for (i__ = igap; i__ <= i__1; ++i__)
3949 if (fabs(x1[j]) > fabs(x1[j + igap]))
3952 x1[j] = x1[j + igap];
3953 x1[j + igap] = temp;
3957 x2[j] = x2[j + igap];
3958 x2[j + igap] = temp;
3983 F77_FUNC(ssesrt, SSESRT) (const char * which,
3991 int a_dim1, a_offset, i__1;
3998 a_offset = 1 + a_dim1 * 0;
4003 if (!strncmp(which, "SA", 2))
4012 for (i__ = igap; i__ <= i__1; ++i__)
4022 if (x[j] < x[j + igap])
4029 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4030 a_dim1 + 1], &c__1);
4046 else if (!strncmp(which, "SM", 2))
4055 for (i__ = igap; i__ <= i__1; ++i__)
4065 if (fabs(x[j]) < fabs(x[j + igap]))
4072 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4073 a_dim1 + 1], &c__1);
4089 else if (!strncmp(which, "LA", 2))
4098 for (i__ = igap; i__ <= i__1; ++i__)
4108 if (x[j] > x[j + igap])
4115 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4116 a_dim1 + 1], &c__1);
4132 else if (!strncmp(which, "LM", 2))
4141 for (i__ = igap; i__ <= i__1; ++i__)
4151 if (fabs(x[j]) > fabs(x[j + igap]))
4158 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4159 a_dim1 + 1], &c__1);
4184 F77_FUNC(ssgets, SSGETS) (int * ishift,
4200 if (!strncmp(which, "BE", 2))
4203 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
4207 i__1 = (kevd2 < *np) ? kevd2 : *np;
4208 i__2 = (kevd2 > *np) ? kevd2 : *np;
4209 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[1], &c__1,
4210 &ritz[i__2 + 1], &c__1);
4211 i__1 = (kevd2 < *np) ? kevd2 : *np;
4212 i__2 = (kevd2 > *np) ? kevd2 : *np;
4213 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[1], &c__1,
4214 &bounds[i__2 + 1], &c__1);
4221 F77_FUNC(ssortr, SSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
4224 if (*ishift == 1 && *np > 0)
4227 F77_FUNC(ssortr, SSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
4228 F77_FUNC(scopy, SCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
4238 F77_FUNC(ssconv, SSCONV) (int * n,
4254 eps23 = GMX_FLOAT_EPS;
4255 eps23 = pow(eps23, c_b3);
4259 for (i__ = 1; i__ <= i__1; ++i__)
4263 d__3 = fabs(ritz[i__]);
4264 temp = (d__2 > d__3) ? d__2 : d__3;
4265 if (bounds[i__] <= *tol * temp)
4276 F77_FUNC(sseigt, SSEIGT) (float * rnorm,
4286 int h_dim1, h_offset, i__1;
4295 h_offset = 1 + h_dim1;
4298 F77_FUNC(scopy, SCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
4300 F77_FUNC(scopy, SCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
4301 F77_FUNC(sstqrb, SSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
4308 for (k = 1; k <= i__1; ++k)
4310 bounds[k] = *rnorm * fabs(bounds[k]);
4323 F77_FUNC(ssaitr, SSAITR) (int * ido,
4347 int h_dim1, h_offset, v_dim1, v_offset, i__1;
4351 float safmin, minval;
4357 v_offset = 1 + v_dim1;
4360 h_offset = 1 + h_dim1;
4364 minval = GMX_FLOAT_MIN;
4365 safmin = minval / GMX_FLOAT_EPS;
4379 iwork[9] = iwork[8] + *n;
4380 iwork[10] = iwork[9] + *n;
4418 F77_FUNC(sgetv0, sgetv0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
4419 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
4432 *info = iwork[12] - 1;
4439 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
4440 if (*rnorm >= safmin)
4442 temp1 = 1. / *rnorm;
4443 F77_FUNC(sscal, SSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
4444 F77_FUNC(sscal, SSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
4449 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
4450 v_dim1 + 1], n, &infol);
4451 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
4456 F77_FUNC(scopy, SCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
4457 ipntr[1] = iwork[10];
4458 ipntr[2] = iwork[9];
4459 ipntr[3] = iwork[8];
4468 F77_FUNC(scopy, SCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
4477 ipntr[1] = iwork[9];
4478 ipntr[2] = iwork[8];
4483 else if (*bmat == 'I')
4485 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4495 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
4497 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
4499 else if (*bmat == 'G')
4501 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4503 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
4505 else if (*bmat == 'I')
4507 workd[*n * 3 + 3] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4512 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
4513 &c__1, &c_b42, &workd[iwork[9]], &c__1);
4515 else if (*mode == 2)
4517 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
4518 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
4521 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4522 c__1, &c_b18, &resid[1], &c__1);
4524 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
4525 if (iwork[12] == 1 || iwork[4] == 1)
4527 h__[iwork[12] + h_dim1] = 0.;
4531 h__[iwork[12] + h_dim1] = *rnorm;
4539 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4540 ipntr[1] = iwork[9];
4541 ipntr[2] = iwork[8];
4546 else if (*bmat == 'I')
4548 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4556 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4557 *rnorm = sqrt(fabs(*rnorm));
4559 else if (*bmat == 'I')
4561 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4564 if (*rnorm > workd[*n * 3 + 3] * .717f)
4571 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
4572 c__1, &c_b42, &workd[iwork[9]], &c__1);
4574 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4575 c__1, &c_b18, &resid[1], &c__1);
4577 if (iwork[12] == 1 || iwork[4] == 1)
4579 h__[iwork[12] + h_dim1] = 0.;
4581 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
4586 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4587 ipntr[1] = iwork[9];
4588 ipntr[2] = iwork[8];
4593 else if (*bmat == 'I')
4595 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4602 workd[*n * 3 + 2] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4604 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
4606 else if (*bmat == 'I')
4608 workd[*n * 3 + 2] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4612 if (workd[*n * 3 + 2] > *rnorm * .717f)
4615 *rnorm = workd[*n * 3 + 2];
4621 *rnorm = workd[*n * 3 + 2];
4629 for (jj = 1; jj <= i__1; ++jj)
4641 if (h__[iwork[12] + h_dim1] < 0.)
4643 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4644 if (iwork[12] < *k + *np)
4646 F77_FUNC(sscal, SSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4650 F77_FUNC(sscal, SSCAL) (n, &c_b50, &resid[1], &c__1);
4655 if (iwork[12] > *k + *np)
4675 F77_FUNC(ssaup2, SSAUP2) (int * ido,
4705 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
4724 v_offset = 1 + v_dim1;
4727 h_offset = 1 + h_dim1;
4730 q_offset = 1 + q_dim1;
4734 eps23 = GMX_FLOAT_EPS;
4735 eps23 = pow(eps23, c_b3);
4748 iwork[7] = iwork[9] + iwork[10];
4771 F77_FUNC(sgetv0, SGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
4772 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
4780 if (workd[*n * 3 + 1] == 0.)
4805 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
4806 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
4832 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[
4833 v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[
4851 F77_FUNC(sseigt, SSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
4852 bounds[1], &workl[1], &ierr);
4860 F77_FUNC(scopy, SCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4861 F77_FUNC(scopy, SCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4865 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4867 F77_FUNC(scopy, SCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4868 F77_FUNC(ssconv, SSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4873 for (j = 1; j <= i__1; ++j)
4875 if (bounds[j] == 0.)
4882 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
4885 if (!strncmp(which, "BE", 2))
4888 strncpy(wprime, "SA", 2);
4889 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4891 nevm2 = *nev - nevd2;
4894 i__1 = (nevd2 < *np) ? nevd2 : *np;
4895 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4896 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
4897 &ritz[((i__2 > i__3) ? i__2 : i__3)],
4899 i__1 = (nevd2 < *np) ? nevd2 : *np;
4900 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4901 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
4902 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
4910 if (!strncmp(which, "LM", 2))
4912 strncpy(wprime, "SM", 2);
4914 if (!strncmp(which, "SM", 2))
4916 strncpy(wprime, "LM", 2);
4918 if (!strncmp(which, "LA", 2))
4920 strncpy(wprime, "SA", 2);
4922 if (!strncmp(which, "SA", 2))
4924 strncpy(wprime, "LA", 2);
4927 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4932 for (j = 1; j <= i__1; ++j)
4935 d__3 = fabs(ritz[j]);
4936 temp = (d__2 > d__3) ? d__2 : d__3;
4940 strncpy(wprime, "LA", 2);
4941 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4944 for (j = 1; j <= i__1; ++j)
4947 d__3 = fabs(ritz[j]);
4948 temp = (d__2 > d__3) ? d__2 : d__3;
4952 if (!strncmp(which, "BE", 2))
4955 strncpy(wprime, "LA", 2);
4956 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4961 F77_FUNC(ssortr, SSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4965 h__[h_dim1 + 1] = workd[*n * 3 + 1];
4968 if (iwork[6] > *mxiter && iwork[8] < *nev)
4972 if (*np == 0 && iwork[8] < iwork[9])
4981 else if (iwork[8] < *nev && *ishift == 1)
4984 i__1 = iwork[8], i__2 = *np / 2;
4985 *nev += (i__1 < i__2) ? i__1 : i__2;
4986 if (*nev == 1 && iwork[7] >= 6)
4988 *nev = iwork[7] / 2;
4990 else if (*nev == 1 && iwork[7] > 2)
4994 *np = iwork[7] - *nev;
4999 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
5019 F77_FUNC(scopy, SCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
5022 F77_FUNC(ssapps, SSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
5023 resid[1], &q[q_offset], ldq, &workd[1]);
5028 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
5035 else if (*bmat == 'I')
5037 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
5044 workd[*n * 3 + 1] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
5045 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
5047 else if (*bmat == 'I')
5049 workd[*n * 3 + 1] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
5071 F77_FUNC(ssaupd, SSAUPD) (int * ido,
5089 int v_dim1, v_offset, i__1, i__2;
5095 v_offset = 1 + v_dim1;
5107 iwork[5] = iparam[1];
5108 iwork[10] = iparam[3];
5109 iwork[12] = iparam[4];
5112 iwork[11] = iparam[7];
5123 else if (*ncv <= *nev || *ncv > *n)
5129 iwork[15] = *ncv - *nev;
5135 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
5136 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
5137 strncmp(which, "BE", 2))
5141 if (*bmat != 'I' && *bmat != 'G')
5147 if (*lworkl < i__1 * i__1 + (*ncv << 3))
5151 if (iwork[11] < 1 || iwork[11] > 5)
5155 else if (iwork[11] == 1 && *bmat == 'G')
5159 else if (iwork[5] < 0 || iwork[5] > 1)
5163 else if (*nev == 1 && !strncmp(which, "BE", 2))
5181 *tol = GMX_FLOAT_EPS;
5184 iwork[15] = *ncv - *nev;
5187 i__1 = i__2 * i__2 + (*ncv << 3);
5188 for (j = 1; j <= i__1; ++j)
5196 iwork[16] = iwork[3] + (iwork[8] << 1);
5197 iwork[1] = iwork[16] + *ncv;
5198 iwork[4] = iwork[1] + *ncv;
5200 iwork[7] = iwork[4] + i__1 * i__1;
5201 iwork[14] = iwork[7] + *ncv * 3;
5203 ipntr[4] = iwork[14];
5204 ipntr[5] = iwork[3];
5205 ipntr[6] = iwork[16];
5206 ipntr[7] = iwork[1];
5207 ipntr[11] = iwork[7];
5210 F77_FUNC(ssaup2, SSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
5211 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
5212 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
5213 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
5218 iparam[8] = iwork[15];
5225 iparam[3] = iwork[10];
5226 iparam[5] = iwork[15];
5246 F77_FUNC(sseupd, SSEUPD) (int * rvec,
5247 const char * howmny,
5272 int v_dim1, v_offset, z_dim1, z_offset, i__1;
5273 float d__1, d__2, d__3;
5275 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
5287 float thres1 = 0, thres2 = 0;
5291 int leftptr, rghtptr;
5297 z_offset = 1 + z_dim1;
5302 v_offset = 1 + v_dim1;
5330 if (*ncv <= *nev || *ncv > *n)
5334 if (strncmp(which, "LM", 2) && strncmp(which, "SM", 2) &&
5335 strncmp(which, "LA", 2) && strncmp(which, "SA", 2) &&
5336 strncmp(which, "BE", 2))
5340 if (*bmat != 'I' && *bmat != 'G')
5344 if (*howmny != 'A' && *howmny != 'P' &&
5345 *howmny != 'S' && *rvec)
5349 if (*rvec && *howmny == 'S')
5354 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
5359 if (mode == 1 || mode == 2)
5361 strncpy(type__, "REGULR", 6);
5365 strncpy(type__, "SHIFTI", 6);
5369 strncpy(type__, "BUCKLE", 6);
5373 strncpy(type__, "CAYLEY", 6);
5379 if (mode == 1 && *bmat == 'G')
5383 if (*nev == 1 && !strncmp(which, "BE", 2))
5402 iw = iq + ldh * *ncv;
5403 next = iw + (*ncv << 1);
5409 irz = ipntr[11] + *ncv;
5413 eps23 = GMX_FLOAT_EPS;
5414 eps23 = pow(eps23, c_b21);
5421 else if (*bmat == 'G')
5423 bnorm2 = F77_FUNC(snrm2, SNRM2) (n, &workd[1], &c__1);
5429 if (!strncmp(which, "LM", 2) || !strncmp(which, "SM", 2) ||
5430 !strncmp(which, "LA", 2) || !strncmp(which, "SA", 2))
5434 else if (!strncmp(which, "BE", 2))
5438 ism = (*nev > nconv) ? *nev : nconv;
5441 thres1 = workl[ism];
5442 thres2 = workl[ilg];
5450 for (j = 0; j <= i__1; ++j)
5453 if (!strncmp(which, "LM", 2))
5455 if (fabs(workl[irz + j]) >= fabs(thres1))
5458 d__3 = fabs(workl[irz + j]);
5459 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5460 if (workl[ibd + j] <= *tol * tempbnd)
5466 else if (!strncmp(which, "SM", 2))
5468 if (fabs(workl[irz + j]) <= fabs(thres1))
5471 d__3 = fabs(workl[irz + j]);
5472 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5473 if (workl[ibd + j] <= *tol * tempbnd)
5479 else if (!strncmp(which, "LA", 2))
5481 if (workl[irz + j] >= thres1)
5484 d__3 = fabs(workl[irz + j]);
5485 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5486 if (workl[ibd + j] <= *tol * tempbnd)
5492 else if (!strncmp(which, "SA", 2))
5494 if (workl[irz + j] <= thres1)
5497 d__3 = fabs(workl[irz + j]);
5498 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5499 if (workl[ibd + j] <= *tol * tempbnd)
5505 else if (!strncmp(which, "BE", 2))
5507 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
5510 d__3 = fabs(workl[irz + j]);
5511 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5512 if (workl[ibd + j] <= *tol * tempbnd)
5520 reord = select[j + 1] || reord;
5529 F77_FUNC(scopy, SCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
5530 F77_FUNC(scopy, SCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
5532 F77_FUNC(ssteqr, SSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
5554 if (select[leftptr])
5560 else if (!select[rghtptr])
5569 temp = workl[ihd + leftptr - 1];
5570 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
5571 workl[ihd + rghtptr - 1] = temp;
5572 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
5574 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
5575 iq + *ncv * (leftptr - 1)], &c__1);
5576 F77_FUNC(scopy, SCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
5583 if (leftptr < rghtptr)
5592 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5598 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
5599 F77_FUNC(scopy, SCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
5602 if (!strncmp(type__, "REGULR", 6))
5607 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5611 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5618 F77_FUNC(scopy, SCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
5619 if (!strncmp(type__, "SHIFTI", 6))
5622 for (k = 1; k <= i__1; ++k)
5624 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
5627 else if (!strncmp(type__, "BUCKLE", 6))
5630 for (k = 1; k <= i__1; ++k)
5632 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
5636 else if (!strncmp(type__, "CAYLEY", 6))
5639 for (k = 1; k <= i__1; ++k)
5641 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
5642 workl[ihd + k - 1] - 1.);
5646 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5647 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
5650 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5654 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5655 d__1 = bnorm2 / rnorm;
5656 F77_FUNC(sscal, SSCAL) (ncv, &d__1, &workl[ihb], &c__1);
5657 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
5662 if (*rvec && *howmny == 'A')
5665 F77_FUNC(sgeqr2, SGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
5668 F77_FUNC(sorm2r, SORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
5669 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
5670 F77_FUNC(slacpy, SLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
5673 for (j = 1; j <= i__1; ++j)
5675 workl[ihb + j - 1] = 0.;
5677 workl[ihb + *ncv - 1] = 1.;
5678 F77_FUNC(sorm2r, SORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
5679 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
5682 else if (*rvec && *howmny == 'S')
5687 if (!strncmp(type__, "REGULR", 6) && *rvec)
5691 for (j = 1; j <= i__1; ++j)
5693 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
5697 else if (strncmp(type__, "REGULR", 6) && *rvec)
5700 F77_FUNC(sscal, SSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
5701 if (!strncmp(type__, "SHIFTI", 6))
5705 for (k = 1; k <= i__1; ++k)
5707 d__2 = workl[iw + k - 1];
5708 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
5712 else if (!strncmp(type__, "BUCKLE", 6))
5716 for (k = 1; k <= i__1; ++k)
5718 d__2 = workl[iw + k - 1] - 1.;
5719 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
5723 else if (!strncmp(type__, "CAYLEY", 6))
5727 for (k = 1; k <= i__1; ++k)
5729 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
5737 if (*rvec && (!strncmp(type__, "SHIFTI", 6) || !strncmp(type__, "CAYLEY", 6)))
5741 for (k = 0; k <= i__1; ++k)
5743 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
5747 else if (*rvec && !strncmp(type__, "BUCKLE", 6))
5751 for (k = 0; k <= i__1; ++k)
5753 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
5759 if (strncmp(type__, "REGULR", 6))
5761 F77_FUNC(sger, SGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[