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 by the GROMACS development team.
6 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gmx_arpack.h"
44 #include "gromacs/utility/basedefinitions.h"
45 #include "gromacs/utility/real.h"
48 #include "gmx_lapack.h"
52 static void F77_FUNC(dstqrb, DSTQRB)(int* n, double* d__, double* e, double* z__, double* work, int* info)
63 int l1, ii, mm, lm1, mm1, nm1;
67 int lend, jtot, lendm1, lendp1, iscale;
69 int lendsv, nmaxit, icompz;
70 double ssfmax, ssfmin, safmin, minval, safmax, anorm;
97 minval = GMX_DOUBLE_MIN;
98 safmin = minval / GMX_DOUBLE_EPS;
100 ssfmax = std::sqrt(safmax) / 3.;
101 ssfmin = std::sqrt(safmin) / eps2;
104 for (j = 1; j <= i__1; ++j)
128 for (m = l1; m <= i__1; ++m)
130 tst = std::abs(e[m]);
135 if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m + 1])) * eps)
156 anorm = F77_FUNC(dlanst, DLANST)("i", &i__1, &d__[l], &e[l]);
166 F77_FUNC(dlascl, DLASCL)
167 ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, info);
169 F77_FUNC(dlascl, DLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, info);
171 else if (anorm < ssfmin)
175 F77_FUNC(dlascl, DLASCL)
176 ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, info);
178 F77_FUNC(dlascl, DLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, info);
181 if (std::abs(d__[lend]) < std::abs(d__[l]))
195 for (m = l; m <= i__1; ++m)
197 d__2 = std::abs(e[m]);
199 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
223 F77_FUNC(dlaev2, DLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
225 work[*n - 1 + l] = s;
228 z__[l + 1] = c__ * tst - s * z__[l];
229 z__[l] = s * tst + c__ * z__[l];
233 F77_FUNC(dlae2, DLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
252 g = (d__[l + 1] - p) / (e[l] * 2.);
253 r__ = F77_FUNC(dlapy2, DLAPY2)(&g, &c_b31);
254 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__));
262 for (i__ = mm1; i__ >= i__1; --i__)
266 F77_FUNC(dlartg, DLARTG)(&g, &f, &c__, &s, &r__);
271 g = d__[i__ + 1] - p;
272 r__ = (d__[i__] - g) * s + c__ * 2. * b;
274 d__[i__ + 1] = g + p;
280 work[*n - 1 + i__] = -s;
288 F77_FUNC(dlasr, DLASR)
289 ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &z__[l], &c__1);
314 for (m = l; m >= i__1; --m)
316 d__2 = std::abs(e[m - 1]);
318 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m - 1]) + safmin)
342 F77_FUNC(dlaev2, DLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s);
345 z__[l] = c__ * tst - s * z__[l - 1];
346 z__[l - 1] = s * tst + c__ * z__[l - 1];
350 F77_FUNC(dlae2, DLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
370 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
371 r__ = F77_FUNC(dlapy2, DLAPY2)(&g, &c_b31);
372 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__));
380 for (i__ = m; i__ <= i__1; ++i__)
384 F77_FUNC(dlartg, DLARTG)(&g, &f, &c__, &s, &r__);
390 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
398 work[*n - 1 + i__] = s;
406 F77_FUNC(dlasr, DLASR)
407 ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &z__[m], &c__1);
428 i__1 = lendsv - lsv + 1;
429 F77_FUNC(dlascl, DLASCL)
430 ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], n, info);
432 F77_FUNC(dlascl, DLASCL)
433 ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n, info);
435 else if (iscale == 2)
437 i__1 = lendsv - lsv + 1;
438 F77_FUNC(dlascl, DLASCL)
439 ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], n, info);
441 F77_FUNC(dlascl, DLASCL)
442 ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n, info);
450 for (i__ = 1; i__ <= i__1; ++i__)
463 F77_FUNC(dlasrt, DLASRT)("i", n, &d__[1], info);
469 for (ii = 2; ii <= i__1; ++ii)
475 for (j = ii; j <= i__2; ++j)
499 static void F77_FUNC(dgetv0, DGETV0)(int* ido,
501 int gmx_unused* itry,
518 int v_dim1, v_offset, i__1;
526 v_offset = 1 + v_dim1;
542 F77_FUNC(dlarnv, DLARNV)(&idist, &iwork[1], n, &resid[1]);
549 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
568 F77_FUNC(dcopy, DCOPY)(n, &workd[*n + 1], &c__1, &resid[1], &c__1);
574 else if (*bmat == 'I')
576 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
585 workd[*n * 3 + 4] = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
586 workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
588 else if (*bmat == 'I')
590 workd[*n * 3 + 4] = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
592 *rnorm = workd[*n * 3 + 4];
602 F77_FUNC(dgemv, DGEMV)
603 ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24, &workd[*n + 1], &c__1);
605 F77_FUNC(dgemv, DGEMV)
606 ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &c_b22, &resid[1], &c__1);
610 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
616 else if (*bmat == 'I')
618 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
625 *rnorm = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
626 *rnorm = std::sqrt(std::abs(*rnorm));
628 else if (*bmat == 'I')
630 *rnorm = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
633 if (*rnorm > workd[*n * 3 + 4] * .717F)
642 workd[*n * 3 + 4] = *rnorm;
649 for (jj = 1; jj <= i__1; ++jj)
666 static void F77_FUNC(dsapps, DSAPPS)(int* n,
683 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;
686 double r__, s, a1, a2, a3, a4;
697 v_offset = 1 + v_dim1;
700 h_offset = 1 + h_dim1;
703 q_offset = 1 + q_dim1;
706 epsmch = GMX_DOUBLE_EPS;
712 F77_FUNC(dlaset, DLASET)("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
720 for (jj = 1; jj <= i__1; ++jj)
728 for (i__ = istart; i__ <= i__2; ++i__)
730 big = std::abs(h__[i__ + (h_dim1 * 2)]) + std::abs(h__[i__ + 1 + (h_dim1 * 2)]);
731 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
733 h__[i__ + 1 + h_dim1] = 0.;
744 f = h__[istart + (h_dim1 << 1)] - shift[jj];
745 g = h__[istart + 1 + h_dim1];
746 F77_FUNC(dlartg, DLARTG)(&f, &g, &c__, &s, &r__);
748 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 + h_dim1];
749 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (h_dim1 << 1)];
750 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 + h_dim1];
751 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 << 1)];
752 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
753 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
754 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
757 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
758 for (j = 1; j <= i__2; ++j)
760 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) * q_dim1];
761 q[j + (istart + 1) * q_dim1] =
762 -s * q[j + istart * q_dim1] + c__ * q[j + (istart + 1) * q_dim1];
763 q[j + istart * q_dim1] = a1;
767 for (i__ = istart + 1; i__ <= i__2; ++i__)
770 f = h__[i__ + h_dim1];
771 g = s * h__[i__ + 1 + h_dim1];
773 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
774 F77_FUNC(dlartg, DLARTG)(&f, &g, &c__, &s, &r__);
783 h__[i__ + h_dim1] = r__;
785 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 + h_dim1];
786 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1 << 1)];
787 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)];
788 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 + h_dim1];
790 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
791 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
792 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
795 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
796 for (j = 1; j <= i__3; ++j)
798 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) * q_dim1];
799 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] + c__ * q[j + (i__ + 1) * q_dim1];
800 q[j + i__ * q_dim1] = a1;
807 if (h__[iend + h_dim1] < 0.)
809 h__[iend + h_dim1] = -h__[iend + h_dim1];
810 F77_FUNC(dscal, DSCAL)(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
819 for (i__ = itop; i__ <= i__2; ++i__)
821 if (h__[i__ + 1 + h_dim1] > 0.)
832 for (i__ = itop; i__ <= i__1; ++i__)
834 big = std::abs(h__[i__ + (h_dim1 * 2)]) + std::abs(h__[i__ + 1 + (h_dim1 * 2)]);
835 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
837 h__[i__ + 1 + h_dim1] = 0.;
841 if (h__[*kev + 1 + h_dim1] > 0.)
843 F77_FUNC(dgemv, DGEMV)
844 ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) * q_dim1 + 1], &c__1, &c_b4,
845 &workd[*n + 1], &c__1);
849 for (i__ = 1; i__ <= i__1; ++i__)
851 i__2 = kplusp - i__ + 1;
852 F77_FUNC(dgemv, DGEMV)
853 ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) * q_dim1 + 1], &c__1, &c_b4,
855 F77_FUNC(dcopy, DCOPY)(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &c__1);
858 F77_FUNC(dlacpy, DLACPY)("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
860 if (h__[*kev + 1 + h_dim1] > 0.)
862 F77_FUNC(dcopy, DCOPY)(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
865 F77_FUNC(dscal, DSCAL)(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
866 if (h__[*kev + 1 + h_dim1] > 0.)
868 F77_FUNC(daxpy, DAXPY)
869 (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1, &resid[1], &c__1);
878 static void F77_FUNC(dsortr, DSORTR)(const char* which, int* apply, int* n, double* x1, double* x2)
888 if (!std::strncmp(which, "SA", 2))
897 for (i__ = igap; i__ <= i__1; ++i__)
907 if (x1[j] < x1[j + igap])
910 x1[j] = x1[j + igap];
915 x2[j] = x2[j + igap];
930 else if (!std::strncmp(which, "SM", 2))
939 for (i__ = igap; i__ <= i__1; ++i__)
949 if (std::abs(x1[j]) < std::abs(x1[j + igap]))
952 x1[j] = x1[j + igap];
957 x2[j] = x2[j + igap];
972 else if (!std::strncmp(which, "LA", 2))
981 for (i__ = igap; i__ <= i__1; ++i__)
991 if (x1[j] > x1[j + igap])
994 x1[j] = x1[j + igap];
999 x2[j] = x2[j + igap];
1000 x2[j + igap] = temp;
1014 else if (!std::strncmp(which, "LM", 2))
1024 for (i__ = igap; i__ <= i__1; ++i__)
1034 if (std::abs(x1[j]) > std::abs(x1[j + igap]))
1037 x1[j] = x1[j + igap];
1038 x1[j + igap] = temp;
1042 x2[j] = x2[j + igap];
1043 x2[j + igap] = temp;
1063 static void F77_FUNC(dsesrt,
1064 DSESRT)(const char* which, int* apply, int* n, double* x, int* na, double* a, int* lda)
1066 int a_dim1, a_offset, i__1;
1073 a_offset = 1 + a_dim1 * 0;
1078 if (!std::strncmp(which, "SA", 2))
1087 for (i__ = igap; i__ <= i__1; ++i__)
1097 if (x[j] < x[j + igap])
1104 F77_FUNC(dswap, DSWAP)
1105 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
1119 else if (!std::strncmp(which, "SM", 2))
1128 for (i__ = igap; i__ <= i__1; ++i__)
1138 if (std::abs(x[j]) < std::abs(x[j + igap]))
1145 F77_FUNC(dswap, DSWAP)
1146 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
1160 else if (!std::strncmp(which, "LA", 2))
1169 for (i__ = igap; i__ <= i__1; ++i__)
1179 if (x[j] > x[j + igap])
1186 F77_FUNC(dswap, DSWAP)
1187 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
1201 else if (!std::strncmp(which, "LM", 2))
1210 for (i__ = igap; i__ <= i__1; ++i__)
1220 if (std::abs(x[j]) > std::abs(x[j + igap]))
1227 F77_FUNC(dswap, DSWAP)
1228 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
1248 static void F77_FUNC(dsgets,
1249 DSGETS)(int* ishift, const char* which, int* kev, int* np, double* ritz, double* bounds, double* shifts)
1259 if (!std::strncmp(which, "BE", 2))
1262 F77_FUNC(dsortr, DSORTR)("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1266 i__1 = (kevd2 < *np) ? kevd2 : *np;
1267 i__2 = (kevd2 > *np) ? kevd2 : *np;
1268 F77_FUNC(dswap, DSWAP)(&i__1, &ritz[1], &c__1, &ritz[i__2 + 1], &c__1);
1269 i__1 = (kevd2 < *np) ? kevd2 : *np;
1270 i__2 = (kevd2 > *np) ? kevd2 : *np;
1271 F77_FUNC(dswap, DSWAP)(&i__1, &bounds[1], &c__1, &bounds[i__2 + 1], &c__1);
1277 F77_FUNC(dsortr, DSORTR)(which, &c__1, &i__1, &ritz[1], &bounds[1]);
1280 if (*ishift == 1 && *np > 0)
1283 F77_FUNC(dsortr, DSORTR)("SM", &c__1, np, &bounds[1], &ritz[1]);
1284 F77_FUNC(dcopy, DCOPY)(np, &ritz[1], &c__1, &shifts[1], &c__1);
1292 static void F77_FUNC(dsconv, DSCONV)(int* n, double* ritz, double* bounds, double* tol, int* nconv)
1294 double c_b3 = 2 / 3.;
1304 eps23 = GMX_DOUBLE_EPS;
1305 eps23 = std::pow(eps23, c_b3);
1309 for (i__ = 1; i__ <= i__1; ++i__)
1313 d__3 = std::abs(ritz[i__]);
1314 temp = (d__2 > d__3) ? d__2 : d__3;
1315 if (bounds[i__] <= *tol * temp)
1325 static void F77_FUNC(dseigt, DSEIGT)(double* rnorm,
1335 int h_dim1, h_offset, i__1;
1344 h_offset = 1 + h_dim1;
1347 F77_FUNC(dcopy, DCOPY)(n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1349 F77_FUNC(dcopy, DCOPY)(&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1350 F77_FUNC(dstqrb, DSTQRB)(n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1357 for (k = 1; k <= i__1; ++k)
1359 bounds[k] = *rnorm * std::abs(bounds[k]);
1368 static void F77_FUNC(dsaitr, DSAITR)(int* ido,
1392 int h_dim1, h_offset, v_dim1, v_offset, i__1;
1396 double safmin, minval;
1402 v_offset = 1 + v_dim1;
1405 h_offset = 1 + h_dim1;
1409 minval = GMX_DOUBLE_MIN;
1410 safmin = minval / GMX_DOUBLE_EPS;
1424 iwork[9] = iwork[8] + *n;
1425 iwork[10] = iwork[9] + *n;
1463 F77_FUNC(dgetv0, DGETV0)
1464 (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv, &resid[1], rnorm, &ipntr[1],
1465 &workd[1], &iwork[21], &iwork[7]);
1478 *info = iwork[12] - 1;
1485 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1486 if (*rnorm >= safmin)
1488 temp1 = 1. / *rnorm;
1489 F77_FUNC(dscal, DSCAL)(n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1490 F77_FUNC(dscal, DSCAL)(n, &temp1, &workd[iwork[8]], &c__1);
1495 F77_FUNC(dlascl, DLASCL)
1496 ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] * v_dim1 + 1], n, &infol);
1497 F77_FUNC(dlascl, DLASCL)
1498 ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[8]], n, &infol);
1502 F77_FUNC(dcopy, DCOPY)(n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1503 ipntr[1] = iwork[10];
1504 ipntr[2] = iwork[9];
1505 ipntr[3] = iwork[8];
1514 F77_FUNC(dcopy, DCOPY)(n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1523 ipntr[1] = iwork[9];
1524 ipntr[2] = iwork[8];
1529 else if (*bmat == 'I')
1531 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1541 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[iwork[10]], &c__1);
1542 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1544 else if (*bmat == 'G')
1546 workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1547 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1549 else if (*bmat == 'I')
1551 workd[*n * 3 + 3] = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
1556 F77_FUNC(dgemv, DGEMV)
1557 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &c__1, &c_b42,
1558 &workd[iwork[9]], &c__1);
1562 F77_FUNC(dgemv, DGEMV)
1563 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]], &c__1, &c_b42,
1564 &workd[iwork[9]], &c__1);
1567 F77_FUNC(dgemv, DGEMV)
1568 ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &c__1, &c_b18, &resid[1], &c__1);
1570 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1571 if (iwork[12] == 1 || iwork[4] == 1)
1573 h__[iwork[12] + h_dim1] = 0.;
1577 h__[iwork[12] + h_dim1] = *rnorm;
1585 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1586 ipntr[1] = iwork[9];
1587 ipntr[2] = iwork[8];
1592 else if (*bmat == 'I')
1594 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1602 *rnorm = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1603 *rnorm = std::sqrt(std::abs(*rnorm));
1605 else if (*bmat == 'I')
1607 *rnorm = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
1610 if (*rnorm > workd[*n * 3 + 3] * .717F)
1617 F77_FUNC(dgemv, DGEMV)
1618 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &c__1, &c_b42,
1619 &workd[iwork[9]], &c__1);
1621 F77_FUNC(dgemv, DGEMV)
1622 ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &c__1, &c_b18, &resid[1], &c__1);
1624 if (iwork[12] == 1 || iwork[4] == 1)
1626 h__[iwork[12] + h_dim1] = 0.;
1628 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1633 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1634 ipntr[1] = iwork[9];
1635 ipntr[2] = iwork[8];
1640 else if (*bmat == 'I')
1642 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1649 workd[*n * 3 + 2] = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1650 workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
1652 else if (*bmat == 'I')
1654 workd[*n * 3 + 2] = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
1658 if (workd[*n * 3 + 2] > *rnorm * .717F)
1661 *rnorm = workd[*n * 3 + 2];
1666 *rnorm = workd[*n * 3 + 2];
1674 for (jj = 1; jj <= i__1; ++jj)
1686 if (h__[iwork[12] + h_dim1] < 0.)
1688 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1689 if (iwork[12] < *k + *np)
1691 F77_FUNC(dscal, DSCAL)(n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1695 F77_FUNC(dscal, DSCAL)(n, &c_b50, &resid[1], &c__1);
1700 if (iwork[12] > *k + *np)
1715 static void F77_FUNC(dsaup2, DSAUP2)(int* ido,
1724 int gmx_unused* iupd,
1741 double c_b3 = 2 / 3.;
1745 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2, i__3;
1763 v_offset = 1 + v_dim1;
1766 h_offset = 1 + h_dim1;
1769 q_offset = 1 + q_dim1;
1773 eps23 = GMX_DOUBLE_EPS;
1774 eps23 = std::pow(eps23, c_b3);
1787 iwork[7] = iwork[9] + iwork[10];
1810 F77_FUNC(dgetv0, DGETV0)
1811 (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &resid[1], &workd[*n * 3 + 1],
1812 &ipntr[1], &workd[1], &iwork[41], info);
1819 if (workd[*n * 3 + 1] == 0.)
1844 F77_FUNC(dsaitr, DSAITR)
1845 (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 + 1], &v[v_offset], ldv,
1846 &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[21], info);
1871 F77_FUNC(dsaitr, DSAITR)
1872 (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[v_offset], ldv, &h__[h_offset],
1873 ldh, &ipntr[1], &workd[1], &iwork[21], info);
1890 F77_FUNC(dseigt, DSEIGT)
1891 (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &bounds[1], &workl[1], &ierr);
1899 F77_FUNC(dcopy, DCOPY)(&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1900 F77_FUNC(dcopy, DCOPY)(&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
1904 F77_FUNC(dsgets, DSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
1906 F77_FUNC(dcopy, DCOPY)(nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
1907 F77_FUNC(dsconv, DSCONV)(nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
1911 for (j = 1; j <= i__1; ++j)
1913 if (bounds[j] == 0.)
1920 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
1923 if (!std::strncmp(which, "BE", 2))
1926 std::strncpy(wprime, "SA", 2);
1927 F77_FUNC(dsortr, DSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
1929 nevm2 = *nev - nevd2;
1932 i__1 = (nevd2 < *np) ? nevd2 : *np;
1933 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
1934 F77_FUNC(dswap, DSWAP)
1935 (&i__1, &ritz[nevm2 + 1], &c__1, &ritz[((i__2 > i__3) ? i__2 : i__3)], &c__1);
1936 i__1 = (nevd2 < *np) ? nevd2 : *np;
1937 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
1938 F77_FUNC(dswap, DSWAP)
1939 (&i__1, &bounds[nevm2 + 1], &c__1, &bounds[((i__2 > i__3) ? i__2 : i__3) + 1], &c__1);
1945 if (!std::strncmp(which, "LM", 2))
1947 std::strncpy(wprime, "SM", 2);
1949 if (!std::strncmp(which, "SM", 2))
1951 std::strncpy(wprime, "LM", 2);
1953 if (!std::strncmp(which, "LA", 2))
1955 std::strncpy(wprime, "SA", 2);
1957 if (!std::strncmp(which, "SA", 2))
1959 std::strncpy(wprime, "LA", 2);
1962 F77_FUNC(dsortr, DSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
1966 for (j = 1; j <= i__1; ++j)
1969 d__3 = std::abs(ritz[j]);
1970 temp = (d__2 > d__3) ? d__2 : d__3;
1974 std::strncpy(wprime, "LA", 2);
1975 F77_FUNC(dsortr, DSORTR)(wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
1978 for (j = 1; j <= i__1; ++j)
1981 d__3 = std::abs(ritz[j]);
1982 temp = (d__2 > d__3) ? d__2 : d__3;
1986 if (!std::strncmp(which, "BE", 2))
1989 std::strncpy(wprime, "LA", 2);
1990 F77_FUNC(dsortr, DSORTR)(wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
1994 F77_FUNC(dsortr, DSORTR)(which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
1997 h__[h_dim1 + 1] = workd[*n * 3 + 1];
2000 if (iwork[6] > *mxiter && iwork[8] < *nev)
2004 if (*np == 0 && iwork[8] < iwork[9])
2012 else if (iwork[8] < *nev && *ishift == 1)
2015 i__1 = iwork[8], i__2 = *np / 2;
2016 *nev += (i__1 < i__2) ? i__1 : i__2;
2017 if (*nev == 1 && iwork[7] >= 6)
2019 *nev = iwork[7] / 2;
2021 else if (*nev == 1 && iwork[7] > 2)
2025 *np = iwork[7] - *nev;
2030 F77_FUNC(dsgets, DSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2049 F77_FUNC(dcopy, DCOPY)(np, &workl[1], &c__1, &ritz[1], &c__1);
2052 F77_FUNC(dsapps, DSAPPS)
2053 (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &resid[1], &q[q_offset], ldq,
2059 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
2066 else if (*bmat == 'I')
2068 F77_FUNC(dcopy, DCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
2075 workd[*n * 3 + 1] = F77_FUNC(ddot, DDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
2076 workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
2078 else if (*bmat == 'I')
2080 workd[*n * 3 + 1] = F77_FUNC(dnrm2, DNRM2)(n, &resid[1], &c__1);
2099 void F77_FUNC(dsaupd, DSAUPD)(int* ido,
2117 int v_dim1, v_offset, i__1, i__2;
2123 v_offset = 1 + v_dim1;
2135 iwork[5] = iparam[1];
2136 iwork[10] = iparam[3];
2137 iwork[12] = iparam[4];
2140 iwork[11] = iparam[7];
2151 else if (*ncv <= *nev || *ncv > *n)
2157 iwork[15] = *ncv - *nev;
2163 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) && std::strncmp(which, "LA", 2)
2164 && std::strncmp(which, "SA", 2) && std::strncmp(which, "BE", 2))
2168 if (*bmat != 'I' && *bmat != 'G')
2174 if (*lworkl < i__1 * i__1 + (*ncv << 3))
2178 if (iwork[11] < 1 || iwork[11] > 5)
2182 else if (iwork[11] == 1 && *bmat == 'G')
2186 else if (iwork[5] < 0 || iwork[5] > 1)
2190 else if (*nev == 1 && !std::strncmp(which, "BE", 2))
2208 *tol = GMX_DOUBLE_EPS;
2211 iwork[15] = *ncv - *nev;
2214 i__1 = i__2 * i__2 + (*ncv << 3);
2215 for (j = 1; j <= i__1; ++j)
2223 iwork[16] = iwork[3] + (iwork[8] << 1);
2224 iwork[1] = iwork[16] + *ncv;
2225 iwork[4] = iwork[1] + *ncv;
2227 iwork[7] = iwork[4] + i__1 * i__1;
2228 iwork[14] = iwork[7] + *ncv * 3;
2230 ipntr[4] = iwork[14];
2231 ipntr[5] = iwork[3];
2232 ipntr[6] = iwork[16];
2233 ipntr[7] = iwork[1];
2234 ipntr[11] = iwork[7];
2237 F77_FUNC(dsaup2, DSAUP2)
2238 (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &iwork[11], &iwork[6], &iwork[5],
2239 &iwork[10], &v[v_offset], ldv, &workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]],
2240 &workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1], &iwork[21], info);
2244 iparam[8] = iwork[15];
2251 iparam[3] = iwork[10];
2252 iparam[5] = iwork[15];
2269 void F77_FUNC(dseupd, DSEUPD)(int* rvec,
2292 double c_b21 = 2 / 3.;
2295 int v_dim1, v_offset, z_dim1, z_offset, i__1;
2296 double d__1, d__2, d__3;
2298 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2310 double thres1 = 0, thres2 = 0;
2314 int leftptr, rghtptr;
2320 z_offset = 1 + z_dim1;
2325 v_offset = 1 + v_dim1;
2353 if (*ncv <= *nev || *ncv > *n)
2357 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) && std::strncmp(which, "LA", 2)
2358 && std::strncmp(which, "SA", 2) && std::strncmp(which, "BE", 2))
2362 if (*bmat != 'I' && *bmat != 'G')
2366 if (*howmny != 'A' && *howmny != 'P' && *howmny != 'S' && *rvec)
2370 if (*rvec && *howmny == 'S')
2375 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
2380 if (mode == 1 || mode == 2)
2382 std::strncpy(type__, "REGULR", 6);
2386 std::strncpy(type__, "SHIFTI", 6);
2390 std::strncpy(type__, "BUCKLE", 6);
2394 std::strncpy(type__, "CAYLEY", 6);
2400 if (mode == 1 && *bmat == 'G')
2404 if (*nev == 1 && !std::strncmp(which, "BE", 2))
2423 iw = iq + ldh * *ncv;
2424 next = iw + (*ncv << 1);
2430 irz = ipntr[11] + *ncv;
2434 eps23 = GMX_DOUBLE_EPS;
2435 eps23 = std::pow(eps23, c_b21);
2442 else if (*bmat == 'G')
2444 bnorm2 = F77_FUNC(dnrm2, DNRM2)(n, &workd[1], &c__1);
2450 if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2)
2451 || !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
2453 else if (!std::strncmp(which, "BE", 2))
2457 ism = (*nev > nconv) ? *nev : nconv;
2460 thres1 = workl[ism];
2461 thres2 = workl[ilg];
2467 for (j = 0; j <= i__1; ++j)
2470 if (!std::strncmp(which, "LM", 2))
2472 if (std::abs(workl[irz + j]) >= std::abs(thres1))
2475 d__3 = std::abs(workl[irz + j]);
2476 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2477 if (workl[ibd + j] <= *tol * tempbnd)
2483 else if (!std::strncmp(which, "SM", 2))
2485 if (std::abs(workl[irz + j]) <= std::abs(thres1))
2488 d__3 = std::abs(workl[irz + j]);
2489 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2490 if (workl[ibd + j] <= *tol * tempbnd)
2496 else if (!std::strncmp(which, "LA", 2))
2498 if (workl[irz + j] >= thres1)
2501 d__3 = std::abs(workl[irz + j]);
2502 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2503 if (workl[ibd + j] <= *tol * tempbnd)
2509 else if (!std::strncmp(which, "SA", 2))
2511 if (workl[irz + j] <= thres1)
2514 d__3 = std::abs(workl[irz + j]);
2515 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2516 if (workl[ibd + j] <= *tol * tempbnd)
2522 else if (!std::strncmp(which, "BE", 2))
2524 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
2527 d__3 = std::abs(workl[irz + j]);
2528 tempbnd = (d__2 > d__3) ? d__2 : d__3;
2529 if (workl[ibd + j] <= *tol * tempbnd)
2537 reord = select[j + 1] || reord;
2546 F77_FUNC(dcopy, DCOPY)(&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2547 F77_FUNC(dcopy, DCOPY)(ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2549 F77_FUNC(dsteqr, DSTEQR)
2550 ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &workl[iw], &ierr);
2571 if (select[leftptr])
2576 else if (!select[rghtptr])
2584 temp = workl[ihd + leftptr - 1];
2585 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2586 workl[ihd + rghtptr - 1] = temp;
2587 F77_FUNC(dcopy, DCOPY)
2588 (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[iw], &c__1);
2589 F77_FUNC(dcopy, DCOPY)
2590 (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[iq + *ncv * (leftptr - 1)], &c__1);
2591 F77_FUNC(dcopy, DCOPY)
2592 (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr - 1)], &c__1);
2597 if (leftptr < rghtptr)
2605 F77_FUNC(dcopy, DCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2610 F77_FUNC(dcopy, DCOPY)(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2611 F77_FUNC(dcopy, DCOPY)(ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2613 if (!std::strncmp(type__, "REGULR", 6))
2618 F77_FUNC(dsesrt, DSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2622 F77_FUNC(dcopy, DCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2628 F77_FUNC(dcopy, DCOPY)(ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2629 if (!std::strncmp(type__, "SHIFTI", 6))
2632 for (k = 1; k <= i__1; ++k)
2634 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2637 else if (!std::strncmp(type__, "BUCKLE", 6))
2640 for (k = 1; k <= i__1; ++k)
2642 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd + k - 1] - 1.);
2645 else if (!std::strncmp(type__, "CAYLEY", 6))
2648 for (k = 1; k <= i__1; ++k)
2650 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (workl[ihd + k - 1] - 1.);
2654 F77_FUNC(dcopy, DCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2655 F77_FUNC(dsortr, DSORTR)("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2658 F77_FUNC(dsesrt, DSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2662 F77_FUNC(dcopy, DCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2663 d__1 = bnorm2 / rnorm;
2664 F77_FUNC(dscal, DSCAL)(ncv, &d__1, &workl[ihb], &c__1);
2665 F77_FUNC(dsortr, DSORTR)("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2669 if (*rvec && *howmny == 'A')
2672 F77_FUNC(dgeqr2, DGEQR2)
2673 (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb], &ierr);
2675 F77_FUNC(dorm2r, DORM2R)
2676 ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &v[v_offset],
2677 ldv, &workd[*n + 1], &ierr);
2678 F77_FUNC(dlacpy, DLACPY)("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2681 for (j = 1; j <= i__1; ++j)
2683 workl[ihb + j - 1] = 0.;
2685 workl[ihb + *ncv - 1] = 1.;
2686 F77_FUNC(dorm2r, DORM2R)
2687 ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2690 else if (*rvec && *howmny == 'S')
2694 if (!std::strncmp(type__, "REGULR", 6) && *rvec)
2698 for (j = 1; j <= i__1; ++j)
2700 workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
2703 else if (std::strncmp(type__, "REGULR", 6) && *rvec)
2706 F77_FUNC(dscal, DSCAL)(ncv, &bnorm2, &workl[ihb], &c__1);
2707 if (!std::strncmp(type__, "SHIFTI", 6))
2711 for (k = 1; k <= i__1; ++k)
2713 d__2 = workl[iw + k - 1];
2714 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1]) / (d__2 * d__2);
2717 else if (!std::strncmp(type__, "BUCKLE", 6))
2721 for (k = 1; k <= i__1; ++k)
2723 d__2 = workl[iw + k - 1] - 1.;
2724 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1]) / (d__2 * d__2);
2727 else if (!std::strncmp(type__, "CAYLEY", 6))
2731 for (k = 1; k <= i__1; ++k)
2733 workl[ihb + k - 1] =
2734 std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2739 if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
2743 for (k = 0; k <= i__1; ++k)
2745 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2748 else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
2752 for (k = 0; k <= i__1; ++k)
2754 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] - 1.);
2758 if (std::strncmp(type__, "REGULR", 6))
2760 F77_FUNC(dger, DGER)
2761 (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[z_offset], ldz);
2770 /* Selected single precision arpack routines */
2773 static void F77_FUNC(sstqrb, SSTQRB)(int* n, float* d__, float* e, float* z__, float* work, int* info)
2782 int i__, j, k, l, m;
2784 int l1, ii, mm, lm1, mm1, nm1;
2785 float rt1, rt2, eps;
2788 int lend, jtot, lendm1, lendp1, iscale;
2790 int lendsv, nmaxit, icompz;
2791 float ssfmax, ssfmin, safmin, minval, safmax, anorm;
2814 eps = GMX_FLOAT_EPS;
2818 minval = GMX_FLOAT_MIN;
2819 safmin = minval / GMX_FLOAT_EPS;
2820 safmax = 1. / safmin;
2821 ssfmax = std::sqrt(safmax) / 3.;
2822 ssfmin = std::sqrt(safmin) / eps2;
2825 for (j = 1; j <= i__1; ++j)
2849 for (m = l1; m <= i__1; ++m)
2851 tst = std::abs(e[m]);
2856 if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m + 1])) * eps)
2876 i__1 = lend - l + 1;
2877 anorm = F77_FUNC(slanst, SLANST)("i", &i__1, &d__[l], &e[l]);
2886 i__1 = lend - l + 1;
2887 F77_FUNC(slascl, SLASCL)
2888 ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, info);
2890 F77_FUNC(slascl, SLASCL)("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, info);
2892 else if (anorm < ssfmin)
2895 i__1 = lend - l + 1;
2896 F77_FUNC(slascl, SLASCL)
2897 ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, info);
2899 F77_FUNC(slascl, SLASCL)("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, info);
2902 if (std::abs(d__[lend]) < std::abs(d__[l]))
2916 for (m = l; m <= i__1; ++m)
2918 d__2 = std::abs(e[m]);
2920 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
2944 F77_FUNC(slaev2, SLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
2946 work[*n - 1 + l] = s;
2949 z__[l + 1] = c__ * tst - s * z__[l];
2950 z__[l] = s * tst + c__ * z__[l];
2954 F77_FUNC(slae2, SLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
2973 g = (d__[l + 1] - p) / (e[l] * 2.);
2974 r__ = F77_FUNC(slapy2, SLAPY2)(&g, &c_b31);
2975 g = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__));
2983 for (i__ = mm1; i__ >= i__1; --i__)
2987 F77_FUNC(slartg, SLARTG)(&g, &f, &c__, &s, &r__);
2992 g = d__[i__ + 1] - p;
2993 r__ = (d__[i__] - g) * s + c__ * 2. * b;
2995 d__[i__ + 1] = g + p;
3001 work[*n - 1 + i__] = -s;
3009 F77_FUNC(slasr, SLASR)
3010 ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &z__[l], &c__1);
3035 for (m = l; m >= i__1; --m)
3037 d__2 = std::abs(e[m - 1]);
3039 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m - 1]) + safmin)
3063 F77_FUNC(slaev2, SLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s);
3066 z__[l] = c__ * tst - s * z__[l - 1];
3067 z__[l - 1] = s * tst + c__ * z__[l - 1];
3071 F77_FUNC(slae2, SLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
3091 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
3092 r__ = F77_FUNC(slapy2, SLAPY2)(&g, &c_b31);
3093 g = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__));
3101 for (i__ = m; i__ <= i__1; ++i__)
3105 F77_FUNC(slartg, SLARTG)(&g, &f, &c__, &s, &r__);
3111 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
3119 work[*n - 1 + i__] = s;
3127 F77_FUNC(slasr, SLASR)
3128 ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &z__[m], &c__1);
3149 i__1 = lendsv - lsv + 1;
3150 F77_FUNC(slascl, SLASCL)
3151 ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], n, info);
3152 i__1 = lendsv - lsv;
3153 F77_FUNC(slascl, SLASCL)
3154 ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n, info);
3156 else if (iscale == 2)
3158 i__1 = lendsv - lsv + 1;
3159 F77_FUNC(slascl, SLASCL)
3160 ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], n, info);
3161 i__1 = lendsv - lsv;
3162 F77_FUNC(slascl, SLASCL)
3163 ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n, info);
3171 for (i__ = 1; i__ <= i__1; ++i__)
3184 F77_FUNC(slasrt, SLASRT)("i", n, &d__[1], info);
3190 for (ii = 2; ii <= i__1; ++ii)
3196 for (j = ii; j <= i__2; ++j)
3220 static void F77_FUNC(sgetv0, SGETV0)(int* ido,
3222 int gmx_unused* itry,
3239 int v_dim1, v_offset, i__1;
3247 v_offset = 1 + v_dim1;
3263 F77_FUNC(slarnv, SLARNV)(&idist, &iwork[1], n, &resid[1]);
3270 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
3289 F77_FUNC(scopy, SCOPY)(n, &workd[*n + 1], &c__1, &resid[1], &c__1);
3295 else if (*bmat == 'I')
3297 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
3306 workd[*n * 3 + 4] = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
3307 workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
3309 else if (*bmat == 'I')
3311 workd[*n * 3 + 4] = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
3313 *rnorm = workd[*n * 3 + 4];
3323 F77_FUNC(sgemv, SGEMV)
3324 ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24, &workd[*n + 1], &c__1);
3326 F77_FUNC(sgemv, SGEMV)
3327 ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &c_b22, &resid[1], &c__1);
3331 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3337 else if (*bmat == 'I')
3339 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
3346 *rnorm = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
3347 *rnorm = std::sqrt(std::abs(*rnorm));
3349 else if (*bmat == 'I')
3351 *rnorm = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
3354 if (*rnorm > workd[*n * 3 + 4] * .717F)
3363 workd[*n * 3 + 4] = *rnorm;
3370 for (jj = 1; jj <= i__1; ++jj)
3387 static void F77_FUNC(ssapps, SSAPPS)(int* n,
3404 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2, i__3, i__4;
3407 float r__, s, a1, a2, a3, a4;
3418 v_offset = 1 + v_dim1;
3421 h_offset = 1 + h_dim1;
3424 q_offset = 1 + q_dim1;
3427 epsmch = GMX_FLOAT_EPS;
3431 kplusp = *kev + *np;
3433 F77_FUNC(slaset, SLASET)("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3441 for (jj = 1; jj <= i__1; ++jj)
3449 for (i__ = istart; i__ <= i__2; ++i__)
3451 big = std::abs(h__[i__ + (h_dim1 * 2)]) + std::abs(h__[i__ + 1 + (h_dim1 * 2)]);
3452 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3454 h__[i__ + 1 + h_dim1] = 0.;
3465 f = h__[istart + (h_dim1 << 1)] - shift[jj];
3466 g = h__[istart + 1 + h_dim1];
3467 F77_FUNC(slartg, SLARTG)(&f, &g, &c__, &s, &r__);
3469 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 + h_dim1];
3470 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (h_dim1 << 1)];
3471 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 + h_dim1];
3472 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 << 1)];
3473 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
3474 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3475 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
3478 i__2 = (i__3 < kplusp) ? i__3 : kplusp;
3479 for (j = 1; j <= i__2; ++j)
3481 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) * q_dim1];
3482 q[j + (istart + 1) * q_dim1] =
3483 -s * q[j + istart * q_dim1] + c__ * q[j + (istart + 1) * q_dim1];
3484 q[j + istart * q_dim1] = a1;
3488 for (i__ = istart + 1; i__ <= i__2; ++i__)
3491 f = h__[i__ + h_dim1];
3492 g = s * h__[i__ + 1 + h_dim1];
3494 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3495 F77_FUNC(slartg, SLARTG)(&f, &g, &c__, &s, &r__);
3504 h__[i__ + h_dim1] = r__;
3506 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 + h_dim1];
3507 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1 << 1)];
3508 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)];
3509 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 + h_dim1];
3511 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
3512 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3513 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
3516 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
3517 for (j = 1; j <= i__3; ++j)
3519 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) * q_dim1];
3520 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] + c__ * q[j + (i__ + 1) * q_dim1];
3521 q[j + i__ * q_dim1] = a1;
3528 if (h__[iend + h_dim1] < 0.)
3530 h__[iend + h_dim1] = -h__[iend + h_dim1];
3531 F77_FUNC(sscal, SSCAL)(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3540 for (i__ = itop; i__ <= i__2; ++i__)
3542 if (h__[i__ + 1 + h_dim1] > 0.)
3553 for (i__ = itop; i__ <= i__1; ++i__)
3555 big = std::abs(h__[i__ + (h_dim1 * 2)]) + std::abs(h__[i__ + 1 + (h_dim1 * 2)]);
3556 if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3558 h__[i__ + 1 + h_dim1] = 0.;
3562 if (h__[*kev + 1 + h_dim1] > 0.)
3564 F77_FUNC(sgemv, SGEMV)
3565 ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) * q_dim1 + 1], &c__1, &c_b4,
3566 &workd[*n + 1], &c__1);
3570 for (i__ = 1; i__ <= i__1; ++i__)
3572 i__2 = kplusp - i__ + 1;
3573 F77_FUNC(sgemv, SGEMV)
3574 ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) * q_dim1 + 1], &c__1, &c_b4,
3576 F77_FUNC(scopy, SCOPY)(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &c__1);
3579 F77_FUNC(slacpy, SLACPY)("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3581 if (h__[*kev + 1 + h_dim1] > 0.)
3583 F77_FUNC(scopy, SCOPY)(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3586 F77_FUNC(sscal, SSCAL)(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3587 if (h__[*kev + 1 + h_dim1] > 0.)
3589 F77_FUNC(saxpy, SAXPY)
3590 (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1, &resid[1], &c__1);
3599 static void F77_FUNC(ssortr, SSORTR)(const char* which, int* apply, int* n, float* x1, float* x2)
3609 if (!std::strncmp(which, "SA", 2))
3618 for (i__ = igap; i__ <= i__1; ++i__)
3628 if (x1[j] < x1[j + igap])
3631 x1[j] = x1[j + igap];
3632 x1[j + igap] = temp;
3636 x2[j] = x2[j + igap];
3637 x2[j + igap] = temp;
3651 else if (!std::strncmp(which, "SM", 2))
3660 for (i__ = igap; i__ <= i__1; ++i__)
3670 if (std::abs(x1[j]) < std::abs(x1[j + igap]))
3673 x1[j] = x1[j + igap];
3674 x1[j + igap] = temp;
3678 x2[j] = x2[j + igap];
3679 x2[j + igap] = temp;
3693 else if (!std::strncmp(which, "LA", 2))
3702 for (i__ = igap; i__ <= i__1; ++i__)
3712 if (x1[j] > x1[j + igap])
3715 x1[j] = x1[j + igap];
3716 x1[j + igap] = temp;
3720 x2[j] = x2[j + igap];
3721 x2[j + igap] = temp;
3735 else if (!std::strncmp(which, "LM", 2))
3745 for (i__ = igap; i__ <= i__1; ++i__)
3755 if (std::abs(x1[j]) > std::abs(x1[j + igap]))
3758 x1[j] = x1[j + igap];
3759 x1[j + igap] = temp;
3763 x2[j] = x2[j + igap];
3764 x2[j + igap] = temp;
3784 static void F77_FUNC(ssesrt,
3785 SSESRT)(const char* which, int* apply, int* n, float* x, int* na, float* a, int* lda)
3787 int a_dim1, a_offset, i__1;
3794 a_offset = 1 + a_dim1 * 0;
3799 if (!std::strncmp(which, "SA", 2))
3808 for (i__ = igap; i__ <= i__1; ++i__)
3818 if (x[j] < x[j + igap])
3825 F77_FUNC(sswap, SSWAP)
3826 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
3840 else if (!std::strncmp(which, "SM", 2))
3849 for (i__ = igap; i__ <= i__1; ++i__)
3859 if (std::abs(x[j]) < std::abs(x[j + igap]))
3866 F77_FUNC(sswap, SSWAP)
3867 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
3881 else if (!std::strncmp(which, "LA", 2))
3890 for (i__ = igap; i__ <= i__1; ++i__)
3900 if (x[j] > x[j + igap])
3907 F77_FUNC(sswap, SSWAP)
3908 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
3922 else if (!std::strncmp(which, "LM", 2))
3931 for (i__ = igap; i__ <= i__1; ++i__)
3941 if (std::abs(x[j]) > std::abs(x[j + igap]))
3948 F77_FUNC(sswap, SSWAP)
3949 (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) * a_dim1 + 1], &c__1);
3969 static void F77_FUNC(ssgets,
3970 SSGETS)(int* ishift, const char* which, int* kev, int* np, float* ritz, float* bounds, float* shifts)
3980 if (!std::strncmp(which, "BE", 2))
3983 F77_FUNC(ssortr, SSORTR)("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
3987 i__1 = (kevd2 < *np) ? kevd2 : *np;
3988 i__2 = (kevd2 > *np) ? kevd2 : *np;
3989 F77_FUNC(sswap, SSWAP)(&i__1, &ritz[1], &c__1, &ritz[i__2 + 1], &c__1);
3990 i__1 = (kevd2 < *np) ? kevd2 : *np;
3991 i__2 = (kevd2 > *np) ? kevd2 : *np;
3992 F77_FUNC(sswap, SSWAP)(&i__1, &bounds[1], &c__1, &bounds[i__2 + 1], &c__1);
3998 F77_FUNC(ssortr, SSORTR)(which, &c__1, &i__1, &ritz[1], &bounds[1]);
4001 if (*ishift == 1 && *np > 0)
4004 F77_FUNC(ssortr, SSORTR)("SM", &c__1, np, &bounds[1], &ritz[1]);
4005 F77_FUNC(scopy, SCOPY)(np, &ritz[1], &c__1, &shifts[1], &c__1);
4013 static void F77_FUNC(ssconv, SSCONV)(int* n, float* ritz, float* bounds, float* tol, int* nconv)
4015 float c_b3 = 2 / 3.;
4025 eps23 = GMX_FLOAT_EPS;
4026 eps23 = std::pow(eps23, c_b3);
4030 for (i__ = 1; i__ <= i__1; ++i__)
4034 d__3 = std::abs(ritz[i__]);
4035 temp = (d__2 > d__3) ? d__2 : d__3;
4036 if (bounds[i__] <= *tol * temp)
4046 static void F77_FUNC(
4048 SSEIGT)(float* rnorm, int* n, float* h__, int* ldh, float* eig, float* bounds, float* workl, int* ierr)
4051 int h_dim1, h_offset, i__1;
4060 h_offset = 1 + h_dim1;
4063 F77_FUNC(scopy, SCOPY)(n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
4065 F77_FUNC(scopy, SCOPY)(&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
4066 F77_FUNC(sstqrb, SSTQRB)(n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
4073 for (k = 1; k <= i__1; ++k)
4075 bounds[k] = *rnorm * std::abs(bounds[k]);
4084 static void F77_FUNC(ssaitr, SSAITR)(int* ido,
4108 int h_dim1, h_offset, v_dim1, v_offset, i__1;
4112 float safmin, minval;
4118 v_offset = 1 + v_dim1;
4121 h_offset = 1 + h_dim1;
4125 minval = GMX_FLOAT_MIN;
4126 safmin = minval / GMX_FLOAT_EPS;
4140 iwork[9] = iwork[8] + *n;
4141 iwork[10] = iwork[9] + *n;
4179 F77_FUNC(sgetv0, sgetv0)
4180 (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv, &resid[1], rnorm, &ipntr[1],
4181 &workd[1], &iwork[21], &iwork[7]);
4194 *info = iwork[12] - 1;
4201 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
4202 if (*rnorm >= safmin)
4204 temp1 = 1. / *rnorm;
4205 F77_FUNC(sscal, SSCAL)(n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
4206 F77_FUNC(sscal, SSCAL)(n, &temp1, &workd[iwork[8]], &c__1);
4211 F77_FUNC(slascl, SLASCL)
4212 ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] * v_dim1 + 1], n, &infol);
4213 F77_FUNC(slascl, SLASCL)
4214 ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[8]], n, &infol);
4218 F77_FUNC(scopy, SCOPY)(n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
4219 ipntr[1] = iwork[10];
4220 ipntr[2] = iwork[9];
4221 ipntr[3] = iwork[8];
4230 F77_FUNC(scopy, SCOPY)(n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
4239 ipntr[1] = iwork[9];
4240 ipntr[2] = iwork[8];
4245 else if (*bmat == 'I')
4247 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4257 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[iwork[10]], &c__1);
4258 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4260 else if (*bmat == 'G')
4262 workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4263 workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4265 else if (*bmat == 'I')
4267 workd[*n * 3 + 3] = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
4272 F77_FUNC(sgemv, SGEMV)
4273 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &c__1, &c_b42,
4274 &workd[iwork[9]], &c__1);
4278 F77_FUNC(sgemv, SGEMV)
4279 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]], &c__1, &c_b42,
4280 &workd[iwork[9]], &c__1);
4283 F77_FUNC(sgemv, SGEMV)
4284 ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &c__1, &c_b18, &resid[1], &c__1);
4286 h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
4287 if (iwork[12] == 1 || iwork[4] == 1)
4289 h__[iwork[12] + h_dim1] = 0.;
4293 h__[iwork[12] + h_dim1] = *rnorm;
4301 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4302 ipntr[1] = iwork[9];
4303 ipntr[2] = iwork[8];
4308 else if (*bmat == 'I')
4310 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4318 *rnorm = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4319 *rnorm = std::sqrt(std::abs(*rnorm));
4321 else if (*bmat == 'I')
4323 *rnorm = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
4326 if (*rnorm > workd[*n * 3 + 3] * .717F)
4333 F77_FUNC(sgemv, SGEMV)
4334 ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &c__1, &c_b42,
4335 &workd[iwork[9]], &c__1);
4337 F77_FUNC(sgemv, SGEMV)
4338 ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &c__1, &c_b18, &resid[1], &c__1);
4340 if (iwork[12] == 1 || iwork[4] == 1)
4342 h__[iwork[12] + h_dim1] = 0.;
4344 h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
4349 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4350 ipntr[1] = iwork[9];
4351 ipntr[2] = iwork[8];
4356 else if (*bmat == 'I')
4358 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4365 workd[*n * 3 + 2] = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4366 workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
4368 else if (*bmat == 'I')
4370 workd[*n * 3 + 2] = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
4374 if (workd[*n * 3 + 2] > *rnorm * .717F)
4377 *rnorm = workd[*n * 3 + 2];
4382 *rnorm = workd[*n * 3 + 2];
4390 for (jj = 1; jj <= i__1; ++jj)
4402 if (h__[iwork[12] + h_dim1] < 0.)
4404 h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4405 if (iwork[12] < *k + *np)
4407 F77_FUNC(sscal, SSCAL)(n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4411 F77_FUNC(sscal, SSCAL)(n, &c_b50, &resid[1], &c__1);
4416 if (iwork[12] > *k + *np)
4431 static void F77_FUNC(ssaup2, SSAUP2)(int* ido,
4440 int gmx_unused* iupd,
4457 float c_b3 = 2 / 3.;
4461 int h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2, i__3;
4479 v_offset = 1 + v_dim1;
4482 h_offset = 1 + h_dim1;
4485 q_offset = 1 + q_dim1;
4489 eps23 = GMX_FLOAT_EPS;
4490 eps23 = std::pow(eps23, c_b3);
4503 iwork[7] = iwork[9] + iwork[10];
4526 F77_FUNC(sgetv0, SGETV0)
4527 (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &resid[1], &workd[*n * 3 + 1],
4528 &ipntr[1], &workd[1], &iwork[41], info);
4535 if (workd[*n * 3 + 1] == 0.)
4560 F77_FUNC(ssaitr, SSAITR)
4561 (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 + 1], &v[v_offset], ldv,
4562 &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[21], info);
4587 F77_FUNC(ssaitr, SSAITR)
4588 (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[v_offset], ldv, &h__[h_offset],
4589 ldh, &ipntr[1], &workd[1], &iwork[21], info);
4606 F77_FUNC(sseigt, SSEIGT)
4607 (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &bounds[1], &workl[1], &ierr);
4615 F77_FUNC(scopy, SCOPY)(&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4616 F77_FUNC(scopy, SCOPY)(&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4620 F77_FUNC(ssgets, SSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4622 F77_FUNC(scopy, SCOPY)(nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4623 F77_FUNC(ssconv, SSCONV)(nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4628 for (j = 1; j <= i__1; ++j)
4630 if (bounds[j] == 0.)
4637 if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
4640 if (!std::strncmp(which, "BE", 2))
4643 std::strncpy(wprime, "SA", 2);
4644 F77_FUNC(ssortr, SSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4646 nevm2 = *nev - nevd2;
4649 i__1 = (nevd2 < *np) ? nevd2 : *np;
4650 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4651 F77_FUNC(sswap, SSWAP)
4652 (&i__1, &ritz[nevm2 + 1], &c__1, &ritz[((i__2 > i__3) ? i__2 : i__3)], &c__1);
4653 i__1 = (nevd2 < *np) ? nevd2 : *np;
4654 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4655 F77_FUNC(sswap, SSWAP)
4656 (&i__1, &bounds[nevm2 + 1], &c__1, &bounds[((i__2 > i__3) ? i__2 : i__3) + 1], &c__1);
4662 if (!std::strncmp(which, "LM", 2))
4664 std::strncpy(wprime, "SM", 2);
4666 if (!std::strncmp(which, "SM", 2))
4668 std::strncpy(wprime, "LM", 2);
4670 if (!std::strncmp(which, "LA", 2))
4672 std::strncpy(wprime, "SA", 2);
4674 if (!std::strncmp(which, "SA", 2))
4676 std::strncpy(wprime, "LA", 2);
4679 F77_FUNC(ssortr, SSORTR)(wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4683 for (j = 1; j <= i__1; ++j)
4686 d__3 = std::abs(ritz[j]);
4687 temp = (d__2 > d__3) ? d__2 : d__3;
4691 std::strncpy(wprime, "LA", 2);
4692 F77_FUNC(ssortr, SSORTR)(wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4695 for (j = 1; j <= i__1; ++j)
4698 d__3 = std::abs(ritz[j]);
4699 temp = (d__2 > d__3) ? d__2 : d__3;
4703 if (!std::strncmp(which, "BE", 2))
4706 std::strncpy(wprime, "LA", 2);
4707 F77_FUNC(ssortr, SSORTR)(wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4711 F77_FUNC(ssortr, SSORTR)(which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4714 h__[h_dim1 + 1] = workd[*n * 3 + 1];
4717 if (iwork[6] > *mxiter && iwork[8] < *nev)
4721 if (*np == 0 && iwork[8] < iwork[9])
4729 else if (iwork[8] < *nev && *ishift == 1)
4732 i__1 = iwork[8], i__2 = *np / 2;
4733 *nev += (i__1 < i__2) ? i__1 : i__2;
4734 if (*nev == 1 && iwork[7] >= 6)
4736 *nev = iwork[7] / 2;
4738 else if (*nev == 1 && iwork[7] > 2)
4742 *np = iwork[7] - *nev;
4747 F77_FUNC(ssgets, SSGETS)(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4766 F77_FUNC(scopy, SCOPY)(np, &workl[1], &c__1, &ritz[1], &c__1);
4769 F77_FUNC(ssapps, SSAPPS)
4770 (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &resid[1], &q[q_offset], ldq,
4776 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
4783 else if (*bmat == 'I')
4785 F77_FUNC(scopy, SCOPY)(n, &resid[1], &c__1, &workd[1], &c__1);
4792 workd[*n * 3 + 1] = F77_FUNC(sdot, SDOT)(n, &resid[1], &c__1, &workd[1], &c__1);
4793 workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
4795 else if (*bmat == 'I')
4797 workd[*n * 3 + 1] = F77_FUNC(snrm2, SNRM2)(n, &resid[1], &c__1);
4816 void F77_FUNC(ssaupd, SSAUPD)(int* ido,
4834 int v_dim1, v_offset, i__1, i__2;
4840 v_offset = 1 + v_dim1;
4852 iwork[5] = iparam[1];
4853 iwork[10] = iparam[3];
4854 iwork[12] = iparam[4];
4857 iwork[11] = iparam[7];
4868 else if (*ncv <= *nev || *ncv > *n)
4874 iwork[15] = *ncv - *nev;
4880 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) && std::strncmp(which, "LA", 2)
4881 && std::strncmp(which, "SA", 2) && std::strncmp(which, "BE", 2))
4885 if (*bmat != 'I' && *bmat != 'G')
4891 if (*lworkl < i__1 * i__1 + (*ncv << 3))
4895 if (iwork[11] < 1 || iwork[11] > 5)
4899 else if (iwork[11] == 1 && *bmat == 'G')
4903 else if (iwork[5] < 0 || iwork[5] > 1)
4907 else if (*nev == 1 && !std::strncmp(which, "BE", 2))
4925 *tol = GMX_FLOAT_EPS;
4928 iwork[15] = *ncv - *nev;
4931 i__1 = i__2 * i__2 + (*ncv << 3);
4932 for (j = 1; j <= i__1; ++j)
4940 iwork[16] = iwork[3] + (iwork[8] << 1);
4941 iwork[1] = iwork[16] + *ncv;
4942 iwork[4] = iwork[1] + *ncv;
4944 iwork[7] = iwork[4] + i__1 * i__1;
4945 iwork[14] = iwork[7] + *ncv * 3;
4947 ipntr[4] = iwork[14];
4948 ipntr[5] = iwork[3];
4949 ipntr[6] = iwork[16];
4950 ipntr[7] = iwork[1];
4951 ipntr[11] = iwork[7];
4954 F77_FUNC(ssaup2, SSAUP2)
4955 (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &iwork[11], &iwork[6], &iwork[5],
4956 &iwork[10], &v[v_offset], ldv, &workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]],
4957 &workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1], &iwork[21], info);
4961 iparam[8] = iwork[15];
4968 iparam[3] = iwork[10];
4969 iparam[5] = iwork[15];
4986 void F77_FUNC(sseupd, SSEUPD)(int* rvec,
5009 float c_b21 = 2 / 3.;
5012 int v_dim1, v_offset, z_dim1, z_offset, i__1;
5013 float d__1, d__2, d__3;
5015 int j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
5027 float thres1 = 0, thres2 = 0;
5031 int leftptr, rghtptr;
5037 z_offset = 1 + z_dim1;
5042 v_offset = 1 + v_dim1;
5070 if (*ncv <= *nev || *ncv > *n)
5074 if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) && std::strncmp(which, "LA", 2)
5075 && std::strncmp(which, "SA", 2) && std::strncmp(which, "BE", 2))
5079 if (*bmat != 'I' && *bmat != 'G')
5083 if (*howmny != 'A' && *howmny != 'P' && *howmny != 'S' && *rvec)
5087 if (*rvec && *howmny == 'S')
5092 if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
5097 if (mode == 1 || mode == 2)
5099 std::strncpy(type__, "REGULR", 6);
5103 std::strncpy(type__, "SHIFTI", 6);
5107 std::strncpy(type__, "BUCKLE", 6);
5111 std::strncpy(type__, "CAYLEY", 6);
5117 if (mode == 1 && *bmat == 'G')
5121 if (*nev == 1 && !std::strncmp(which, "BE", 2))
5140 iw = iq + ldh * *ncv;
5141 next = iw + (*ncv << 1);
5147 irz = ipntr[11] + *ncv;
5151 eps23 = GMX_FLOAT_EPS;
5152 eps23 = std::pow(eps23, c_b21);
5159 else if (*bmat == 'G')
5161 bnorm2 = F77_FUNC(snrm2, SNRM2)(n, &workd[1], &c__1);
5167 if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2)
5168 || !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
5170 else if (!std::strncmp(which, "BE", 2))
5174 ism = (*nev > nconv) ? *nev : nconv;
5177 thres1 = workl[ism];
5178 thres2 = workl[ilg];
5184 for (j = 0; j <= i__1; ++j)
5187 if (!std::strncmp(which, "LM", 2))
5189 if (std::abs(workl[irz + j]) >= std::abs(thres1))
5192 d__3 = std::abs(workl[irz + j]);
5193 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5194 if (workl[ibd + j] <= *tol * tempbnd)
5200 else if (!std::strncmp(which, "SM", 2))
5202 if (std::abs(workl[irz + j]) <= std::abs(thres1))
5205 d__3 = std::abs(workl[irz + j]);
5206 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5207 if (workl[ibd + j] <= *tol * tempbnd)
5213 else if (!std::strncmp(which, "LA", 2))
5215 if (workl[irz + j] >= thres1)
5218 d__3 = std::abs(workl[irz + j]);
5219 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5220 if (workl[ibd + j] <= *tol * tempbnd)
5226 else if (!std::strncmp(which, "SA", 2))
5228 if (workl[irz + j] <= thres1)
5231 d__3 = std::abs(workl[irz + j]);
5232 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5233 if (workl[ibd + j] <= *tol * tempbnd)
5239 else if (!std::strncmp(which, "BE", 2))
5241 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
5244 d__3 = std::abs(workl[irz + j]);
5245 tempbnd = (d__2 > d__3) ? d__2 : d__3;
5246 if (workl[ibd + j] <= *tol * tempbnd)
5254 reord = select[j + 1] || reord;
5263 F77_FUNC(scopy, SCOPY)(&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
5264 F77_FUNC(scopy, SCOPY)(ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
5266 F77_FUNC(ssteqr, SSTEQR)
5267 ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &workl[iw], &ierr);
5288 if (select[leftptr])
5293 else if (!select[rghtptr])
5301 temp = workl[ihd + leftptr - 1];
5302 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
5303 workl[ihd + rghtptr - 1] = temp;
5304 F77_FUNC(scopy, SCOPY)
5305 (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[iw], &c__1);
5306 F77_FUNC(scopy, SCOPY)
5307 (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[iq + *ncv * (leftptr - 1)], &c__1);
5308 F77_FUNC(scopy, SCOPY)
5309 (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr - 1)], &c__1);
5314 if (leftptr < rghtptr)
5322 F77_FUNC(scopy, SCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5327 F77_FUNC(scopy, SCOPY)(&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
5328 F77_FUNC(scopy, SCOPY)(ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
5330 if (!std::strncmp(type__, "REGULR", 6))
5335 F77_FUNC(ssesrt, SSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5339 F77_FUNC(scopy, SCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5345 F77_FUNC(scopy, SCOPY)(ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
5346 if (!std::strncmp(type__, "SHIFTI", 6))
5349 for (k = 1; k <= i__1; ++k)
5351 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
5354 else if (!std::strncmp(type__, "BUCKLE", 6))
5357 for (k = 1; k <= i__1; ++k)
5359 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd + k - 1] - 1.);
5362 else if (!std::strncmp(type__, "CAYLEY", 6))
5365 for (k = 1; k <= i__1; ++k)
5367 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (workl[ihd + k - 1] - 1.);
5371 F77_FUNC(scopy, SCOPY)(&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5372 F77_FUNC(ssortr, SSORTR)("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
5375 F77_FUNC(ssesrt, SSESRT)("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5379 F77_FUNC(scopy, SCOPY)(ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5380 d__1 = bnorm2 / rnorm;
5381 F77_FUNC(sscal, SSCAL)(ncv, &d__1, &workl[ihb], &c__1);
5382 F77_FUNC(ssortr, SSORTR)("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
5386 if (*rvec && *howmny == 'A')
5389 F77_FUNC(sgeqr2, SGEQR2)
5390 (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb], &ierr);
5392 F77_FUNC(sorm2r, SORM2R)
5393 ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &v[v_offset],
5394 ldv, &workd[*n + 1], &ierr);
5395 F77_FUNC(slacpy, SLACPY)("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
5398 for (j = 1; j <= i__1; ++j)
5400 workl[ihb + j - 1] = 0.;
5402 workl[ihb + *ncv - 1] = 1.;
5403 F77_FUNC(sorm2r, SORM2R)
5404 ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
5407 else if (*rvec && *howmny == 'S')
5411 if (!std::strncmp(type__, "REGULR", 6) && *rvec)
5415 for (j = 1; j <= i__1; ++j)
5417 workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
5420 else if (std::strncmp(type__, "REGULR", 6) && *rvec)
5423 F77_FUNC(sscal, SSCAL)(ncv, &bnorm2, &workl[ihb], &c__1);
5424 if (!std::strncmp(type__, "SHIFTI", 6))
5428 for (k = 1; k <= i__1; ++k)
5430 d__2 = workl[iw + k - 1];
5431 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1]) / (d__2 * d__2);
5434 else if (!std::strncmp(type__, "BUCKLE", 6))
5438 for (k = 1; k <= i__1; ++k)
5440 d__2 = workl[iw + k - 1] - 1.;
5441 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1]) / (d__2 * d__2);
5444 else if (!std::strncmp(type__, "CAYLEY", 6))
5448 for (k = 1; k <= i__1; ++k)
5450 workl[ihb + k - 1] =
5451 std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
5456 if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
5460 for (k = 0; k <= i__1; ++k)
5462 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
5465 else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
5469 for (k = 0; k <= i__1; ++k)
5471 workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] - 1.);
5475 if (std::strncmp(type__, "REGULR", 6))
5477 F77_FUNC(sger, SGER)
5478 (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[z_offset], ldz);