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 "gmx_lapack.h"
37 #include "lapack_limits.h"
39 #include <types/simple.h>
42 F77_FUNC(ssterf,SSTERF)(int *n,
54 float bb, rt1, rt2, eps, rte;
58 float gamma, alpha, sigma, anorm;
69 const float safmin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
89 ssfmax = sqrt(safmax) / 3.;
90 ssfmin = sqrt(safmin) / eps2;
100 F77_FUNC(slasrt,SLASRT)("I", n, &d__[1], info);
107 for (m = l1; m <= i__1; ++m) {
108 if (fabs(e[m]) <= sqrt(fabs(d__[m])) *
109 sqrt(fabs(d__[m + 1])) * eps) {
127 anorm = F77_FUNC(slanst,SLANST)("I", &i__1, &d__[l], &e[l]);
129 if (anorm > ssfmax) {
132 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
135 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
137 } else if (anorm < ssfmin) {
140 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
143 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
148 for (i__ = l; i__ <= i__1; ++i__) {
150 e[i__] = d__1 * d__1;
153 if (fabs(d__[lend]) < fabs(d__[l])) {
163 for (m = l; m <= i__1; ++m) {
164 if (fabs(e[m]) <= eps2 * fabs(d__[m] * d__[m + 1])) {
181 F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
192 if (jtot == nmaxit) {
198 sigma = (d__[l + 1] - p) / (rte * 2.);
199 r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
200 sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
204 gamma = d__[m] - sigma;
208 for (i__ = m - 1; i__ >= i__1; --i__) {
212 e[i__ + 1] = s * r__;
219 gamma = c__ * (alpha - sigma) - s * oldgam;
220 d__[i__ + 1] = oldgam + (alpha - gamma);
221 if (fabs(c__)>GMX_FLOAT_MIN) {
222 p = gamma * gamma / c__;
229 d__[l] = sigma + gamma;
245 for (m = l; m >= i__1; --m) {
246 if (fabs(e[m - 1]) <= eps2 * fabs(d__[m] * d__[m - 1])) {
262 rte = sqrt(e[l - 1]);
263 F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
274 if (jtot == nmaxit) {
279 rte = sqrt(e[l - 1]);
280 sigma = (d__[l - 1] - p) / (rte * 2.);
281 r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
282 sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
286 gamma = d__[m] - sigma;
290 for (i__ = m; i__ <= i__1; ++i__) {
294 e[i__ - 1] = s * r__;
300 alpha = d__[i__ + 1];
301 gamma = c__ * (alpha - sigma) - s * oldgam;
302 d__[i__] = oldgam + (alpha - gamma);
303 if (fabs(c__)>GMX_FLOAT_MIN) {
304 p = gamma * gamma / c__;
311 d__[l] = sigma + gamma;
327 i__1 = lendsv - lsv + 1;
328 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
332 i__1 = lendsv - lsv + 1;
333 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
341 for (i__ = 1; i__ <= i__1; ++i__) {
342 if (fabs(e[i__])>GMX_FLOAT_MIN) {