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"
37 #include "gmx_arpack.h"
39 #include "gmx_lapack.h"
42 F77_FUNC(dstqrb,DSTQRB)(int * n,
58 int l1, ii, mm, lm1, mm1, nm1;
62 int lend, jtot, lendm1, lendp1, iscale;
64 int lendsv, nmaxit, icompz;
65 double ssfmax, ssfmin,safmin,minval,safmax,anorm;
92 minval = GMX_DOUBLE_MIN;
93 safmin = minval / GMX_DOUBLE_EPS;
95 ssfmax = sqrt(safmax) / 3.;
96 ssfmin = sqrt(safmin) / eps2;
100 for (j = 1; j <= i__1; ++j) {
122 for (m = l1; m <= i__1; ++m) {
127 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps) {
146 anorm =F77_FUNC(dlanst,DLANST)("i", &i__1, &d__[l], &e[l]);
151 if (anorm > ssfmax) {
154 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
157 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
159 } else if (anorm < ssfmin) {
162 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
165 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
169 if (fabs(d__[lend]) < fabs(d__[l])) {
180 for (m = l; m <= i__1; ++m) {
183 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin) {
202 F77_FUNC(dlaev2,DLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
204 work[*n - 1 + l] = s;
207 z__[l + 1] = c__ * tst - s * z__[l];
208 z__[l] = s * tst + c__ * z__[l];
210 F77_FUNC(dlae2,DLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
222 if (jtot == nmaxit) {
227 g = (d__[l + 1] - p) / (e[l] * 2.);
228 r__ =F77_FUNC(dlapy2,DLAPY2)(&g, &c_b31);
229 g = d__[m] - p + e[l] / (g + ((g>0) ? r__ : -r__ ));
237 for (i__ = mm1; i__ >= i__1; --i__) {
240 F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
244 g = d__[i__ + 1] - p;
245 r__ = (d__[i__] - g) * s + c__ * 2. * b;
247 d__[i__ + 1] = g + p;
252 work[*n - 1 + i__] = -s;
260 F77_FUNC(dlasr,DLASR)("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
283 for (m = l; m >= i__1; --m) {
284 d__2 = fabs(e[m - 1]);
286 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin) {
305 F77_FUNC(dlaev2,DLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
309 z__[l] = c__ * tst - s * z__[l - 1];
310 z__[l - 1] = s * tst + c__ * z__[l - 1];
312 F77_FUNC(dlae2,DLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
324 if (jtot == nmaxit) {
330 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
331 r__ =F77_FUNC(dlapy2,DLAPY2)(&g, &c_b31);
332 g = d__[m] - p + e[l - 1] / (g + ((g>0) ? r__ : -r__ ));
340 for (i__ = m; i__ <= i__1; ++i__) {
343 F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
348 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
355 work[*n - 1 + i__] = s;
363 F77_FUNC(dlasr,DLASR)("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
384 i__1 = lendsv - lsv + 1;
385 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
388 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
390 } else if (iscale == 2) {
391 i__1 = lendsv - lsv + 1;
392 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
395 F77_FUNC(dlascl,DLASCL)("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
403 for (i__ = 1; i__ <= i__1; ++i__) {
413 F77_FUNC(dlasrt,DLASRT)("i", n, &d__[1], info);
418 for (ii = 2; ii <= i__1; ++ii) {
423 for (j = ii; j <= i__2; ++j) {
446 F77_FUNC(dgetv0,DGETV0)(int * ido,
465 int v_dim1, v_offset, i__1;
473 v_offset = 1 + v_dim1;
487 F77_FUNC(dlarnv,DLARNV)(&idist, &iwork[1], n, &resid[1]);
493 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
509 F77_FUNC(dcopy,DCOPY)(n, &workd[*n + 1], &c__1, &resid[1], &c__1);
514 } else if (*bmat == 'I') {
515 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
523 workd[*n * 3 + 4] =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
524 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
525 } else if (*bmat == 'I') {
526 workd[*n * 3 + 4] =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
528 *rnorm = workd[*n * 3 + 4];
537 F77_FUNC(dgemv,DGEMV)("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
538 &workd[*n + 1], &c__1);
540 F77_FUNC(dgemv,DGEMV)("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
541 c_b22, &resid[1], &c__1);
544 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
549 } else if (*bmat == 'I') {
550 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
556 *rnorm =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
557 *rnorm = sqrt(fabs(*rnorm));
558 } else if (*bmat == 'I') {
559 *rnorm =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
562 if (*rnorm > workd[*n * 3 + 4] * .717f) {
569 workd[*n * 3 + 4] = *rnorm;
574 for (jj = 1; jj <= i__1; ++jj) {
594 F77_FUNC(dsapps,DSAPPS)(int * n,
611 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
615 double r__, s, a1, a2, a3, a4;
626 v_offset = 1 + v_dim1;
629 h_offset = 1 + h_dim1;
632 q_offset = 1 + q_dim1;
635 epsmch = GMX_DOUBLE_EPS;
641 F77_FUNC(dlaset,DLASET)("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
648 for (jj = 1; jj <= i__1; ++jj) {
655 for (i__ = istart; i__ <= i__2; ++i__) {
656 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
657 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
658 h__[i__ + 1 + h_dim1] = 0.;
668 f = h__[istart + (h_dim1 << 1)] - shift[jj];
669 g = h__[istart + 1 + h_dim1];
670 F77_FUNC(dlartg,DLARTG)(&f, &g, &c__, &s, &r__);
672 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
674 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
676 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
678 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
680 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
681 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
682 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
685 i__2 = (i__3<kplusp) ? i__3 : kplusp;
686 for (j = 1; j <= i__2; ++j) {
687 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
689 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
690 c__ * q[j + (istart + 1) * q_dim1];
691 q[j + istart * q_dim1] = a1;
696 for (i__ = istart + 1; i__ <= i__2; ++i__) {
698 f = h__[i__ + h_dim1];
699 g = s * h__[i__ + 1 + h_dim1];
701 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
702 F77_FUNC(dlartg,DLARTG)(&f, &g, &c__, &s, &r__);
710 h__[i__ + h_dim1] = r__;
712 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
714 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
716 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
718 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
721 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
722 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
723 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
726 i__3 = (i__4<kplusp) ? i__4 : kplusp;
727 for (j = 1; j <= i__3; ++j) {
728 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
730 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
731 c__ * q[j + (i__ + 1) * q_dim1];
732 q[j + i__ * q_dim1] = a1;
741 if (h__[iend + h_dim1] < 0.) {
742 h__[iend + h_dim1] = -h__[iend + h_dim1];
743 F77_FUNC(dscal,DSCAL)(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
751 for (i__ = itop; i__ <= i__2; ++i__) {
752 if (h__[i__ + 1 + h_dim1] > 0.) {
763 for (i__ = itop; i__ <= i__1; ++i__) {
764 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
765 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
766 h__[i__ + 1 + h_dim1] = 0.;
771 if (h__[*kev + 1 + h_dim1] > 0.) {
772 F77_FUNC(dgemv,DGEMV)("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
773 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
777 for (i__ = 1; i__ <= i__1; ++i__) {
778 i__2 = kplusp - i__ + 1;
779 F77_FUNC(dgemv,DGEMV)("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
780 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
781 F77_FUNC(dcopy,DCOPY)(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
786 F77_FUNC(dlacpy,DLACPY)("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
788 if (h__[*kev + 1 + h_dim1] > 0.) {
789 F77_FUNC(dcopy,DCOPY)(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
792 F77_FUNC(dscal,DSCAL)(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
793 if (h__[*kev + 1 + h_dim1] > 0.) {
794 F77_FUNC(daxpy,DAXPY)(n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
808 F77_FUNC(dsortr,DSORTR)(const char * which,
823 if (!strncmp(which, "SA", 2)) {
830 for (i__ = igap; i__ <= i__1; ++i__) {
838 if (x1[j] < x1[j + igap]) {
840 x1[j] = x1[j + igap];
844 x2[j] = x2[j + igap];
858 } else if (!strncmp(which, "SM", 2)) {
865 for (i__ = igap; i__ <= i__1; ++i__) {
873 if (fabs(x1[j]) < fabs(x1[j + igap])) {
875 x1[j] = x1[j + igap];
879 x2[j] = x2[j + igap];
893 } else if (!strncmp(which, "LA", 2)) {
900 for (i__ = igap; i__ <= i__1; ++i__) {
908 if (x1[j] > x1[j + igap]) {
910 x1[j] = x1[j + igap];
914 x2[j] = x2[j + igap];
928 } else if (!strncmp(which, "LM", 2)) {
936 for (i__ = igap; i__ <= i__1; ++i__) {
944 if (fabs(x1[j]) > fabs(x1[j + igap])) {
946 x1[j] = x1[j + igap];
950 x2[j] = x2[j + igap];
974 F77_FUNC(dsesrt,DSESRT)(const char * which,
982 int a_dim1, a_offset, i__1;
989 a_offset = 1 + a_dim1 * 0;
994 if (!strncmp(which, "SA", 2)) {
1001 for (i__ = igap; i__ <= i__1; ++i__) {
1009 if (x[j] < x[j + igap]) {
1014 F77_FUNC(dswap,DSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1015 a_dim1 + 1], &c__1);
1028 } else if (!strncmp(which, "SM", 2)) {
1035 for (i__ = igap; i__ <= i__1; ++i__) {
1043 if (fabs(x[j]) < fabs(x[j + igap])) {
1048 F77_FUNC(dswap,DSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1049 a_dim1 + 1], &c__1);
1062 } else if (!strncmp(which, "LA", 2)) {
1069 for (i__ = igap; i__ <= i__1; ++i__) {
1077 if (x[j] > x[j + igap]) {
1082 F77_FUNC(dswap,DSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1083 a_dim1 + 1], &c__1);
1096 } else if (!strncmp(which, "LM", 2)) {
1103 for (i__ = igap; i__ <= i__1; ++i__) {
1111 if (fabs(x[j]) > fabs(x[j + igap])) {
1116 F77_FUNC(dswap,DSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1117 a_dim1 + 1], &c__1);
1140 F77_FUNC(dsgets,DSGETS)(int * ishift,
1156 if (!strncmp(which, "BE", 2)) {
1158 F77_FUNC(dsortr,DSORTR)("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1161 i__1 = (kevd2<*np) ? kevd2 : *np;
1162 i__2 = (kevd2>*np) ? kevd2 : *np;
1163 F77_FUNC(dswap,DSWAP)(&i__1, &ritz[1], &c__1,
1164 &ritz[i__2 + 1], &c__1);
1165 i__1 = (kevd2<*np) ? kevd2 : *np;
1166 i__2 = (kevd2>*np) ? kevd2 : *np;
1167 F77_FUNC(dswap,DSWAP)(&i__1, &bounds[1], &c__1,
1168 &bounds[i__2 + 1], &c__1);
1173 F77_FUNC(dsortr,DSORTR)(which, &c__1, &i__1, &ritz[1], &bounds[1]);
1176 if (*ishift == 1 && *np > 0) {
1178 F77_FUNC(dsortr,DSORTR)("SM", &c__1, np, &bounds[1], &ritz[1]);
1179 F77_FUNC(dcopy,DCOPY)(np, &ritz[1], &c__1, &shifts[1], &c__1);
1189 F77_FUNC(dsconv,DSCONV)(int * n,
1205 eps23 = GMX_DOUBLE_EPS;
1206 eps23 = pow(eps23, c_b3);
1210 for (i__ = 1; i__ <= i__1; ++i__) {
1213 d__3 = fabs(ritz[i__]);
1214 temp = (d__2 > d__3) ? d__2 : d__3;
1215 if (bounds[i__] <= *tol * temp) {
1225 F77_FUNC(dseigt,DSEIGT)(double * rnorm,
1235 int h_dim1, h_offset, i__1;
1244 h_offset = 1 + h_dim1;
1247 F77_FUNC(dcopy,DCOPY)(n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1249 F77_FUNC(dcopy,DCOPY)(&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1250 F77_FUNC(dstqrb,DSTQRB)(n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1256 for (k = 1; k <= i__1; ++k) {
1257 bounds[k] = *rnorm * fabs(bounds[k]);
1270 F77_FUNC(dsaitr,DSAITR)(int * ido,
1294 int h_dim1, h_offset, v_dim1, v_offset, i__1;
1298 double safmin,minval;
1304 v_offset = 1 + v_dim1;
1307 h_offset = 1 + h_dim1;
1311 minval = GMX_DOUBLE_MIN;
1312 safmin = minval / GMX_DOUBLE_EPS;
1325 iwork[9] = iwork[8] + *n;
1326 iwork[10] = iwork[9] + *n;
1329 if (iwork[5] == 1) {
1332 if (iwork[6] == 1) {
1335 if (iwork[2] == 1) {
1338 if (iwork[3] == 1) {
1341 if (iwork[4] == 1) {
1358 F77_FUNC(dgetv0,DGETV0)(ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
1359 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
1365 if (iwork[11] <= 3) {
1369 *info = iwork[12] - 1;
1376 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1377 if (*rnorm >= safmin) {
1378 temp1 = 1. / *rnorm;
1379 F77_FUNC(dscal,DSCAL)(n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1380 F77_FUNC(dscal,DSCAL)(n, &temp1, &workd[iwork[8]], &c__1);
1383 F77_FUNC(dlascl,DLASCL)("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
1384 v_dim1 + 1], n, &infol);
1385 F77_FUNC(dlascl,DLASCL)("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
1390 F77_FUNC(dcopy,DCOPY)(n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1391 ipntr[1] = iwork[10];
1392 ipntr[2] = iwork[9];
1393 ipntr[3] = iwork[8];
1402 F77_FUNC(dcopy,DCOPY)(n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1409 ipntr[1] = iwork[9];
1410 ipntr[2] = iwork[8];
1414 } else if (*bmat == 'I') {
1415 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1424 workd[*n * 3 + 3] =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[iwork[10]], &
1426 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1427 } else if (*bmat == 'G') {
1428 workd[*n * 3 + 3] =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &
1430 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
1431 } else if (*bmat == 'I') {
1432 workd[*n * 3 + 3] =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
1436 F77_FUNC(dgemv,DGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]]
1437 , &c__1, &c_b42, &workd[iwork[9]], &c__1);
1438 } else if (*mode == 2) {
1439 F77_FUNC(dgemv,DGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
1440 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
1443 F77_FUNC(dgemv,DGEMV)("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1444 c__1, &c_b18, &resid[1], &c__1);
1446 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1447 if (iwork[12] == 1 || iwork[4] == 1) {
1448 h__[iwork[12] + h_dim1] = 0.;
1450 h__[iwork[12] + h_dim1] = *rnorm;
1457 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1458 ipntr[1] = iwork[9];
1459 ipntr[2] = iwork[8];
1463 } else if (*bmat == 'I') {
1464 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1471 *rnorm =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1472 *rnorm = sqrt(fabs(*rnorm));
1473 } else if (*bmat == 'I') {
1474 *rnorm =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
1477 if (*rnorm > workd[*n * 3 + 3] * .717f) {
1483 F77_FUNC(dgemv,DGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
1484 c__1, &c_b42, &workd[iwork[9]], &c__1);
1486 F77_FUNC(dgemv,DGEMV)("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1487 c__1, &c_b18, &resid[1], &c__1);
1489 if (iwork[12] == 1 || iwork[4] == 1) {
1490 h__[iwork[12] + h_dim1] = 0.;
1492 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1496 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1497 ipntr[1] = iwork[9];
1498 ipntr[2] = iwork[8];
1502 } else if (*bmat == 'I') {
1503 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1509 workd[*n * 3 + 2] =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &
1511 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
1512 } else if (*bmat == 'I') {
1513 workd[*n * 3 + 2] =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
1517 if (workd[*n * 3 + 2] > *rnorm * .717f) {
1519 *rnorm = workd[*n * 3 + 2];
1523 *rnorm = workd[*n * 3 + 2];
1525 if (iwork[1] <= 1) {
1530 for (jj = 1; jj <= i__1; ++jj) {
1541 if (h__[iwork[12] + h_dim1] < 0.) {
1542 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1543 if (iwork[12] < *k + *np) {
1544 F77_FUNC(dscal,DSCAL)(n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1546 F77_FUNC(dscal,DSCAL)(n, &c_b50, &resid[1], &c__1);
1551 if (iwork[12] > *k + *np) {
1570 F77_FUNC(dsaup2,DSAUP2)(int * ido,
1600 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
1619 v_offset = 1 + v_dim1;
1622 h_offset = 1 + h_dim1;
1625 q_offset = 1 + q_dim1;
1629 eps23 = GMX_DOUBLE_EPS;
1630 eps23 = pow(eps23, c_b3);
1642 iwork[7] = iwork[9] + iwork[10];
1660 if (iwork[2] == 1) {
1661 F77_FUNC(dgetv0,DGETV0)(ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
1662 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41]
1669 if (workd[*n * 3 + 1] == 0.) {
1678 if (iwork[4] == 1) {
1682 if (iwork[5] == 1) {
1686 if (iwork[1] == 1) {
1690 F77_FUNC(dsaitr,DSAITR)(ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
1691 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1715 F77_FUNC(dsaitr,DSAITR)(ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1],
1716 &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1732 F77_FUNC(dseigt,DSEIGT)(&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
1733 bounds[1], &workl[1], &ierr);
1740 F77_FUNC(dcopy,DCOPY)(&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1741 F77_FUNC(dcopy,DCOPY)(&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
1745 F77_FUNC(dsgets,DSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
1747 F77_FUNC(dcopy,DCOPY)(nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
1748 F77_FUNC(dsconv,DSCONV)(nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
1752 for (j = 1; j <= i__1; ++j) {
1753 if (bounds[j] == 0.) {
1759 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0) {
1761 if (!strncmp(which, "BE", 2)) {
1763 strncpy(wprime, "SA",2);
1764 F77_FUNC(dsortr,DSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
1766 nevm2 = *nev - nevd2;
1768 i__1 = (nevd2 < *np) ? nevd2 : *np;
1769 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
1770 F77_FUNC(dswap,DSWAP)(&i__1, &ritz[nevm2 + 1], &c__1,
1771 &ritz[((i__2>i__3) ? i__2 : i__3)],
1773 i__1 = (nevd2 < *np) ? nevd2 : *np;
1774 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
1775 F77_FUNC(dswap,DSWAP)(&i__1, &bounds[nevm2 + 1], &c__1,
1776 &bounds[((i__2>i__3) ? i__2 : i__3) + 1],
1782 if (!strncmp(which, "LM", 2)) {
1783 strncpy(wprime, "SM", 2);
1785 if (!strncmp(which, "SM", 2)) {
1786 strncpy(wprime, "LM", 2);
1788 if (!strncmp(which, "LA", 2)) {
1789 strncpy(wprime, "SA", 2);
1791 if (!strncmp(which, "SA", 2)) {
1792 strncpy(wprime, "LA", 2);
1795 F77_FUNC(dsortr,DSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
1800 for (j = 1; j <= i__1; ++j) {
1802 d__3 = fabs(ritz[j]);
1803 temp = (d__2>d__3) ? d__2 : d__3;
1807 strncpy(wprime, "LA",2);
1808 F77_FUNC(dsortr,DSORTR)(wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
1811 for (j = 1; j <= i__1; ++j) {
1813 d__3 = fabs(ritz[j]);
1814 temp = (d__2>d__3) ? d__2 : d__3;
1818 if (!strncmp(which, "BE", 2)) {
1820 strncpy(wprime, "LA", 2);
1821 F77_FUNC(dsortr,DSORTR)(wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
1824 F77_FUNC(dsortr,DSORTR)(which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
1828 h__[h_dim1 + 1] = workd[*n * 3 + 1];
1831 if (iwork[6] > *mxiter && iwork[8] < *nev) {
1834 if (*np == 0 && iwork[8] < iwork[9]) {
1841 } else if (iwork[8] < *nev && *ishift == 1) {
1843 i__1 = iwork[8], i__2 = *np / 2;
1844 *nev += (i__1 < i__2) ? i__1 : i__2;
1845 if (*nev == 1 && iwork[7] >= 6) {
1846 *nev = iwork[7] / 2;
1847 } else if (*nev == 1 && iwork[7] > 2) {
1850 *np = iwork[7] - *nev;
1853 if (nevbef < *nev) {
1854 F77_FUNC(dsgets,DSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
1872 F77_FUNC(dcopy,DCOPY)(np, &workl[1], &c__1, &ritz[1], &c__1);
1875 F77_FUNC(dsapps,DSAPPS)(n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
1876 resid[1], &q[q_offset], ldq, &workd[1]);
1880 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
1886 } else if (*bmat == 'I') {
1887 F77_FUNC(dcopy,DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
1893 workd[*n * 3 + 1] =F77_FUNC(ddot,DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
1894 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
1895 } else if (*bmat == 'I') {
1896 workd[*n * 3 + 1] =F77_FUNC(dnrm2,DNRM2)(n, &resid[1], &c__1);
1918 F77_FUNC(dsaupd,DSAUPD)(int * ido,
1936 int v_dim1, v_offset, i__1, i__2;
1942 v_offset = 1 + v_dim1;
1953 iwork[5] = iparam[1];
1954 iwork[10] = iparam[3];
1955 iwork[12] = iparam[4];
1958 iwork[11] = iparam[7];
1963 } else if (*nev <= 0) {
1965 } else if (*ncv <= *nev || *ncv > *n) {
1970 iwork[15] = *ncv - *nev;
1972 if (iwork[10] <= 0) {
1975 if (strncmp(which,"LM",2) && strncmp(which,"SM",2) &&
1976 strncmp(which,"LA",2) && strncmp(which,"SA",2) &&
1977 strncmp(which,"BE",2)) {
1980 if (*bmat != 'I' && *bmat != 'G') {
1985 if (*lworkl < i__1 * i__1 + (*ncv << 3)) {
1988 if (iwork[11] < 1 || iwork[11] > 5) {
1990 } else if (iwork[11] == 1 && *bmat == 'G') {
1992 } else if (iwork[5] < 0 || iwork[5] > 1) {
1994 } else if (*nev == 1 && !strncmp(which, "BE", 2)) {
1998 if (iwork[2] != 0) {
2004 if (iwork[12] <= 0) {
2008 *tol = GMX_DOUBLE_EPS;
2011 iwork[15] = *ncv - *nev;
2014 i__1 = i__2 * i__2 + (*ncv << 3);
2015 for (j = 1; j <= i__1; ++j) {
2022 iwork[16] = iwork[3] + (iwork[8] << 1);
2023 iwork[1] = iwork[16] + *ncv;
2024 iwork[4] = iwork[1] + *ncv;
2026 iwork[7] = iwork[4] + i__1 * i__1;
2027 iwork[14] = iwork[7] + *ncv * 3;
2029 ipntr[4] = iwork[14];
2030 ipntr[5] = iwork[3];
2031 ipntr[6] = iwork[16];
2032 ipntr[7] = iwork[1];
2033 ipntr[11] = iwork[7];
2036 F77_FUNC(dsaup2,DSAUP2)(ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
2037 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
2038 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
2039 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1]
2040 , &iwork[21], info);
2043 iparam[8] = iwork[15];
2049 iparam[3] = iwork[10];
2050 iparam[5] = iwork[15];
2068 F77_FUNC(dseupd,DSEUPD)(int * rvec,
2069 const char * howmny,
2094 int v_dim1, v_offset, z_dim1, z_offset, i__1;
2095 double d__1, d__2, d__3;
2097 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2109 double thres1=0, thres2=0;
2113 int leftptr, rghtptr;
2119 z_offset = 1 + z_dim1;
2124 v_offset = 1 + v_dim1;
2148 if (*ncv <= *nev || *ncv > *n) {
2151 if (strncmp(which,"LM",2) && strncmp(which,"SM",2) &&
2152 strncmp(which,"LA",2) && strncmp(which,"SA",2) &&
2153 strncmp(which,"BE",2)) {
2156 if (*bmat != 'I' && *bmat != 'G') {
2159 if (*howmny != 'A' && *howmny != 'P' &&
2160 *howmny != 'S' && *rvec) {
2163 if (*rvec && *howmny == 'S') {
2167 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3)) {
2171 if (mode == 1 || mode == 2) {
2172 strncpy(type__, "REGULR",6);
2173 } else if (mode == 3) {
2174 strncpy(type__, "SHIFTI",6);
2175 } else if (mode == 4) {
2176 strncpy(type__, "BUCKLE",6);
2177 } else if (mode == 5) {
2178 strncpy(type__, "CAYLEY",6);
2182 if (mode == 1 && *bmat == 'G') {
2185 if (*nev == 1 && !strncmp(which, "BE",2)) {
2202 iw = iq + ldh * *ncv;
2203 next = iw + (*ncv << 1);
2209 irz = ipntr[11] + *ncv;
2213 eps23 = GMX_DOUBLE_EPS;
2214 eps23 = pow(eps23, c_b21);
2219 } else if (*bmat == 'G') {
2220 bnorm2 =F77_FUNC(dnrm2,DNRM2)(n, &workd[1], &c__1);
2225 if (!strncmp(which,"LM",2) || !strncmp(which,"SM",2) ||
2226 !strncmp(which,"LA",2) || !strncmp(which,"SA",2)) {
2228 } else if (!strncmp(which,"BE",2)) {
2231 ism = (*nev>nconv) ? *nev : nconv;
2234 thres1 = workl[ism];
2235 thres2 = workl[ilg];
2243 for (j = 0; j <= i__1; ++j) {
2245 if (!strncmp(which,"LM",2)) {
2246 if (fabs(workl[irz + j]) >= fabs(thres1)) {
2248 d__3 = fabs(workl[irz + j]);
2249 tempbnd = (d__2>d__3) ? d__2 : d__3;
2250 if (workl[ibd + j] <= *tol * tempbnd) {
2254 } else if (!strncmp(which,"SM",2)) {
2255 if (fabs(workl[irz + j]) <= fabs(thres1)) {
2257 d__3 = fabs(workl[irz + j]);
2258 tempbnd = (d__2>d__3) ? d__2 : d__3;
2259 if (workl[ibd + j] <= *tol * tempbnd) {
2263 } else if (!strncmp(which,"LA",2)) {
2264 if (workl[irz + j] >= thres1) {
2266 d__3 = fabs(workl[irz + j]);
2267 tempbnd = (d__2>d__3) ? d__2 : d__3;
2268 if (workl[ibd + j] <= *tol * tempbnd) {
2272 } else if (!strncmp(which,"SA",2)) {
2273 if (workl[irz + j] <= thres1) {
2275 d__3 = fabs(workl[irz + j]);
2276 tempbnd = (d__2>d__3) ? d__2 : d__3;
2277 if (workl[ibd + j] <= *tol * tempbnd) {
2281 } else if (!strncmp(which,"BE",2)) {
2282 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2) {
2284 d__3 = fabs(workl[irz + j]);
2285 tempbnd = (d__2>d__3) ? d__2 : d__3;
2286 if (workl[ibd + j] <= *tol * tempbnd) {
2291 if (j + 1 > nconv) {
2292 reord = select[j + 1] || reord;
2294 if (select[j + 1]) {
2300 F77_FUNC(dcopy,DCOPY)(&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2301 F77_FUNC(dcopy,DCOPY)(ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2303 F77_FUNC(dsteqr,DSTEQR)("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
2322 if (select[leftptr]) {
2326 } else if (! select[rghtptr]) {
2332 temp = workl[ihd + leftptr - 1];
2333 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2334 workl[ihd + rghtptr - 1] = temp;
2335 F77_FUNC(dcopy,DCOPY)(ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
2337 F77_FUNC(dcopy,DCOPY)(ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
2338 iq + *ncv * (leftptr - 1)], &c__1);
2339 F77_FUNC(dcopy,DCOPY)(ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
2346 if (leftptr < rghtptr) {
2354 F77_FUNC(dcopy,DCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2358 F77_FUNC(dcopy,DCOPY)(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2359 F77_FUNC(dcopy,DCOPY)(ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2362 if (!strncmp(type__, "REGULR",6)) {
2365 F77_FUNC(dsesrt,DSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2367 F77_FUNC(dcopy,DCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2372 F77_FUNC(dcopy,DCOPY)(ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2373 if (!strncmp(type__, "SHIFTI", 6)) {
2375 for (k = 1; k <= i__1; ++k) {
2376 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2378 } else if (!strncmp(type__, "BUCKLE",6)) {
2380 for (k = 1; k <= i__1; ++k) {
2381 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
2384 } else if (!strncmp(type__, "CAYLEY",6)) {
2386 for (k = 1; k <= i__1; ++k) {
2387 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
2388 workl[ihd + k - 1] - 1.);
2392 F77_FUNC(dcopy,DCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2393 F77_FUNC(dsortr,DSORTR)("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2395 F77_FUNC(dsesrt,DSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2397 F77_FUNC(dcopy,DCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2398 d__1 = bnorm2 / rnorm;
2399 F77_FUNC(dscal,DSCAL)(ncv, &d__1, &workl[ihb], &c__1);
2400 F77_FUNC(dsortr,DSORTR)("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2405 if (*rvec && *howmny == 'A') {
2407 F77_FUNC(dgeqr2,DGEQR2)(ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2410 F77_FUNC(dorm2r,DORM2R)("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
2411 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
2412 F77_FUNC(dlacpy,DLACPY)("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2415 for (j = 1; j <= i__1; ++j) {
2416 workl[ihb + j - 1] = 0.;
2418 workl[ihb + *ncv - 1] = 1.;
2419 F77_FUNC(dorm2r,DORM2R)("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
2420 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
2422 } else if (*rvec && *howmny == 'S') {
2426 if (!strncmp(type__, "REGULR",6) && *rvec) {
2429 for (j = 1; j <= i__1; ++j) {
2430 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
2433 } else if (strncmp(type__, "REGULR",6) && *rvec) {
2435 F77_FUNC(dscal,DSCAL)(ncv, &bnorm2, &workl[ihb], &c__1);
2436 if (!strncmp(type__, "SHIFTI",6)) {
2439 for (k = 1; k <= i__1; ++k) {
2440 d__2 = workl[iw + k - 1];
2441 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
2444 } else if (!strncmp(type__, "BUCKLE",6)) {
2447 for (k = 1; k <= i__1; ++k) {
2448 d__2 = workl[iw + k - 1] - 1.;
2449 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
2452 } else if (!strncmp(type__, "CAYLEY",6)) {
2455 for (k = 1; k <= i__1; ++k) {
2456 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2464 if (*rvec && (!strncmp(type__, "SHIFTI",6) || !strncmp(type__, "CAYLEY",6))) {
2467 for (k = 0; k <= i__1; ++k) {
2468 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2471 } else if (*rvec && !strncmp(type__, "BUCKLE", 6)) {
2474 for (k = 0; k <= i__1; ++k) {
2475 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
2481 if (strncmp(type__, "REGULR",6)) {
2482 F77_FUNC(dger,DGER)(n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
2496 /* Selected single precision arpack routines */
2500 F77_FUNC(sstqrb,SSTQRB)(int * n,
2514 int i__, j, k, l, m;
2516 int l1, ii, mm, lm1, mm1, nm1;
2517 float rt1, rt2, eps;
2520 int lend, jtot, lendm1, lendp1, iscale;
2522 int lendsv, nmaxit, icompz;
2523 float ssfmax, ssfmin,safmin,minval,safmax,anorm;
2546 eps = GMX_FLOAT_EPS;
2550 minval = GMX_FLOAT_MIN;
2551 safmin = minval / GMX_FLOAT_EPS;
2552 safmax = 1. / safmin;
2553 ssfmax = sqrt(safmax) / 3.;
2554 ssfmin = sqrt(safmin) / eps2;
2558 for (j = 1; j <= i__1; ++j) {
2580 for (m = l1; m <= i__1; ++m) {
2585 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m+1])) * eps) {
2603 i__1 = lend - l + 1;
2604 anorm =F77_FUNC(slanst,SLANST)("i", &i__1, &d__[l], &e[l]);
2609 if (anorm > ssfmax) {
2611 i__1 = lend - l + 1;
2612 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
2615 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
2617 } else if (anorm < ssfmin) {
2619 i__1 = lend - l + 1;
2620 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
2623 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
2627 if (fabs(d__[lend]) < fabs(d__[l])) {
2638 for (m = l; m <= i__1; ++m) {
2641 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m + 1]) + safmin) {
2660 F77_FUNC(slaev2,SLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
2662 work[*n - 1 + l] = s;
2665 z__[l + 1] = c__ * tst - s * z__[l];
2666 z__[l] = s * tst + c__ * z__[l];
2668 F77_FUNC(slae2,SLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
2680 if (jtot == nmaxit) {
2685 g = (d__[l + 1] - p) / (e[l] * 2.);
2686 r__ =F77_FUNC(slapy2,SLAPY2)(&g, &c_b31);
2687 g = d__[m] - p + e[l] / (g + ((g>0) ? r__ : -r__ ));
2695 for (i__ = mm1; i__ >= i__1; --i__) {
2698 F77_FUNC(slartg,SLARTG)(&g, &f, &c__, &s, &r__);
2702 g = d__[i__ + 1] - p;
2703 r__ = (d__[i__] - g) * s + c__ * 2. * b;
2705 d__[i__ + 1] = g + p;
2710 work[*n - 1 + i__] = -s;
2718 F77_FUNC(slasr,SLASR)("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
2741 for (m = l; m >= i__1; --m) {
2742 d__2 = fabs(e[m - 1]);
2744 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin) {
2763 F77_FUNC(slaev2,SLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
2767 z__[l] = c__ * tst - s * z__[l - 1];
2768 z__[l - 1] = s * tst + c__ * z__[l - 1];
2770 F77_FUNC(slae2,SLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
2782 if (jtot == nmaxit) {
2788 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
2789 r__ =F77_FUNC(slapy2,SLAPY2)(&g, &c_b31);
2790 g = d__[m] - p + e[l - 1] / (g + ((g>0) ? r__ : -r__ ));
2798 for (i__ = m; i__ <= i__1; ++i__) {
2801 F77_FUNC(slartg,SLARTG)(&g, &f, &c__, &s, &r__);
2806 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
2813 work[*n - 1 + i__] = s;
2821 F77_FUNC(slasr,SLASR)("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
2842 i__1 = lendsv - lsv + 1;
2843 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
2845 i__1 = lendsv - lsv;
2846 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
2848 } else if (iscale == 2) {
2849 i__1 = lendsv - lsv + 1;
2850 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
2852 i__1 = lendsv - lsv;
2853 F77_FUNC(slascl,SLASCL)("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
2857 if (jtot < nmaxit) {
2861 for (i__ = 1; i__ <= i__1; ++i__) {
2871 F77_FUNC(slasrt,SLASRT)("i", n, &d__[1], info);
2876 for (ii = 2; ii <= i__1; ++ii) {
2881 for (j = ii; j <= i__2; ++j) {
2904 F77_FUNC(sgetv0,SGETV0)(int * ido,
2923 int v_dim1, v_offset, i__1;
2931 v_offset = 1 + v_dim1;
2945 F77_FUNC(slarnv,SLARNV)(&idist, &iwork[1], n, &resid[1]);
2951 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
2957 if (iwork[5] == 1) {
2961 if (iwork[6] == 1) {
2967 F77_FUNC(scopy,SCOPY)(n, &workd[*n + 1], &c__1, &resid[1], &c__1);
2972 } else if (*bmat == 'I') {
2973 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
2981 workd[*n * 3 + 4] =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
2982 workd[*n * 3 + 4] = sqrt(fabs(workd[*n * 3 + 4]));
2983 } else if (*bmat == 'I') {
2984 workd[*n * 3 + 4] =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
2986 *rnorm = workd[*n * 3 + 4];
2995 F77_FUNC(sgemv,SGEMV)("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
2996 &workd[*n + 1], &c__1);
2998 F77_FUNC(sgemv,SGEMV)("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
2999 c_b22, &resid[1], &c__1);
3002 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3007 } else if (*bmat == 'I') {
3008 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
3014 *rnorm =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
3015 *rnorm = sqrt(fabs(*rnorm));
3016 } else if (*bmat == 'I') {
3017 *rnorm =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
3020 if (*rnorm > workd[*n * 3 + 4] * .717f) {
3025 if (iwork[7] <= 1) {
3027 workd[*n * 3 + 4] = *rnorm;
3032 for (jj = 1; jj <= i__1; ++jj) {
3052 F77_FUNC(ssapps,SSAPPS)(int * n,
3069 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
3073 float r__, s, a1, a2, a3, a4;
3084 v_offset = 1 + v_dim1;
3087 h_offset = 1 + h_dim1;
3090 q_offset = 1 + q_dim1;
3093 epsmch = GMX_FLOAT_EPS;
3097 kplusp = *kev + *np;
3099 F77_FUNC(slaset,SLASET)("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3106 for (jj = 1; jj <= i__1; ++jj) {
3113 for (i__ = istart; i__ <= i__2; ++i__) {
3114 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__ + 1 + (h_dim1*2)]);
3115 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
3116 h__[i__ + 1 + h_dim1] = 0.;
3124 if (istart < iend) {
3126 f = h__[istart + (h_dim1 << 1)] - shift[jj];
3127 g = h__[istart + 1 + h_dim1];
3128 F77_FUNC(slartg,SLARTG)(&f, &g, &c__, &s, &r__);
3130 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
3132 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
3134 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
3136 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
3138 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
3139 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3140 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
3143 i__2 = (i__3<kplusp) ? i__3 : kplusp;
3144 for (j = 1; j <= i__2; ++j) {
3145 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
3147 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
3148 c__ * q[j + (istart + 1) * q_dim1];
3149 q[j + istart * q_dim1] = a1;
3154 for (i__ = istart + 1; i__ <= i__2; ++i__) {
3156 f = h__[i__ + h_dim1];
3157 g = s * h__[i__ + 1 + h_dim1];
3159 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3160 F77_FUNC(slartg,SLARTG)(&f, &g, &c__, &s, &r__);
3168 h__[i__ + h_dim1] = r__;
3170 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
3172 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
3174 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
3176 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
3179 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
3180 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3181 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
3184 i__3 = (i__4<kplusp) ? i__4 : kplusp;
3185 for (j = 1; j <= i__3; ++j) {
3186 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
3188 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
3189 c__ * q[j + (i__ + 1) * q_dim1];
3190 q[j + i__ * q_dim1] = a1;
3199 if (h__[iend + h_dim1] < 0.) {
3200 h__[iend + h_dim1] = -h__[iend + h_dim1];
3201 F77_FUNC(sscal,SSCAL)(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3204 if (iend < kplusp) {
3209 for (i__ = itop; i__ <= i__2; ++i__) {
3210 if (h__[i__ + 1 + h_dim1] > 0.) {
3221 for (i__ = itop; i__ <= i__1; ++i__) {
3222 big = fabs(h__[i__ + (h_dim1*2)]) + fabs(h__[i__+ 1 + (h_dim1*2)]);
3223 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
3224 h__[i__ + 1 + h_dim1] = 0.;
3229 if (h__[*kev + 1 + h_dim1] > 0.) {
3230 F77_FUNC(sgemv,SGEMV)("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
3231 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
3235 for (i__ = 1; i__ <= i__1; ++i__) {
3236 i__2 = kplusp - i__ + 1;
3237 F77_FUNC(sgemv,SGEMV)("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
3238 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
3239 F77_FUNC(scopy,SCOPY)(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
3244 F77_FUNC(slacpy,SLACPY)("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3246 if (h__[*kev + 1 + h_dim1] > 0.) {
3247 F77_FUNC(scopy,SCOPY)(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3250 F77_FUNC(sscal,SSCAL)(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3251 if (h__[*kev + 1 + h_dim1] > 0.) {
3252 F77_FUNC(saxpy,SAXPY)(n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
3266 F77_FUNC(ssortr,SSORTR)(const char * which,
3281 if (!strncmp(which, "SA", 2)) {
3288 for (i__ = igap; i__ <= i__1; ++i__) {
3296 if (x1[j] < x1[j + igap]) {
3298 x1[j] = x1[j + igap];
3299 x1[j + igap] = temp;
3302 x2[j] = x2[j + igap];
3303 x2[j + igap] = temp;
3316 } else if (!strncmp(which, "SM", 2)) {
3323 for (i__ = igap; i__ <= i__1; ++i__) {
3331 if (fabs(x1[j]) < fabs(x1[j + igap])) {
3333 x1[j] = x1[j + igap];
3334 x1[j + igap] = temp;
3337 x2[j] = x2[j + igap];
3338 x2[j + igap] = temp;
3351 } else if (!strncmp(which, "LA", 2)) {
3358 for (i__ = igap; i__ <= i__1; ++i__) {
3366 if (x1[j] > x1[j + igap]) {
3368 x1[j] = x1[j + igap];
3369 x1[j + igap] = temp;
3372 x2[j] = x2[j + igap];
3373 x2[j + igap] = temp;
3386 } else if (!strncmp(which, "LM", 2)) {
3394 for (i__ = igap; i__ <= i__1; ++i__) {
3402 if (fabs(x1[j]) > fabs(x1[j + igap])) {
3404 x1[j] = x1[j + igap];
3405 x1[j + igap] = temp;
3408 x2[j] = x2[j + igap];
3409 x2[j + igap] = temp;
3432 F77_FUNC(ssesrt,SSESRT)(const char * which,
3440 int a_dim1, a_offset, i__1;
3447 a_offset = 1 + a_dim1 * 0;
3452 if (!strncmp(which, "SA", 2)) {
3459 for (i__ = igap; i__ <= i__1; ++i__) {
3467 if (x[j] < x[j + igap]) {
3472 F77_FUNC(sswap,SSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
3473 a_dim1 + 1], &c__1);
3486 } else if (!strncmp(which, "SM", 2)) {
3493 for (i__ = igap; i__ <= i__1; ++i__) {
3501 if (fabs(x[j]) < fabs(x[j + igap])) {
3506 F77_FUNC(sswap,SSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
3507 a_dim1 + 1], &c__1);
3520 } else if (!strncmp(which, "LA", 2)) {
3527 for (i__ = igap; i__ <= i__1; ++i__) {
3535 if (x[j] > x[j + igap]) {
3540 F77_FUNC(sswap,SSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
3541 a_dim1 + 1], &c__1);
3554 } else if (!strncmp(which, "LM", 2)) {
3561 for (i__ = igap; i__ <= i__1; ++i__) {
3569 if (fabs(x[j]) > fabs(x[j + igap])) {
3574 F77_FUNC(sswap,SSWAP)(na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
3575 a_dim1 + 1], &c__1);
3598 F77_FUNC(ssgets,SSGETS)(int * ishift,
3614 if (!strncmp(which, "BE", 2)) {
3616 F77_FUNC(ssortr,SSORTR)("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
3619 i__1 = (kevd2<*np) ? kevd2 : *np;
3620 i__2 = (kevd2>*np) ? kevd2 : *np;
3621 F77_FUNC(sswap,SSWAP)(&i__1, &ritz[1], &c__1,
3622 &ritz[i__2 + 1], &c__1);
3623 i__1 = (kevd2<*np) ? kevd2 : *np;
3624 i__2 = (kevd2>*np) ? kevd2 : *np;
3625 F77_FUNC(sswap,SSWAP)(&i__1, &bounds[1], &c__1,
3626 &bounds[i__2 + 1], &c__1);
3631 F77_FUNC(ssortr,SSORTR)(which, &c__1, &i__1, &ritz[1], &bounds[1]);
3634 if (*ishift == 1 && *np > 0) {
3636 F77_FUNC(ssortr,SSORTR)("SM", &c__1, np, &bounds[1], &ritz[1]);
3637 F77_FUNC(scopy,SCOPY)(np, &ritz[1], &c__1, &shifts[1], &c__1);
3647 F77_FUNC(ssconv,SSCONV)(int * n,
3663 eps23 = GMX_FLOAT_EPS;
3664 eps23 = pow(eps23, c_b3);
3668 for (i__ = 1; i__ <= i__1; ++i__) {
3671 d__3 = fabs(ritz[i__]);
3672 temp = (d__2 > d__3) ? d__2 : d__3;
3673 if (bounds[i__] <= *tol * temp) {
3683 F77_FUNC(sseigt,SSEIGT)(float * rnorm,
3693 int h_dim1, h_offset, i__1;
3702 h_offset = 1 + h_dim1;
3705 F77_FUNC(scopy,SCOPY)(n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
3707 F77_FUNC(scopy,SCOPY)(&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
3708 F77_FUNC(sstqrb,SSTQRB)(n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
3714 for (k = 1; k <= i__1; ++k) {
3715 bounds[k] = *rnorm * fabs(bounds[k]);
3728 F77_FUNC(ssaitr,SSAITR)(int * ido,
3752 int h_dim1, h_offset, v_dim1, v_offset, i__1;
3756 float safmin,minval;
3762 v_offset = 1 + v_dim1;
3765 h_offset = 1 + h_dim1;
3769 minval = GMX_FLOAT_MIN;
3770 safmin = minval / GMX_FLOAT_EPS;
3783 iwork[9] = iwork[8] + *n;
3784 iwork[10] = iwork[9] + *n;
3787 if (iwork[5] == 1) {
3790 if (iwork[6] == 1) {
3793 if (iwork[2] == 1) {
3796 if (iwork[3] == 1) {
3799 if (iwork[4] == 1) {
3816 F77_FUNC(sgetv0,sgetv0)(ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
3817 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
3823 if (iwork[11] <= 3) {
3827 *info = iwork[12] - 1;
3834 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
3835 if (*rnorm >= safmin) {
3836 temp1 = 1. / *rnorm;
3837 F77_FUNC(sscal,SSCAL)(n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
3838 F77_FUNC(sscal,SSCAL)(n, &temp1, &workd[iwork[8]], &c__1);
3841 F77_FUNC(slascl,SLASCL)("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
3842 v_dim1 + 1], n, &infol);
3843 F77_FUNC(slascl,SLASCL)("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
3848 F77_FUNC(scopy,SCOPY)(n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
3849 ipntr[1] = iwork[10];
3850 ipntr[2] = iwork[9];
3851 ipntr[3] = iwork[8];
3860 F77_FUNC(scopy,SCOPY)(n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
3867 ipntr[1] = iwork[9];
3868 ipntr[2] = iwork[8];
3872 } else if (*bmat == 'I') {
3873 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
3882 workd[*n * 3 + 3] =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[iwork[10]], &
3884 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
3885 } else if (*bmat == 'G') {
3886 workd[*n * 3 + 3] =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &
3888 workd[*n * 3 + 3] = sqrt(fabs(workd[*n * 3 + 3]));
3889 } else if (*bmat == 'I') {
3890 workd[*n * 3 + 3] =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
3894 F77_FUNC(sgemv,SGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]]
3895 , &c__1, &c_b42, &workd[iwork[9]], &c__1);
3896 } else if (*mode == 2) {
3897 F77_FUNC(sgemv,SGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
3898 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
3901 F77_FUNC(sgemv,SGEMV)("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
3902 c__1, &c_b18, &resid[1], &c__1);
3904 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
3905 if (iwork[12] == 1 || iwork[4] == 1) {
3906 h__[iwork[12] + h_dim1] = 0.;
3908 h__[iwork[12] + h_dim1] = *rnorm;
3915 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
3916 ipntr[1] = iwork[9];
3917 ipntr[2] = iwork[8];
3921 } else if (*bmat == 'I') {
3922 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
3929 *rnorm =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
3930 *rnorm = sqrt(fabs(*rnorm));
3931 } else if (*bmat == 'I') {
3932 *rnorm =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
3935 if (*rnorm > workd[*n * 3 + 3] * .717f) {
3941 F77_FUNC(sgemv,SGEMV)("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
3942 c__1, &c_b42, &workd[iwork[9]], &c__1);
3944 F77_FUNC(sgemv,SGEMV)("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
3945 c__1, &c_b18, &resid[1], &c__1);
3947 if (iwork[12] == 1 || iwork[4] == 1) {
3948 h__[iwork[12] + h_dim1] = 0.;
3950 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
3954 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
3955 ipntr[1] = iwork[9];
3956 ipntr[2] = iwork[8];
3960 } else if (*bmat == 'I') {
3961 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
3967 workd[*n * 3 + 2] =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &
3969 workd[*n * 3 + 2] = sqrt(fabs(workd[*n * 3 + 2]));
3970 } else if (*bmat == 'I') {
3971 workd[*n * 3 + 2] =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
3975 if (workd[*n * 3 + 2] > *rnorm * .717f) {
3977 *rnorm = workd[*n * 3 + 2];
3981 *rnorm = workd[*n * 3 + 2];
3983 if (iwork[1] <= 1) {
3988 for (jj = 1; jj <= i__1; ++jj) {
3999 if (h__[iwork[12] + h_dim1] < 0.) {
4000 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4001 if (iwork[12] < *k + *np) {
4002 F77_FUNC(sscal,SSCAL)(n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4004 F77_FUNC(sscal,SSCAL)(n, &c_b50, &resid[1], &c__1);
4009 if (iwork[12] > *k + *np) {
4028 F77_FUNC(ssaup2,SSAUP2)(int * ido,
4058 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
4077 v_offset = 1 + v_dim1;
4080 h_offset = 1 + h_dim1;
4083 q_offset = 1 + q_dim1;
4087 eps23 = GMX_FLOAT_EPS;
4088 eps23 = pow(eps23, c_b3);
4100 iwork[7] = iwork[9] + iwork[10];
4118 if (iwork[2] == 1) {
4119 F77_FUNC(sgetv0,SGETV0)(ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
4120 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41]
4127 if (workd[*n * 3 + 1] == 0.) {
4136 if (iwork[4] == 1) {
4140 if (iwork[5] == 1) {
4144 if (iwork[1] == 1) {
4148 F77_FUNC(ssaitr,SSAITR)(ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
4149 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
4173 F77_FUNC(ssaitr,SSAITR)(ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[
4174 v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[
4190 F77_FUNC(sseigt,SSEIGT)(&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
4191 bounds[1], &workl[1], &ierr);
4198 F77_FUNC(scopy,SCOPY)(&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4199 F77_FUNC(scopy,SCOPY)(&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4203 F77_FUNC(ssgets,SSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4205 F77_FUNC(scopy,SCOPY)(nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4206 F77_FUNC(ssconv,SSCONV)(nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4211 for (j = 1; j <= i__1; ++j) {
4212 if (bounds[j] == 0.) {
4218 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0) {
4220 if (!strncmp(which, "BE", 2)) {
4222 strncpy(wprime, "SA",2);
4223 F77_FUNC(ssortr,SSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4225 nevm2 = *nev - nevd2;
4227 i__1 = (nevd2 < *np) ? nevd2 : *np;
4228 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4229 F77_FUNC(sswap,SSWAP)(&i__1, &ritz[nevm2 + 1], &c__1,
4230 &ritz[((i__2>i__3) ? i__2 : i__3)],
4232 i__1 = (nevd2 < *np) ? nevd2 : *np;
4233 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4234 F77_FUNC(sswap,SSWAP)(&i__1, &bounds[nevm2 + 1], &c__1,
4235 &bounds[((i__2>i__3) ? i__2 : i__3) + 1],
4241 if (!strncmp(which, "LM", 2)) {
4242 strncpy(wprime, "SM", 2);
4244 if (!strncmp(which, "SM", 2)) {
4245 strncpy(wprime, "LM", 2);
4247 if (!strncmp(which, "LA", 2)) {
4248 strncpy(wprime, "SA", 2);
4250 if (!strncmp(which, "SA", 2)) {
4251 strncpy(wprime, "LA", 2);
4254 F77_FUNC(ssortr,SSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4259 for (j = 1; j <= i__1; ++j) {
4261 d__3 = fabs(ritz[j]);
4262 temp = (d__2>d__3) ? d__2 : d__3;
4266 strncpy(wprime, "LA",2);
4267 F77_FUNC(ssortr,SSORTR)(wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4270 for (j = 1; j <= i__1; ++j) {
4272 d__3 = fabs(ritz[j]);
4273 temp = (d__2>d__3) ? d__2 : d__3;
4277 if (!strncmp(which, "BE", 2)) {
4279 strncpy(wprime, "LA", 2);
4280 F77_FUNC(ssortr,SSORTR)(wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4283 F77_FUNC(ssortr,SSORTR)(which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4287 h__[h_dim1 + 1] = workd[*n * 3 + 1];
4290 if (iwork[6] > *mxiter && iwork[8] < *nev) {
4293 if (*np == 0 && iwork[8] < iwork[9]) {
4300 } else if (iwork[8] < *nev && *ishift == 1) {
4302 i__1 = iwork[8], i__2 = *np / 2;
4303 *nev += (i__1 < i__2) ? i__1 : i__2;
4304 if (*nev == 1 && iwork[7] >= 6) {
4305 *nev = iwork[7] / 2;
4306 } else if (*nev == 1 && iwork[7] > 2) {
4309 *np = iwork[7] - *nev;
4312 if (nevbef < *nev) {
4313 F77_FUNC(ssgets,SSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4331 F77_FUNC(scopy,SCOPY)(np, &workl[1], &c__1, &ritz[1], &c__1);
4334 F77_FUNC(ssapps,SSAPPS)(n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
4335 resid[1], &q[q_offset], ldq, &workd[1]);
4339 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
4345 } else if (*bmat == 'I') {
4346 F77_FUNC(scopy,SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
4352 workd[*n * 3 + 1] =F77_FUNC(sdot,SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
4353 workd[*n * 3 + 1] = sqrt(fabs(workd[*n * 3 + 1]));
4354 } else if (*bmat == 'I') {
4355 workd[*n * 3 + 1] =F77_FUNC(snrm2,SNRM2)(n, &resid[1], &c__1);
4377 F77_FUNC(ssaupd,SSAUPD)(int * ido,
4395 int v_dim1, v_offset, i__1, i__2;
4401 v_offset = 1 + v_dim1;
4412 iwork[5] = iparam[1];
4413 iwork[10] = iparam[3];
4414 iwork[12] = iparam[4];
4417 iwork[11] = iparam[7];
4422 } else if (*nev <= 0) {
4424 } else if (*ncv <= *nev || *ncv > *n) {
4429 iwork[15] = *ncv - *nev;
4431 if (iwork[10] <= 0) {
4434 if (strncmp(which,"LM",2) && strncmp(which,"SM",2) &&
4435 strncmp(which,"LA",2) && strncmp(which,"SA",2) &&
4436 strncmp(which,"BE",2)) {
4439 if (*bmat != 'I' && *bmat != 'G') {
4444 if (*lworkl < i__1 * i__1 + (*ncv << 3)) {
4447 if (iwork[11] < 1 || iwork[11] > 5) {
4449 } else if (iwork[11] == 1 && *bmat == 'G') {
4451 } else if (iwork[5] < 0 || iwork[5] > 1) {
4453 } else if (*nev == 1 && !strncmp(which, "BE", 2)) {
4457 if (iwork[2] != 0) {
4463 if (iwork[12] <= 0) {
4467 *tol = GMX_FLOAT_EPS;
4470 iwork[15] = *ncv - *nev;
4473 i__1 = i__2 * i__2 + (*ncv << 3);
4474 for (j = 1; j <= i__1; ++j) {
4481 iwork[16] = iwork[3] + (iwork[8] << 1);
4482 iwork[1] = iwork[16] + *ncv;
4483 iwork[4] = iwork[1] + *ncv;
4485 iwork[7] = iwork[4] + i__1 * i__1;
4486 iwork[14] = iwork[7] + *ncv * 3;
4488 ipntr[4] = iwork[14];
4489 ipntr[5] = iwork[3];
4490 ipntr[6] = iwork[16];
4491 ipntr[7] = iwork[1];
4492 ipntr[11] = iwork[7];
4495 F77_FUNC(ssaup2,SSAUP2)(ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
4496 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
4497 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
4498 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1]
4499 , &iwork[21], info);
4502 iparam[8] = iwork[15];
4508 iparam[3] = iwork[10];
4509 iparam[5] = iwork[15];
4527 F77_FUNC(sseupd,SSEUPD)(int * rvec,
4528 const char * howmny,
4553 int v_dim1, v_offset, z_dim1, z_offset, i__1;
4554 float d__1, d__2, d__3;
4556 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
4568 float thres1=0, thres2=0;
4572 int leftptr, rghtptr;
4578 z_offset = 1 + z_dim1;
4583 v_offset = 1 + v_dim1;
4607 if (*ncv <= *nev || *ncv > *n) {
4610 if (strncmp(which,"LM",2) && strncmp(which,"SM",2) &&
4611 strncmp(which,"LA",2) && strncmp(which,"SA",2) &&
4612 strncmp(which,"BE",2)) {
4615 if (*bmat != 'I' && *bmat != 'G') {
4618 if (*howmny != 'A' && *howmny != 'P' &&
4619 *howmny != 'S' && *rvec) {
4622 if (*rvec && *howmny == 'S') {
4626 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3)) {
4630 if (mode == 1 || mode == 2) {
4631 strncpy(type__, "REGULR",6);
4632 } else if (mode == 3) {
4633 strncpy(type__, "SHIFTI",6);
4634 } else if (mode == 4) {
4635 strncpy(type__, "BUCKLE",6);
4636 } else if (mode == 5) {
4637 strncpy(type__, "CAYLEY",6);
4641 if (mode == 1 && *bmat == 'G') {
4644 if (*nev == 1 && !strncmp(which, "BE",2)) {
4661 iw = iq + ldh * *ncv;
4662 next = iw + (*ncv << 1);
4668 irz = ipntr[11] + *ncv;
4672 eps23 = GMX_FLOAT_EPS;
4673 eps23 = pow(eps23, c_b21);
4678 } else if (*bmat == 'G') {
4679 bnorm2 =F77_FUNC(snrm2,SNRM2)(n, &workd[1], &c__1);
4684 if (!strncmp(which,"LM",2) || !strncmp(which,"SM",2) ||
4685 !strncmp(which,"LA",2) || !strncmp(which,"SA",2)) {
4687 } else if (!strncmp(which,"BE",2)) {
4690 ism = (*nev>nconv) ? *nev : nconv;
4693 thres1 = workl[ism];
4694 thres2 = workl[ilg];
4702 for (j = 0; j <= i__1; ++j) {
4704 if (!strncmp(which,"LM",2)) {
4705 if (fabs(workl[irz + j]) >= fabs(thres1)) {
4707 d__3 = fabs(workl[irz + j]);
4708 tempbnd = (d__2>d__3) ? d__2 : d__3;
4709 if (workl[ibd + j] <= *tol * tempbnd) {
4713 } else if (!strncmp(which,"SM",2)) {
4714 if (fabs(workl[irz + j]) <= fabs(thres1)) {
4716 d__3 = fabs(workl[irz + j]);
4717 tempbnd = (d__2>d__3) ? d__2 : d__3;
4718 if (workl[ibd + j] <= *tol * tempbnd) {
4722 } else if (!strncmp(which,"LA",2)) {
4723 if (workl[irz + j] >= thres1) {
4725 d__3 = fabs(workl[irz + j]);
4726 tempbnd = (d__2>d__3) ? d__2 : d__3;
4727 if (workl[ibd + j] <= *tol * tempbnd) {
4731 } else if (!strncmp(which,"SA",2)) {
4732 if (workl[irz + j] <= thres1) {
4734 d__3 = fabs(workl[irz + j]);
4735 tempbnd = (d__2>d__3) ? d__2 : d__3;
4736 if (workl[ibd + j] <= *tol * tempbnd) {
4740 } else if (!strncmp(which,"BE",2)) {
4741 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2) {
4743 d__3 = fabs(workl[irz + j]);
4744 tempbnd = (d__2>d__3) ? d__2 : d__3;
4745 if (workl[ibd + j] <= *tol * tempbnd) {
4750 if (j + 1 > nconv) {
4751 reord = select[j + 1] || reord;
4753 if (select[j + 1]) {
4759 F77_FUNC(scopy,SCOPY)(&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
4760 F77_FUNC(scopy,SCOPY)(ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
4762 F77_FUNC(ssteqr,SSTEQR)("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
4781 if (select[leftptr]) {
4785 } else if (! select[rghtptr]) {
4791 temp = workl[ihd + leftptr - 1];
4792 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
4793 workl[ihd + rghtptr - 1] = temp;
4794 F77_FUNC(scopy,SCOPY)(ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
4796 F77_FUNC(scopy,SCOPY)(ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
4797 iq + *ncv * (leftptr - 1)], &c__1);
4798 F77_FUNC(scopy,SCOPY)(ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
4805 if (leftptr < rghtptr) {
4813 F77_FUNC(scopy,SCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
4817 F77_FUNC(scopy,SCOPY)(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
4818 F77_FUNC(scopy,SCOPY)(ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
4821 if (!strncmp(type__, "REGULR",6)) {
4824 F77_FUNC(ssesrt,SSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
4826 F77_FUNC(scopy,SCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
4831 F77_FUNC(scopy,SCOPY)(ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
4832 if (!strncmp(type__, "SHIFTI", 6)) {
4834 for (k = 1; k <= i__1; ++k) {
4835 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
4837 } else if (!strncmp(type__, "BUCKLE",6)) {
4839 for (k = 1; k <= i__1; ++k) {
4840 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
4843 } else if (!strncmp(type__, "CAYLEY",6)) {
4845 for (k = 1; k <= i__1; ++k) {
4846 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
4847 workl[ihd + k - 1] - 1.);
4851 F77_FUNC(scopy,SCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
4852 F77_FUNC(ssortr,SSORTR)("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
4854 F77_FUNC(ssesrt,SSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
4856 F77_FUNC(scopy,SCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
4857 d__1 = bnorm2 / rnorm;
4858 F77_FUNC(sscal,SSCAL)(ncv, &d__1, &workl[ihb], &c__1);
4859 F77_FUNC(ssortr,SSORTR)("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
4864 if (*rvec && *howmny == 'A') {
4866 F77_FUNC(sgeqr2,SGEQR2)(ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
4869 F77_FUNC(sorm2r,SORM2R)("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
4870 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
4871 F77_FUNC(slacpy,SLACPY)("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
4874 for (j = 1; j <= i__1; ++j) {
4875 workl[ihb + j - 1] = 0.;
4877 workl[ihb + *ncv - 1] = 1.;
4878 F77_FUNC(sorm2r,SORM2R)("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
4879 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
4881 } else if (*rvec && *howmny == 'S') {
4885 if (!strncmp(type__, "REGULR",6) && *rvec) {
4888 for (j = 1; j <= i__1; ++j) {
4889 workl[ihb + j - 1] = rnorm * fabs(workl[ihb + j - 1]);
4892 } else if (strncmp(type__, "REGULR",6) && *rvec) {
4894 F77_FUNC(sscal,SSCAL)(ncv, &bnorm2, &workl[ihb], &c__1);
4895 if (!strncmp(type__, "SHIFTI",6)) {
4898 for (k = 1; k <= i__1; ++k) {
4899 d__2 = workl[iw + k - 1];
4900 workl[ihb + k - 1] = fabs(workl[ihb + k - 1])/(d__2 * d__2);
4903 } else if (!strncmp(type__, "BUCKLE",6)) {
4906 for (k = 1; k <= i__1; ++k) {
4907 d__2 = workl[iw + k - 1] - 1.;
4908 workl[ihb + k - 1] = *sigma * fabs(workl[ihb + k - 1])/(d__2 * d__2);
4911 } else if (!strncmp(type__, "CAYLEY",6)) {
4914 for (k = 1; k <= i__1; ++k) {
4915 workl[ihb + k - 1] = fabs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
4923 if (*rvec && (!strncmp(type__, "SHIFTI",6) || !strncmp(type__, "CAYLEY",6))) {
4926 for (k = 0; k <= i__1; ++k) {
4927 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
4930 } else if (*rvec && !strncmp(type__, "BUCKLE", 6)) {
4933 for (k = 0; k <= i__1; ++k) {
4934 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
4940 if (strncmp(type__, "REGULR",6)) {
4941 F77_FUNC(sger,SGER)(n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[