2 #include "gromacs/utility/real.h"
4 #include "../gmx_blas.h"
5 #include "../gmx_lapack.h"
6 #include "lapack_limits.h"
9 F77_FUNC(dsteqr,DSTEQR)(const char * compz,
23 int z_dim1, z_offset, i__1, i__2;
29 int l1, ii, mm, lm1, mm1, nm1;
48 z_offset = 1 + z_dim1;
54 if (*compz=='N' || *compz=='n') {
56 } else if (*compz=='V' || *compz=='v') {
58 } else if (*compz=='I' || *compz=='i') {
67 } else if (*ldz < 1 || (icompz > 0 && *ldz < ((*n>1) ? *n : 1))) {
89 minval = GMX_DOUBLE_MIN;
90 safmin = minval*(1.0+GMX_DOUBLE_EPS);
93 ssfmax = sqrt(safmax) / 3.;
94 ssfmin = sqrt(safmin) / eps2;
97 F77_FUNC(dlaset,DLASET)("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
115 for (m = l1; m <= i__1; ++m) {
117 if (fabs(tst)<GMX_DOUBLE_MIN) {
120 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m + 1])) * eps) {
139 anorm = F77_FUNC(dlanst,DLANST)("I", &i__1, &d__[l], &e[l]);
141 if (fabs(anorm)<GMX_DOUBLE_MIN) {
144 if (anorm > ssfmax) {
147 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
150 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
152 } else if (anorm < ssfmin) {
155 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
158 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
162 if (fabs(d__[lend]) < fabs(d__[l])) {
173 for (m = l; m <= i__1; ++m) {
176 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m+ 1]) + safmin) {
195 F77_FUNC(dlaev2,DLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
197 work[*n - 1 + l] = s;
198 F77_FUNC(dlasr,DLASR)("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
199 z__[l * z_dim1 + 1], ldz);
201 F77_FUNC(dlae2,DLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
213 if (jtot == nmaxit) {
218 g = (d__[l + 1] - p) / (e[l] * 2.);
219 r__ = F77_FUNC(dlapy2,DLAPY2)(&g, &c_b10);
220 g = d__[m] - p + e[l] / (g + ( (g>0) ? r__ : -r__ ) );
228 for (i__ = mm1; i__ >= i__1; --i__) {
231 F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
235 g = d__[i__ + 1] - p;
236 r__ = (d__[i__] - g) * s + c__ * 2. * b;
238 d__[i__ + 1] = g + p;
243 work[*n - 1 + i__] = -s;
249 F77_FUNC(dlasr,DLASR)("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
272 for (m = l; m >= i__1; --m) {
273 d__2 = fabs(e[m - 1]);
275 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin) {
293 F77_FUNC(dlaev2,DLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
296 work[*n - 1 + m] = s;
297 F77_FUNC(dlasr,DLASR)("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
298 z__[(l - 1) * z_dim1 + 1], ldz);
300 F77_FUNC(dlae2,DLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
312 if (jtot == nmaxit) {
317 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
318 r__ = F77_FUNC(dlapy2,DLAPY2)(&g, &c_b10);
319 g = d__[m] - p + e[l - 1] / (g + ( (g>0) ? r__ : -r__ ));
327 for (i__ = m; i__ <= i__1; ++i__) {
330 F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
335 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
342 work[*n - 1 + i__] = s;
348 F77_FUNC(dlasr,DLASR)("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
369 i__1 = lendsv - lsv + 1;
370 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
373 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
375 } else if (iscale == 2) {
376 i__1 = lendsv - lsv + 1;
377 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
380 F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
388 for (i__ = 1; i__ <= i__1; ++i__) {
389 if (fabs(e[i__])>GMX_DOUBLE_MIN) {
398 F77_FUNC(dlasrt,DLASRT)("I", n, &d__[1], info);
403 for (ii = 2; ii <= i__1; ++ii) {
408 for (j = ii; j <= i__2; ++j) {
417 F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],