2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
36 #include <types/simple.h>
39 #include "gmx_lapack.h"
40 #include "lapack_limits.h"
43 F77_FUNC(ssteqr,SSTEQR)(const char * compz,
57 int z_dim1, z_offset, i__1, i__2;
63 int l1, ii, mm, lm1, mm1, nm1;
82 z_offset = 1 + z_dim1;
88 if (*compz=='N' || *compz=='n') {
90 } else if (*compz=='V' || *compz=='v') {
92 } else if (*compz=='I' || *compz=='i') {
101 } else if (*ldz < 1 || (icompz > 0 && *ldz < ((*n>1) ? *n : 1))) {
115 z__[z_dim1 + 1] = 1.;
123 minval = GMX_FLOAT_MIN;
124 safmin = minval*(1.0+GMX_FLOAT_EPS);
126 safmax = 1. / safmin;
127 ssfmax = sqrt(safmax) / 3.;
128 ssfmin = sqrt(safmin) / eps2;
131 F77_FUNC(slaset,SLASET)("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
149 for (m = l1; m <= i__1; ++m) {
151 if (fabs(tst)<GMX_FLOAT_MIN) {
154 if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m + 1])) * eps) {
173 anorm = F77_FUNC(slanst,SLANST)("I", &i__1, &d__[l], &e[l]);
175 if (fabs(anorm)<GMX_FLOAT_MIN) {
178 if (anorm > ssfmax) {
181 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
184 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
186 } else if (anorm < ssfmin) {
189 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
192 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
196 if (fabs(d__[lend]) < fabs(d__[l])) {
207 for (m = l; m <= i__1; ++m) {
210 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m+ 1]) + safmin) {
229 F77_FUNC(slaev2,SLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
231 work[*n - 1 + l] = s;
232 F77_FUNC(slasr,SLASR)("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
233 z__[l * z_dim1 + 1], ldz);
235 F77_FUNC(slae2,SLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
247 if (jtot == nmaxit) {
252 g = (d__[l + 1] - p) / (e[l] * 2.);
253 r__ = F77_FUNC(slapy2,SLAPY2)(&g, &c_b10);
254 g = d__[m] - p + e[l] / (g + ( (g>0) ? r__ : -r__ ) );
262 for (i__ = mm1; i__ >= i__1; --i__) {
265 F77_FUNC(slartg,SLARTG)(&g, &f, &c__, &s, &r__);
269 g = d__[i__ + 1] - p;
270 r__ = (d__[i__] - g) * s + c__ * 2. * b;
272 d__[i__ + 1] = g + p;
277 work[*n - 1 + i__] = -s;
283 F77_FUNC(slasr,SLASR)("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
306 for (m = l; m >= i__1; --m) {
307 d__2 = fabs(e[m - 1]);
309 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin) {
327 F77_FUNC(slaev2,SLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
330 work[*n - 1 + m] = s;
331 F77_FUNC(slasr,SLASR)("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
332 z__[(l - 1) * z_dim1 + 1], ldz);
334 F77_FUNC(slae2,SLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
346 if (jtot == nmaxit) {
351 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
352 r__ = F77_FUNC(slapy2,SLAPY2)(&g, &c_b10);
353 g = d__[m] - p + e[l - 1] / (g + ( (g>0) ? r__ : -r__ ));
361 for (i__ = m; i__ <= i__1; ++i__) {
364 F77_FUNC(slartg,SLARTG)(&g, &f, &c__, &s, &r__);
369 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
376 work[*n - 1 + i__] = s;
382 F77_FUNC(slasr,SLASR)("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
403 i__1 = lendsv - lsv + 1;
404 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
407 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
409 } else if (iscale == 2) {
410 i__1 = lendsv - lsv + 1;
411 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
414 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
422 for (i__ = 1; i__ <= i__1; ++i__) {
423 if (fabs(e[i__])>GMX_FLOAT_MIN) {
432 F77_FUNC(slasrt,SLASRT)("I", n, &d__[1], info);
437 for (ii = 2; ii <= i__1; ++ii) {
442 for (j = ii; j <= i__2; ++j) {
451 F77_FUNC(sswap,SSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],