2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2004 David van der Spoel, Erik Lindahl, University of Groningen.
5 * Copyright (c) 2012,2013,2014,2015,2016,2018, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 #include "gmx_arpack.h"
43 #include "gromacs/utility/basedefinitions.h"
44 #include "gromacs/utility/real.h"
47 #include "gmx_lapack.h"
50 F77_FUNC(dstqrb, DSTQRB) (int * n,
66 int l1, ii, mm, lm1, mm1, nm1;
70 int lend, jtot, lendm1, lendp1, iscale;
72 int lendsv, nmaxit, icompz;
73 double ssfmax, ssfmin, safmin, minval, safmax, anorm;
100 minval = GMX_DOUBLE_MIN;
101 safmin = minval / GMX_DOUBLE_EPS;
102 safmax = 1. / safmin;
103 ssfmax = std::sqrt(safmax) / 3.;
104 ssfmin = std::sqrt(safmin) / eps2;
107 for (j = 1; j <= i__1; ++j)
132 for (m = l1; m <= i__1; ++m)
134 tst = std::abs(e[m]);
139 if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m+1])) * eps)
160 anorm = F77_FUNC(dlanst, DLANST) ("i", &i__1, &d__[l], &e[l]);
170 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
173 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
176 else if (anorm < ssfmin)
180 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
183 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
187 if (std::abs(d__[lend]) < std::abs(d__[l]))
201 for (m = l; m <= i__1; ++m)
203 d__2 = std::abs(e[m]);
205 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
229 F77_FUNC(dlaev2, DLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
231 work[*n - 1 + l] = s;
234 z__[l + 1] = c__ * tst - s * z__[l];
235 z__[l] = s * tst + c__ * z__[l];
239 F77_FUNC(dlae2, DLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
258 g = (d__[l + 1] - p) / (e[l] * 2.);
259 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
260 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
268 for (i__ = mm1; i__ >= i__1; --i__)
272 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
277 g = d__[i__ + 1] - p;
278 r__ = (d__[i__] - g) * s + c__ * 2. * b;
280 d__[i__ + 1] = g + p;
286 work[*n - 1 + i__] = -s;
295 F77_FUNC(dlasr, DLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
322 for (m = l; m >= i__1; --m)
324 d__2 = std::abs(e[m - 1]);
326 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m- 1]) + safmin)
350 F77_FUNC(dlaev2, DLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
354 z__[l] = c__ * tst - s * z__[l - 1];
355 z__[l - 1] = s * tst + c__ * z__[l - 1];
359 F77_FUNC(dlae2, DLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
379 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
380 r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
381 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
389 for (i__ = m; i__ <= i__1; ++i__)
393 F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
399 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
407 work[*n - 1 + i__] = s;
416 F77_FUNC(dlasr, DLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
439 i__1 = lendsv - lsv + 1;
440 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
443 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
446 else if (iscale == 2)
448 i__1 = lendsv - lsv + 1;
449 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
452 F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
461 for (i__ = 1; i__ <= i__1; ++i__)
474 F77_FUNC(dlasrt, DLASRT) ("i", n, &d__[1], info);
481 for (ii = 2; ii <= i__1; ++ii)
487 for (j = ii; j <= i__2; ++j)
513 F77_FUNC(dgetv0, DGETV0) (int * ido,
515 int gmx_unused * itry,
532 int v_dim1, v_offset, i__1;
540 v_offset = 1 + v_dim1;
556 F77_FUNC(dlarnv, DLARNV) (&idist, &iwork[1], n, &resid[1]);
563 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
582 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
588 else if (*bmat == 'I')
590 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
599 workd[*n * 3 + 4] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
600 workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
602 else if (*bmat == 'I')
604 workd[*n * 3 + 4] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
606 *rnorm = workd[*n * 3 + 4];
616 F77_FUNC(dgemv, DGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
617 &workd[*n + 1], &c__1);
619 F77_FUNC(dgemv, DGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
620 c_b22, &resid[1], &c__1);
624 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
630 else if (*bmat == 'I')
632 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
639 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
640 *rnorm = std::sqrt(std::abs(*rnorm));
642 else if (*bmat == 'I')
644 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
647 if (*rnorm > workd[*n * 3 + 4] * .717f)
656 workd[*n * 3 + 4] = *rnorm;
663 for (jj = 1; jj <= i__1; ++jj)
684 F77_FUNC(dsapps, DSAPPS) (int * n,
701 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
705 double r__, s, a1, a2, a3, a4;
716 v_offset = 1 + v_dim1;
719 h_offset = 1 + h_dim1;
722 q_offset = 1 + q_dim1;
725 epsmch = GMX_DOUBLE_EPS;
731 F77_FUNC(dlaset, DLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
739 for (jj = 1; jj <= i__1; ++jj)
747 for (i__ = istart; i__ <= i__2; ++i__)
749 big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__ + 1 + (h_dim1*2)]);
750 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
752 h__[i__ + 1 + h_dim1] = 0.;
763 f = h__[istart + (h_dim1 << 1)] - shift[jj];
764 g = h__[istart + 1 + h_dim1];
765 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
767 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
769 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
771 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
773 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
775 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
776 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
777 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
780 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
781 for (j = 1; j <= i__2; ++j)
783 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
785 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
786 c__ * q[j + (istart + 1) * q_dim1];
787 q[j + istart * q_dim1] = a1;
792 for (i__ = istart + 1; i__ <= i__2; ++i__)
795 f = h__[i__ + h_dim1];
796 g = s * h__[i__ + 1 + h_dim1];
798 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
799 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
808 h__[i__ + h_dim1] = r__;
810 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
812 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
814 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
816 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
819 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
820 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
821 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
824 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
825 for (j = 1; j <= i__3; ++j)
827 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
829 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
830 c__ * q[j + (i__ + 1) * q_dim1];
831 q[j + i__ * q_dim1] = a1;
840 if (h__[iend + h_dim1] < 0.)
842 h__[iend + h_dim1] = -h__[iend + h_dim1];
843 F77_FUNC(dscal, DSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
852 for (i__ = itop; i__ <= i__2; ++i__)
854 if (h__[i__ + 1 + h_dim1] > 0.)
866 for (i__ = itop; i__ <= i__1; ++i__)
868 big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__+ 1 + (h_dim1*2)]);
869 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
871 h__[i__ + 1 + h_dim1] = 0.;
876 if (h__[*kev + 1 + h_dim1] > 0.)
878 F77_FUNC(dgemv, DGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
879 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
883 for (i__ = 1; i__ <= i__1; ++i__)
885 i__2 = kplusp - i__ + 1;
886 F77_FUNC(dgemv, DGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
887 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
888 F77_FUNC(dcopy, DCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
893 F77_FUNC(dlacpy, DLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
895 if (h__[*kev + 1 + h_dim1] > 0.)
897 F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
900 F77_FUNC(dscal, DSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
901 if (h__[*kev + 1 + h_dim1] > 0.)
903 F77_FUNC(daxpy, DAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
917 F77_FUNC(dsortr, DSORTR) (const char * which,
932 if (!std::strncmp(which, "SA", 2))
941 for (i__ = igap; i__ <= i__1; ++i__)
951 if (x1[j] < x1[j + igap])
954 x1[j] = x1[j + igap];
959 x2[j] = x2[j + igap];
976 else if (!std::strncmp(which, "SM", 2))
985 for (i__ = igap; i__ <= i__1; ++i__)
995 if (std::abs(x1[j]) < std::abs(x1[j + igap]))
998 x1[j] = x1[j + igap];
1003 x2[j] = x2[j + igap];
1004 x2[j + igap] = temp;
1020 else if (!std::strncmp(which, "LA", 2))
1029 for (i__ = igap; i__ <= i__1; ++i__)
1039 if (x1[j] > x1[j + igap])
1042 x1[j] = x1[j + igap];
1043 x1[j + igap] = temp;
1047 x2[j] = x2[j + igap];
1048 x2[j + igap] = temp;
1064 else if (!std::strncmp(which, "LM", 2))
1074 for (i__ = igap; i__ <= i__1; ++i__)
1084 if (std::abs(x1[j]) > std::abs(x1[j + igap]))
1087 x1[j] = x1[j + igap];
1088 x1[j + igap] = temp;
1092 x2[j] = x2[j + igap];
1093 x2[j + igap] = temp;
1118 F77_FUNC(dsesrt, DSESRT) (const char * which,
1126 int a_dim1, a_offset, i__1;
1133 a_offset = 1 + a_dim1 * 0;
1138 if (!std::strncmp(which, "SA", 2))
1147 for (i__ = igap; i__ <= i__1; ++i__)
1157 if (x[j] < x[j + igap])
1164 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1165 a_dim1 + 1], &c__1);
1181 else if (!std::strncmp(which, "SM", 2))
1190 for (i__ = igap; i__ <= i__1; ++i__)
1200 if (std::abs(x[j]) < std::abs(x[j + igap]))
1207 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1208 a_dim1 + 1], &c__1);
1224 else if (!std::strncmp(which, "LA", 2))
1233 for (i__ = igap; i__ <= i__1; ++i__)
1243 if (x[j] > x[j + igap])
1250 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1251 a_dim1 + 1], &c__1);
1267 else if (!std::strncmp(which, "LM", 2))
1276 for (i__ = igap; i__ <= i__1; ++i__)
1286 if (std::abs(x[j]) > std::abs(x[j + igap]))
1293 F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1294 a_dim1 + 1], &c__1);
1319 F77_FUNC(dsgets, DSGETS) (int * ishift,
1335 if (!std::strncmp(which, "BE", 2))
1338 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1342 i__1 = (kevd2 < *np) ? kevd2 : *np;
1343 i__2 = (kevd2 > *np) ? kevd2 : *np;
1344 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[1], &c__1,
1345 &ritz[i__2 + 1], &c__1);
1346 i__1 = (kevd2 < *np) ? kevd2 : *np;
1347 i__2 = (kevd2 > *np) ? kevd2 : *np;
1348 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[1], &c__1,
1349 &bounds[i__2 + 1], &c__1);
1356 F77_FUNC(dsortr, DSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
1359 if (*ishift == 1 && *np > 0)
1362 F77_FUNC(dsortr, DSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
1363 F77_FUNC(dcopy, DCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
1373 F77_FUNC(dsconv, DSCONV) (int * n,
1389 eps23 = GMX_DOUBLE_EPS;
1390 eps23 = std::pow(eps23, c_b3);
1394 for (i__ = 1; i__ <= i__1; ++i__)
1398 d__3 = std::abs(ritz[i__]);
1399 temp = (d__2 > d__3) ? d__2 : d__3;
1400 if (bounds[i__] <= *tol * temp)
1411 F77_FUNC(dseigt, DSEIGT) (double * rnorm,
1421 int h_dim1, h_offset, i__1;
1430 h_offset = 1 + h_dim1;
1433 F77_FUNC(dcopy, DCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1435 F77_FUNC(dcopy, DCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1436 F77_FUNC(dstqrb, DSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1443 for (k = 1; k <= i__1; ++k)
1445 bounds[k] = *rnorm * std::abs(bounds[k]);
1458 F77_FUNC(dsaitr, DSAITR) (int * ido,
1482 int h_dim1, h_offset, v_dim1, v_offset, i__1;
1486 double safmin, minval;
1492 v_offset = 1 + v_dim1;
1495 h_offset = 1 + h_dim1;
1499 minval = GMX_DOUBLE_MIN;
1500 safmin = minval / GMX_DOUBLE_EPS;
1514 iwork[9] = iwork[8] + *n;
1515 iwork[10] = iwork[9] + *n;
1553 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
1554 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
1567 *info = iwork[12] - 1;
1574 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1575 if (*rnorm >= safmin)
1577 temp1 = 1. / *rnorm;
1578 F77_FUNC(dscal, DSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1579 F77_FUNC(dscal, DSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
1584 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
1585 v_dim1 + 1], n, &infol);
1586 F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
1591 F77_FUNC(dcopy, DCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1592 ipntr[1] = iwork[10];
1593 ipntr[2] = iwork[9];
1594 ipntr[3] = iwork[8];
1603 F77_FUNC(dcopy, DCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1612 ipntr[1] = iwork[9];
1613 ipntr[2] = iwork[8];
1618 else if (*bmat == 'I')
1620 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1630 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
1632 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1634 else if (*bmat == 'G')
1636 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1638 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1640 else if (*bmat == 'I')
1642 workd[*n * 3 + 3] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1647 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
1648 &c__1, &c_b42, &workd[iwork[9]], &c__1);
1652 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
1653 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
1656 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1657 c__1, &c_b18, &resid[1], &c__1);
1659 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1660 if (iwork[12] == 1 || iwork[4] == 1)
1662 h__[iwork[12] + h_dim1] = 0.;
1666 h__[iwork[12] + h_dim1] = *rnorm;
1674 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1675 ipntr[1] = iwork[9];
1676 ipntr[2] = iwork[8];
1681 else if (*bmat == 'I')
1683 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1691 *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1692 *rnorm = std::sqrt(std::abs(*rnorm));
1694 else if (*bmat == 'I')
1696 *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1699 if (*rnorm > workd[*n * 3 + 3] * .717f)
1706 F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
1707 c__1, &c_b42, &workd[iwork[9]], &c__1);
1709 F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1710 c__1, &c_b18, &resid[1], &c__1);
1712 if (iwork[12] == 1 || iwork[4] == 1)
1714 h__[iwork[12] + h_dim1] = 0.;
1716 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1721 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1722 ipntr[1] = iwork[9];
1723 ipntr[2] = iwork[8];
1728 else if (*bmat == 'I')
1730 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1737 workd[*n * 3 + 2] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1739 workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
1741 else if (*bmat == 'I')
1743 workd[*n * 3 + 2] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1747 if (workd[*n * 3 + 2] > *rnorm * .717f)
1750 *rnorm = workd[*n * 3 + 2];
1756 *rnorm = workd[*n * 3 + 2];
1764 for (jj = 1; jj <= i__1; ++jj)
1776 if (h__[iwork[12] + h_dim1] < 0.)
1778 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1779 if (iwork[12] < *k + *np)
1781 F77_FUNC(dscal, DSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1785 F77_FUNC(dscal, DSCAL) (n, &c_b50, &resid[1], &c__1);
1790 if (iwork[12] > *k + *np)
1810 F77_FUNC(dsaup2, DSAUP2) (int * ido,
1819 int gmx_unused * iupd,
1840 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
1859 v_offset = 1 + v_dim1;
1862 h_offset = 1 + h_dim1;
1865 q_offset = 1 + q_dim1;
1869 eps23 = GMX_DOUBLE_EPS;
1870 eps23 = std::pow(eps23, c_b3);
1883 iwork[7] = iwork[9] + iwork[10];
1906 F77_FUNC(dgetv0, DGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
1907 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
1915 if (workd[*n * 3 + 1] == 0.)
1940 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
1941 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1967 F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1],
1968 &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1986 F77_FUNC(dseigt, DSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
1987 bounds[1], &workl[1], &ierr);
1995 F77_FUNC(dcopy, DCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1996 F77_FUNC(dcopy, DCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
2000 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2002 F77_FUNC(dcopy, DCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
2003 F77_FUNC(dsconv, DSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
2007 for (j = 1; j <= i__1; ++j)
2009 if (bounds[j] == 0.)
2016 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
2019 if (!std::strncmp(which, "BE", 2))
2022 std::strncpy(wprime, "SA", 2);
2023 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2025 nevm2 = *nev - nevd2;
2028 i__1 = (nevd2 < *np) ? nevd2 : *np;
2029 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
2030 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
2031 &ritz[((i__2 > i__3) ? i__2 : i__3)],
2033 i__1 = (nevd2 < *np) ? nevd2 : *np;
2034 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
2035 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
2036 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
2044 if (!std::strncmp(which, "LM", 2))
2046 std::strncpy(wprime, "SM", 2);
2048 if (!std::strncmp(which, "SM", 2))
2050 std::strncpy(wprime, "LM", 2);
2052 if (!std::strncmp(which, "LA", 2))
2054 std::strncpy(wprime, "SA", 2);
2056 if (!std::strncmp(which, "SA", 2))
2058 std::strncpy(wprime, "LA", 2);
2061 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2066 for (j = 1; j <= i__1; ++j)
2069 d__3 = std::abs(ritz[j]);
2070 temp = (d__2 > d__3) ? d__2 : d__3;
2074 std::strncpy(wprime, "LA", 2);
2075 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
2078 for (j = 1; j <= i__1; ++j)
2081 d__3 = std::abs(ritz[j]);
2082 temp = (d__2 > d__3) ? d__2 : d__3;
2086 if (!std::strncmp(which, "BE", 2))
2089 std::strncpy(wprime, "LA", 2);
2090 F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2095 F77_FUNC(dsortr, DSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2099 h__[h_dim1 + 1] = workd[*n * 3 + 1];
2102 if (iwork[6] > *mxiter && iwork[8] < *nev)
2106 if (*np == 0 && iwork[8] < iwork[9])
2115 else if (iwork[8] < *nev && *ishift == 1)
2118 i__1 = iwork[8], i__2 = *np / 2;
2119 *nev += (i__1 < i__2) ? i__1 : i__2;
2120 if (*nev == 1 && iwork[7] >= 6)
2122 *nev = iwork[7] / 2;
2124 else if (*nev == 1 && iwork[7] > 2)
2128 *np = iwork[7] - *nev;
2133 F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2153 F77_FUNC(dcopy, DCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
2156 F77_FUNC(dsapps, DSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
2157 resid[1], &q[q_offset], ldq, &workd[1]);
2162 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
2169 else if (*bmat == 'I')
2171 F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
2178 workd[*n * 3 + 1] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
2179 workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
2181 else if (*bmat == 'I')
2183 workd[*n * 3 + 1] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
2205 F77_FUNC(dsaupd, DSAUPD) (int * ido,
2223 int v_dim1, v_offset, i__1, i__2;
2229 v_offset = 1 + v_dim1;
2241 iwork[5] = iparam[1];
2242 iwork[10] = iparam[3];
2243 iwork[12] = iparam[4];
2246 iwork[11] = iparam[7];
2257 else if (*ncv <= *nev || *ncv > *n)
2263 iwork[15] = *ncv - *nev;
2269 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
2270 std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
2271 std::strncmp(which, "BE", 2))
2275 if (*bmat != 'I' && *bmat != 'G')
2281 if (*lworkl < i__1 * i__1 + (*ncv << 3))
2285 if (iwork[11] < 1 || iwork[11] > 5)
2289 else if (iwork[11] == 1 && *bmat == 'G')
2293 else if (iwork[5] < 0 || iwork[5] > 1)
2297 else if (*nev == 1 && !std::strncmp(which, "BE", 2))
2315 *tol = GMX_DOUBLE_EPS;
2318 iwork[15] = *ncv - *nev;
2321 i__1 = i__2 * i__2 + (*ncv << 3);
2322 for (j = 1; j <= i__1; ++j)
2330 iwork[16] = iwork[3] + (iwork[8] << 1);
2331 iwork[1] = iwork[16] + *ncv;
2332 iwork[4] = iwork[1] + *ncv;
2334 iwork[7] = iwork[4] + i__1 * i__1;
2335 iwork[14] = iwork[7] + *ncv * 3;
2337 ipntr[4] = iwork[14];
2338 ipntr[5] = iwork[3];
2339 ipntr[6] = iwork[16];
2340 ipntr[7] = iwork[1];
2341 ipntr[11] = iwork[7];
2344 F77_FUNC(dsaup2, DSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
2345 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
2346 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
2347 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
2352 iparam[8] = iwork[15];
2359 iparam[3] = iwork[10];
2360 iparam[5] = iwork[15];
2380 F77_FUNC(dseupd, DSEUPD) (int * rvec,
2381 const char * howmny,
2403 double c_b21 = 2/3.;
2406 int v_dim1, v_offset, z_dim1, z_offset, i__1;
2407 double d__1, d__2, d__3;
2409 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2421 double thres1 = 0, thres2 = 0;
2425 int leftptr, rghtptr;
2431 z_offset = 1 + z_dim1;
2436 v_offset = 1 + v_dim1;
2464 if (*ncv <= *nev || *ncv > *n)
2468 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
2469 std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
2470 std::strncmp(which, "BE", 2))
2474 if (*bmat != 'I' && *bmat != 'G')
2478 if (*howmny != 'A' && *howmny != 'P' &&
2479 *howmny != 'S' && *rvec)
2483 if (*rvec && *howmny == 'S')
2488 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
2493 if (mode == 1 || mode == 2)
2495 std::strncpy(type__, "REGULR", 6);
2499 std::strncpy(type__, "SHIFTI", 6);
2503 std::strncpy(type__, "BUCKLE", 6);
2507 std::strncpy(type__, "CAYLEY", 6);
2513 if (mode == 1 && *bmat == 'G')
2517 if (*nev == 1 && !std::strncmp(which, "BE", 2))
2536 iw = iq + ldh * *ncv;
2537 next = iw + (*ncv << 1);
2543 irz = ipntr[11] + *ncv;
2547 eps23 = GMX_DOUBLE_EPS;
2548 eps23 = std::pow(eps23, c_b21);
2555 else if (*bmat == 'G')
2557 bnorm2 = F77_FUNC(dnrm2, DNRM2) (n, &workd[1], &c__1);
2563 if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2) ||
2564 !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
2568 else if (!std::strncmp(which, "BE", 2))
2572 ism = (*nev > nconv) ? *nev : nconv;
2575 thres1 = workl[ism];
2576 thres2 = workl[ilg];
2584 for (j = 0; j <= i__1; ++j)
2587 if (!std::strncmp(which, "LM", 2))
2589 if (std::abs(workl[irz + j]) >= std::abs(thres1))
2592 d__3 = std::abs(workl[irz + j]);
2593 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2594 if (workl[ibd + j] <= *tol * tempbnd)
2600 else if (!std::strncmp(which, "SM", 2))
2602 if (std::abs(workl[irz + j]) <= std::abs(thres1))
2605 d__3 = std::abs(workl[irz + j]);
2606 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2607 if (workl[ibd + j] <= *tol * tempbnd)
2613 else if (!std::strncmp(which, "LA", 2))
2615 if (workl[irz + j] >= thres1)
2618 d__3 = std::abs(workl[irz + j]);
2619 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2620 if (workl[ibd + j] <= *tol * tempbnd)
2626 else if (!std::strncmp(which, "SA", 2))
2628 if (workl[irz + j] <= thres1)
2631 d__3 = std::abs(workl[irz + j]);
2632 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2633 if (workl[ibd + j] <= *tol * tempbnd)
2639 else if (!std::strncmp(which, "BE", 2))
2641 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
2644 d__3 = std::abs(workl[irz + j]);
2645 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2646 if (workl[ibd + j] <= *tol * tempbnd)
2654 reord = select[j + 1] || reord;
2663 F77_FUNC(dcopy, DCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2664 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2666 F77_FUNC(dsteqr, DSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
2688 if (select[leftptr])
2694 else if (!select[rghtptr])
2703 temp = workl[ihd + leftptr - 1];
2704 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2705 workl[ihd + rghtptr - 1] = temp;
2706 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
2708 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
2709 iq + *ncv * (leftptr - 1)], &c__1);
2710 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
2717 if (leftptr < rghtptr)
2726 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2732 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2733 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2736 if (!std::strncmp(type__, "REGULR", 6))
2741 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2745 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2752 F77_FUNC(dcopy, DCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2753 if (!std::strncmp(type__, "SHIFTI", 6))
2756 for (k = 1; k <= i__1; ++k)
2758 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2761 else if (!std::strncmp(type__, "BUCKLE", 6))
2764 for (k = 1; k <= i__1; ++k)
2766 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
2770 else if (!std::strncmp(type__, "CAYLEY", 6))
2773 for (k = 1; k <= i__1; ++k)
2775 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
2776 workl[ihd + k - 1] - 1.);
2780 F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2781 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2784 F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2788 F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2789 d__1 = bnorm2 / rnorm;
2790 F77_FUNC(dscal, DSCAL) (ncv, &d__1, &workl[ihb], &c__1);
2791 F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2796 if (*rvec && *howmny == 'A')
2799 F77_FUNC(dgeqr2, DGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2802 F77_FUNC(dorm2r, DORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
2803 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
2804 F77_FUNC(dlacpy, DLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2807 for (j = 1; j <= i__1; ++j)
2809 workl[ihb + j - 1] = 0.;
2811 workl[ihb + *ncv - 1] = 1.;
2812 F77_FUNC(dorm2r, DORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
2813 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
2816 else if (*rvec && *howmny == 'S')
2821 if (!std::strncmp(type__, "REGULR", 6) && *rvec)
2825 for (j = 1; j <= i__1; ++j)
2827 workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
2831 else if (std::strncmp(type__, "REGULR", 6) && *rvec)
2834 F77_FUNC(dscal, DSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
2835 if (!std::strncmp(type__, "SHIFTI", 6))
2839 for (k = 1; k <= i__1; ++k)
2841 d__2 = workl[iw + k - 1];
2842 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1])/(d__2 * d__2);
2846 else if (!std::strncmp(type__, "BUCKLE", 6))
2850 for (k = 1; k <= i__1; ++k)
2852 d__2 = workl[iw + k - 1] - 1.;
2853 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1])/(d__2 * d__2);
2857 else if (!std::strncmp(type__, "CAYLEY", 6))
2861 for (k = 1; k <= i__1; ++k)
2863 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2871 if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
2875 for (k = 0; k <= i__1; ++k)
2877 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2881 else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
2885 for (k = 0; k <= i__1; ++k)
2887 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
2893 if (std::strncmp(type__, "REGULR", 6))
2895 F77_FUNC(dger, DGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
2909 /* Selected single precision arpack routines */
2913 F77_FUNC(sstqrb, SSTQRB) (int * n,
2927 int i__, j, k, l, m;
2929 int l1, ii, mm, lm1, mm1, nm1;
2930 float rt1, rt2, eps;
2933 int lend, jtot, lendm1, lendp1, iscale;
2935 int lendsv, nmaxit, icompz;
2936 float ssfmax, ssfmin, safmin, minval, safmax, anorm;
2959 eps = GMX_FLOAT_EPS;
2963 minval = GMX_FLOAT_MIN;
2964 safmin = minval / GMX_FLOAT_EPS;
2965 safmax = 1. / safmin;
2966 ssfmax = std::sqrt(safmax) / 3.;
2967 ssfmin = std::sqrt(safmin) / eps2;
2970 for (j = 1; j <= i__1; ++j)
2995 for (m = l1; m <= i__1; ++m)
2997 tst = std::abs(e[m]);
3002 if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m+1])) * eps)
3022 i__1 = lend - l + 1;
3023 anorm = F77_FUNC(slanst, SLANST) ("i", &i__1, &d__[l], &e[l]);
3032 i__1 = lend - l + 1;
3033 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
3036 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
3039 else if (anorm < ssfmin)
3042 i__1 = lend - l + 1;
3043 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
3046 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
3050 if (std::abs(d__[lend]) < std::abs(d__[l]))
3064 for (m = l; m <= i__1; ++m)
3066 d__2 = std::abs(e[m]);
3068 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
3092 F77_FUNC(slaev2, SLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
3094 work[*n - 1 + l] = s;
3097 z__[l + 1] = c__ * tst - s * z__[l];
3098 z__[l] = s * tst + c__ * z__[l];
3102 F77_FUNC(slae2, SLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
3121 g = (d__[l + 1] - p) / (e[l] * 2.);
3122 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3123 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
3131 for (i__ = mm1; i__ >= i__1; --i__)
3135 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3140 g = d__[i__ + 1] - p;
3141 r__ = (d__[i__] - g) * s + c__ * 2. * b;
3143 d__[i__ + 1] = g + p;
3149 work[*n - 1 + i__] = -s;
3158 F77_FUNC(slasr, SLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
3185 for (m = l; m >= i__1; --m)
3187 d__2 = std::abs(e[m - 1]);
3189 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m- 1]) + safmin)
3213 F77_FUNC(slaev2, SLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
3217 z__[l] = c__ * tst - s * z__[l - 1];
3218 z__[l - 1] = s * tst + c__ * z__[l - 1];
3222 F77_FUNC(slae2, SLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
3242 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
3243 r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3244 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
3252 for (i__ = m; i__ <= i__1; ++i__)
3256 F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3262 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
3270 work[*n - 1 + i__] = s;
3279 F77_FUNC(slasr, SLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
3302 i__1 = lendsv - lsv + 1;
3303 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
3305 i__1 = lendsv - lsv;
3306 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
3309 else if (iscale == 2)
3311 i__1 = lendsv - lsv + 1;
3312 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
3314 i__1 = lendsv - lsv;
3315 F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
3324 for (i__ = 1; i__ <= i__1; ++i__)
3337 F77_FUNC(slasrt, SLASRT) ("i", n, &d__[1], info);
3344 for (ii = 2; ii <= i__1; ++ii)
3350 for (j = ii; j <= i__2; ++j)
3376 F77_FUNC(sgetv0, SGETV0) (int * ido,
3378 int gmx_unused * itry,
3395 int v_dim1, v_offset, i__1;
3403 v_offset = 1 + v_dim1;
3419 F77_FUNC(slarnv, SLARNV) (&idist, &iwork[1], n, &resid[1]);
3426 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3445 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
3451 else if (*bmat == 'I')
3453 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3462 workd[*n * 3 + 4] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3463 workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
3465 else if (*bmat == 'I')
3467 workd[*n * 3 + 4] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3469 *rnorm = workd[*n * 3 + 4];
3479 F77_FUNC(sgemv, SGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
3480 &workd[*n + 1], &c__1);
3482 F77_FUNC(sgemv, SGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
3483 c_b22, &resid[1], &c__1);
3487 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3493 else if (*bmat == 'I')
3495 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3502 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3503 *rnorm = std::sqrt(std::abs(*rnorm));
3505 else if (*bmat == 'I')
3507 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3510 if (*rnorm > workd[*n * 3 + 4] * .717f)
3519 workd[*n * 3 + 4] = *rnorm;
3526 for (jj = 1; jj <= i__1; ++jj)
3547 F77_FUNC(ssapps, SSAPPS) (int * n,
3564 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
3568 float r__, s, a1, a2, a3, a4;
3579 v_offset = 1 + v_dim1;
3582 h_offset = 1 + h_dim1;
3585 q_offset = 1 + q_dim1;
3588 epsmch = GMX_FLOAT_EPS;
3592 kplusp = *kev + *np;
3594 F77_FUNC(slaset, SLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3602 for (jj = 1; jj <= i__1; ++jj)
3610 for (i__ = istart; i__ <= i__2; ++i__)
3612 big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__ + 1 + (h_dim1*2)]);
3613 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3615 h__[i__ + 1 + h_dim1] = 0.;
3626 f = h__[istart + (h_dim1 << 1)] - shift[jj];
3627 g = h__[istart + 1 + h_dim1];
3628 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3630 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
3632 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
3634 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
3636 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
3638 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
3639 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3640 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
3643 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
3644 for (j = 1; j <= i__2; ++j)
3646 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
3648 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
3649 c__ * q[j + (istart + 1) * q_dim1];
3650 q[j + istart * q_dim1] = a1;
3655 for (i__ = istart + 1; i__ <= i__2; ++i__)
3658 f = h__[i__ + h_dim1];
3659 g = s * h__[i__ + 1 + h_dim1];
3661 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3662 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3671 h__[i__ + h_dim1] = r__;
3673 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
3675 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
3677 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
3679 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
3682 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
3683 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3684 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
3687 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
3688 for (j = 1; j <= i__3; ++j)
3690 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
3692 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
3693 c__ * q[j + (i__ + 1) * q_dim1];
3694 q[j + i__ * q_dim1] = a1;
3703 if (h__[iend + h_dim1] < 0.)
3705 h__[iend + h_dim1] = -h__[iend + h_dim1];
3706 F77_FUNC(sscal, SSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3715 for (i__ = itop; i__ <= i__2; ++i__)
3717 if (h__[i__ + 1 + h_dim1] > 0.)
3729 for (i__ = itop; i__ <= i__1; ++i__)
3731 big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__+ 1 + (h_dim1*2)]);
3732 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3734 h__[i__ + 1 + h_dim1] = 0.;
3739 if (h__[*kev + 1 + h_dim1] > 0.)
3741 F77_FUNC(sgemv, SGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
3742 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
3746 for (i__ = 1; i__ <= i__1; ++i__)
3748 i__2 = kplusp - i__ + 1;
3749 F77_FUNC(sgemv, SGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
3750 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
3751 F77_FUNC(scopy, SCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
3756 F77_FUNC(slacpy, SLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3758 if (h__[*kev + 1 + h_dim1] > 0.)
3760 F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3763 F77_FUNC(sscal, SSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3764 if (h__[*kev + 1 + h_dim1] > 0.)
3766 F77_FUNC(saxpy, SAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
3780 F77_FUNC(ssortr, SSORTR) (const char * which,
3795 if (!std::strncmp(which, "SA", 2))
3804 for (i__ = igap; i__ <= i__1; ++i__)
3814 if (x1[j] < x1[j + igap])
3817 x1[j] = x1[j + igap];
3818 x1[j + igap] = temp;
3822 x2[j] = x2[j + igap];
3823 x2[j + igap] = temp;
3839 else if (!std::strncmp(which, "SM", 2))
3848 for (i__ = igap; i__ <= i__1; ++i__)
3858 if (std::abs(x1[j]) < std::abs(x1[j + igap]))
3861 x1[j] = x1[j + igap];
3862 x1[j + igap] = temp;
3866 x2[j] = x2[j + igap];
3867 x2[j + igap] = temp;
3883 else if (!std::strncmp(which, "LA", 2))
3892 for (i__ = igap; i__ <= i__1; ++i__)
3902 if (x1[j] > x1[j + igap])
3905 x1[j] = x1[j + igap];
3906 x1[j + igap] = temp;
3910 x2[j] = x2[j + igap];
3911 x2[j + igap] = temp;
3927 else if (!std::strncmp(which, "LM", 2))
3937 for (i__ = igap; i__ <= i__1; ++i__)
3947 if (std::abs(x1[j]) > std::abs(x1[j + igap]))
3950 x1[j] = x1[j + igap];
3951 x1[j + igap] = temp;
3955 x2[j] = x2[j + igap];
3956 x2[j + igap] = temp;
3981 F77_FUNC(ssesrt, SSESRT) (const char * which,
3989 int a_dim1, a_offset, i__1;
3996 a_offset = 1 + a_dim1 * 0;
4001 if (!std::strncmp(which, "SA", 2))
4010 for (i__ = igap; i__ <= i__1; ++i__)
4020 if (x[j] < x[j + igap])
4027 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4028 a_dim1 + 1], &c__1);
4044 else if (!std::strncmp(which, "SM", 2))
4053 for (i__ = igap; i__ <= i__1; ++i__)
4063 if (std::abs(x[j]) < std::abs(x[j + igap]))
4070 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4071 a_dim1 + 1], &c__1);
4087 else if (!std::strncmp(which, "LA", 2))
4096 for (i__ = igap; i__ <= i__1; ++i__)
4106 if (x[j] > x[j + igap])
4113 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4114 a_dim1 + 1], &c__1);
4130 else if (!std::strncmp(which, "LM", 2))
4139 for (i__ = igap; i__ <= i__1; ++i__)
4149 if (std::abs(x[j]) > std::abs(x[j + igap]))
4156 F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4157 a_dim1 + 1], &c__1);
4182 F77_FUNC(ssgets, SSGETS) (int * ishift,
4198 if (!std::strncmp(which, "BE", 2))
4201 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
4205 i__1 = (kevd2 < *np) ? kevd2 : *np;
4206 i__2 = (kevd2 > *np) ? kevd2 : *np;
4207 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[1], &c__1,
4208 &ritz[i__2 + 1], &c__1);
4209 i__1 = (kevd2 < *np) ? kevd2 : *np;
4210 i__2 = (kevd2 > *np) ? kevd2 : *np;
4211 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[1], &c__1,
4212 &bounds[i__2 + 1], &c__1);
4219 F77_FUNC(ssortr, SSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
4222 if (*ishift == 1 && *np > 0)
4225 F77_FUNC(ssortr, SSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
4226 F77_FUNC(scopy, SCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
4236 F77_FUNC(ssconv, SSCONV) (int * n,
4252 eps23 = GMX_FLOAT_EPS;
4253 eps23 = std::pow(eps23, c_b3);
4257 for (i__ = 1; i__ <= i__1; ++i__)
4261 d__3 = std::abs(ritz[i__]);
4262 temp = (d__2 > d__3) ? d__2 : d__3;
4263 if (bounds[i__] <= *tol * temp)
4274 F77_FUNC(sseigt, SSEIGT) (float * rnorm,
4284 int h_dim1, h_offset, i__1;
4293 h_offset = 1 + h_dim1;
4296 F77_FUNC(scopy, SCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
4298 F77_FUNC(scopy, SCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
4299 F77_FUNC(sstqrb, SSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
4306 for (k = 1; k <= i__1; ++k)
4308 bounds[k] = *rnorm * std::abs(bounds[k]);
4321 F77_FUNC(ssaitr, SSAITR) (int * ido,
4345 int h_dim1, h_offset, v_dim1, v_offset, i__1;
4349 float safmin, minval;
4355 v_offset = 1 + v_dim1;
4358 h_offset = 1 + h_dim1;
4362 minval = GMX_FLOAT_MIN;
4363 safmin = minval / GMX_FLOAT_EPS;
4377 iwork[9] = iwork[8] + *n;
4378 iwork[10] = iwork[9] + *n;
4416 F77_FUNC(sgetv0, sgetv0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
4417 &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
4430 *info = iwork[12] - 1;
4437 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
4438 if (*rnorm >= safmin)
4440 temp1 = 1. / *rnorm;
4441 F77_FUNC(sscal, SSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
4442 F77_FUNC(sscal, SSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
4447 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
4448 v_dim1 + 1], n, &infol);
4449 F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
4454 F77_FUNC(scopy, SCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
4455 ipntr[1] = iwork[10];
4456 ipntr[2] = iwork[9];
4457 ipntr[3] = iwork[8];
4466 F77_FUNC(scopy, SCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
4475 ipntr[1] = iwork[9];
4476 ipntr[2] = iwork[8];
4481 else if (*bmat == 'I')
4483 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4493 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
4495 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4497 else if (*bmat == 'G')
4499 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4501 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4503 else if (*bmat == 'I')
4505 workd[*n * 3 + 3] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4510 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
4511 &c__1, &c_b42, &workd[iwork[9]], &c__1);
4515 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
4516 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
4519 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4520 c__1, &c_b18, &resid[1], &c__1);
4522 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
4523 if (iwork[12] == 1 || iwork[4] == 1)
4525 h__[iwork[12] + h_dim1] = 0.;
4529 h__[iwork[12] + h_dim1] = *rnorm;
4537 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4538 ipntr[1] = iwork[9];
4539 ipntr[2] = iwork[8];
4544 else if (*bmat == 'I')
4546 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4554 *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4555 *rnorm = std::sqrt(std::abs(*rnorm));
4557 else if (*bmat == 'I')
4559 *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4562 if (*rnorm > workd[*n * 3 + 3] * .717f)
4569 F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
4570 c__1, &c_b42, &workd[iwork[9]], &c__1);
4572 F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4573 c__1, &c_b18, &resid[1], &c__1);
4575 if (iwork[12] == 1 || iwork[4] == 1)
4577 h__[iwork[12] + h_dim1] = 0.;
4579 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
4584 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4585 ipntr[1] = iwork[9];
4586 ipntr[2] = iwork[8];
4591 else if (*bmat == 'I')
4593 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4600 workd[*n * 3 + 2] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4602 workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
4604 else if (*bmat == 'I')
4606 workd[*n * 3 + 2] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4610 if (workd[*n * 3 + 2] > *rnorm * .717f)
4613 *rnorm = workd[*n * 3 + 2];
4619 *rnorm = workd[*n * 3 + 2];
4627 for (jj = 1; jj <= i__1; ++jj)
4639 if (h__[iwork[12] + h_dim1] < 0.)
4641 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4642 if (iwork[12] < *k + *np)
4644 F77_FUNC(sscal, SSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4648 F77_FUNC(sscal, SSCAL) (n, &c_b50, &resid[1], &c__1);
4653 if (iwork[12] > *k + *np)
4673 F77_FUNC(ssaup2, SSAUP2) (int * ido,
4682 int gmx_unused * iupd,
4703 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
4722 v_offset = 1 + v_dim1;
4725 h_offset = 1 + h_dim1;
4728 q_offset = 1 + q_dim1;
4732 eps23 = GMX_FLOAT_EPS;
4733 eps23 = std::pow(eps23, c_b3);
4746 iwork[7] = iwork[9] + iwork[10];
4769 F77_FUNC(sgetv0, SGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
4770 resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
4778 if (workd[*n * 3 + 1] == 0.)
4803 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
4804 1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
4830 F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[
4831 v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[
4849 F77_FUNC(sseigt, SSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
4850 bounds[1], &workl[1], &ierr);
4858 F77_FUNC(scopy, SCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4859 F77_FUNC(scopy, SCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4863 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4865 F77_FUNC(scopy, SCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4866 F77_FUNC(ssconv, SSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4871 for (j = 1; j <= i__1; ++j)
4873 if (bounds[j] == 0.)
4880 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
4883 if (!std::strncmp(which, "BE", 2))
4886 std::strncpy(wprime, "SA", 2);
4887 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4889 nevm2 = *nev - nevd2;
4892 i__1 = (nevd2 < *np) ? nevd2 : *np;
4893 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4894 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
4895 &ritz[((i__2 > i__3) ? i__2 : i__3)],
4897 i__1 = (nevd2 < *np) ? nevd2 : *np;
4898 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4899 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
4900 &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
4908 if (!std::strncmp(which, "LM", 2))
4910 std::strncpy(wprime, "SM", 2);
4912 if (!std::strncmp(which, "SM", 2))
4914 std::strncpy(wprime, "LM", 2);
4916 if (!std::strncmp(which, "LA", 2))
4918 std::strncpy(wprime, "SA", 2);
4920 if (!std::strncmp(which, "SA", 2))
4922 std::strncpy(wprime, "LA", 2);
4925 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4930 for (j = 1; j <= i__1; ++j)
4933 d__3 = std::abs(ritz[j]);
4934 temp = (d__2 > d__3) ? d__2 : d__3;
4938 std::strncpy(wprime, "LA", 2);
4939 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4942 for (j = 1; j <= i__1; ++j)
4945 d__3 = std::abs(ritz[j]);
4946 temp = (d__2 > d__3) ? d__2 : d__3;
4950 if (!std::strncmp(which, "BE", 2))
4953 std::strncpy(wprime, "LA", 2);
4954 F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4959 F77_FUNC(ssortr, SSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4963 h__[h_dim1 + 1] = workd[*n * 3 + 1];
4966 if (iwork[6] > *mxiter && iwork[8] < *nev)
4970 if (*np == 0 && iwork[8] < iwork[9])
4979 else if (iwork[8] < *nev && *ishift == 1)
4982 i__1 = iwork[8], i__2 = *np / 2;
4983 *nev += (i__1 < i__2) ? i__1 : i__2;
4984 if (*nev == 1 && iwork[7] >= 6)
4986 *nev = iwork[7] / 2;
4988 else if (*nev == 1 && iwork[7] > 2)
4992 *np = iwork[7] - *nev;
4997 F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
5017 F77_FUNC(scopy, SCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
5020 F77_FUNC(ssapps, SSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
5021 resid[1], &q[q_offset], ldq, &workd[1]);
5026 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
5033 else if (*bmat == 'I')
5035 F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
5042 workd[*n * 3 + 1] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
5043 workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
5045 else if (*bmat == 'I')
5047 workd[*n * 3 + 1] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
5069 F77_FUNC(ssaupd, SSAUPD) (int * ido,
5087 int v_dim1, v_offset, i__1, i__2;
5093 v_offset = 1 + v_dim1;
5105 iwork[5] = iparam[1];
5106 iwork[10] = iparam[3];
5107 iwork[12] = iparam[4];
5110 iwork[11] = iparam[7];
5121 else if (*ncv <= *nev || *ncv > *n)
5127 iwork[15] = *ncv - *nev;
5133 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
5134 std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
5135 std::strncmp(which, "BE", 2))
5139 if (*bmat != 'I' && *bmat != 'G')
5145 if (*lworkl < i__1 * i__1 + (*ncv << 3))
5149 if (iwork[11] < 1 || iwork[11] > 5)
5153 else if (iwork[11] == 1 && *bmat == 'G')
5157 else if (iwork[5] < 0 || iwork[5] > 1)
5161 else if (*nev == 1 && !std::strncmp(which, "BE", 2))
5179 *tol = GMX_FLOAT_EPS;
5182 iwork[15] = *ncv - *nev;
5185 i__1 = i__2 * i__2 + (*ncv << 3);
5186 for (j = 1; j <= i__1; ++j)
5194 iwork[16] = iwork[3] + (iwork[8] << 1);
5195 iwork[1] = iwork[16] + *ncv;
5196 iwork[4] = iwork[1] + *ncv;
5198 iwork[7] = iwork[4] + i__1 * i__1;
5199 iwork[14] = iwork[7] + *ncv * 3;
5201 ipntr[4] = iwork[14];
5202 ipntr[5] = iwork[3];
5203 ipntr[6] = iwork[16];
5204 ipntr[7] = iwork[1];
5205 ipntr[11] = iwork[7];
5208 F77_FUNC(ssaup2, SSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
5209 iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
5210 workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
5211 workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
5216 iparam[8] = iwork[15];
5223 iparam[3] = iwork[10];
5224 iparam[5] = iwork[15];
5244 F77_FUNC(sseupd, SSEUPD) (int * rvec,
5245 const char * howmny,
5270 int v_dim1, v_offset, z_dim1, z_offset, i__1;
5271 float d__1, d__2, d__3;
5273 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
5285 float thres1 = 0, thres2 = 0;
5289 int leftptr, rghtptr;
5295 z_offset = 1 + z_dim1;
5300 v_offset = 1 + v_dim1;
5328 if (*ncv <= *nev || *ncv > *n)
5332 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
5333 std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
5334 std::strncmp(which, "BE", 2))
5338 if (*bmat != 'I' && *bmat != 'G')
5342 if (*howmny != 'A' && *howmny != 'P' &&
5343 *howmny != 'S' && *rvec)
5347 if (*rvec && *howmny == 'S')
5352 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
5357 if (mode == 1 || mode == 2)
5359 std::strncpy(type__, "REGULR", 6);
5363 std::strncpy(type__, "SHIFTI", 6);
5367 std::strncpy(type__, "BUCKLE", 6);
5371 std::strncpy(type__, "CAYLEY", 6);
5377 if (mode == 1 && *bmat == 'G')
5381 if (*nev == 1 && !std::strncmp(which, "BE", 2))
5400 iw = iq + ldh * *ncv;
5401 next = iw + (*ncv << 1);
5407 irz = ipntr[11] + *ncv;
5411 eps23 = GMX_FLOAT_EPS;
5412 eps23 = std::pow(eps23, c_b21);
5419 else if (*bmat == 'G')
5421 bnorm2 = F77_FUNC(snrm2, SNRM2) (n, &workd[1], &c__1);
5427 if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2) ||
5428 !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
5432 else if (!std::strncmp(which, "BE", 2))
5436 ism = (*nev > nconv) ? *nev : nconv;
5439 thres1 = workl[ism];
5440 thres2 = workl[ilg];
5448 for (j = 0; j <= i__1; ++j)
5451 if (!std::strncmp(which, "LM", 2))
5453 if (std::abs(workl[irz + j]) >= std::abs(thres1))
5456 d__3 = std::abs(workl[irz + j]);
5457 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5458 if (workl[ibd + j] <= *tol * tempbnd)
5464 else if (!std::strncmp(which, "SM", 2))
5466 if (std::abs(workl[irz + j]) <= std::abs(thres1))
5469 d__3 = std::abs(workl[irz + j]);
5470 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5471 if (workl[ibd + j] <= *tol * tempbnd)
5477 else if (!std::strncmp(which, "LA", 2))
5479 if (workl[irz + j] >= thres1)
5482 d__3 = std::abs(workl[irz + j]);
5483 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5484 if (workl[ibd + j] <= *tol * tempbnd)
5490 else if (!std::strncmp(which, "SA", 2))
5492 if (workl[irz + j] <= thres1)
5495 d__3 = std::abs(workl[irz + j]);
5496 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5497 if (workl[ibd + j] <= *tol * tempbnd)
5503 else if (!std::strncmp(which, "BE", 2))
5505 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
5508 d__3 = std::abs(workl[irz + j]);
5509 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5510 if (workl[ibd + j] <= *tol * tempbnd)
5518 reord = select[j + 1] || reord;
5527 F77_FUNC(scopy, SCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
5528 F77_FUNC(scopy, SCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
5530 F77_FUNC(ssteqr, SSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
5552 if (select[leftptr])
5558 else if (!select[rghtptr])
5567 temp = workl[ihd + leftptr - 1];
5568 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
5569 workl[ihd + rghtptr - 1] = temp;
5570 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
5572 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
5573 iq + *ncv * (leftptr - 1)], &c__1);
5574 F77_FUNC(scopy, SCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
5581 if (leftptr < rghtptr)
5590 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5596 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
5597 F77_FUNC(scopy, SCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
5600 if (!std::strncmp(type__, "REGULR", 6))
5605 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5609 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5616 F77_FUNC(scopy, SCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
5617 if (!std::strncmp(type__, "SHIFTI", 6))
5620 for (k = 1; k <= i__1; ++k)
5622 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
5625 else if (!std::strncmp(type__, "BUCKLE", 6))
5628 for (k = 1; k <= i__1; ++k)
5630 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
5634 else if (!std::strncmp(type__, "CAYLEY", 6))
5637 for (k = 1; k <= i__1; ++k)
5639 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
5640 workl[ihd + k - 1] - 1.);
5644 F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5645 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
5648 F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5652 F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5653 d__1 = bnorm2 / rnorm;
5654 F77_FUNC(sscal, SSCAL) (ncv, &d__1, &workl[ihb], &c__1);
5655 F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
5660 if (*rvec && *howmny == 'A')
5663 F77_FUNC(sgeqr2, SGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
5666 F77_FUNC(sorm2r, SORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
5667 workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
5668 F77_FUNC(slacpy, SLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
5671 for (j = 1; j <= i__1; ++j)
5673 workl[ihb + j - 1] = 0.;
5675 workl[ihb + *ncv - 1] = 1.;
5676 F77_FUNC(sorm2r, SORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
5677 workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
5680 else if (*rvec && *howmny == 'S')
5685 if (!std::strncmp(type__, "REGULR", 6) && *rvec)
5689 for (j = 1; j <= i__1; ++j)
5691 workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
5695 else if (std::strncmp(type__, "REGULR", 6) && *rvec)
5698 F77_FUNC(sscal, SSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
5699 if (!std::strncmp(type__, "SHIFTI", 6))
5703 for (k = 1; k <= i__1; ++k)
5705 d__2 = workl[iw + k - 1];
5706 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1])/(d__2 * d__2);
5710 else if (!std::strncmp(type__, "BUCKLE", 6))
5714 for (k = 1; k <= i__1; ++k)
5716 d__2 = workl[iw + k - 1] - 1.;
5717 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1])/(d__2 * d__2);
5721 else if (!std::strncmp(type__, "CAYLEY", 6))
5725 for (k = 1; k <= i__1; ++k)
5727 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
5735 if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
5739 for (k = 0; k <= i__1; ++k)
5741 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
5745 else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
5749 for (k = 0; k <= i__1; ++k)
5751 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
5757 if (std::strncmp(type__, "REGULR", 6))
5759 F77_FUNC(sger, SGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[