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.
37 #include "gmx_lapack.h"
39 #include <types/simple.h>
42 F77_FUNC(slaed6,SLAED6)(int *kniter,
52 float r__1, r__2, r__3, r__4;
56 float fc, df, ddf, eta, eps, base;
58 float temp, temp1, temp2, temp3, temp4;
61 float small1, small2, sminv1, sminv2, dscale[3], sclfac;
62 float zscale[3], erretm;
75 temp = (d__[3] - d__[2]) / 2.f;
76 c__ = *rho + z__[1] / (d__[1] - d__[2] - temp);
77 a = c__ * (d__[2] + d__[3]) + z__[2] + z__[3];
78 b = c__ * d__[2] * d__[3] + z__[2] * d__[3] + z__[3] * d__[2];
80 temp = (d__[1] - d__[2]) / 2.f;
81 c__ = *rho + z__[3] / (d__[3] - d__[2] - temp);
82 a = c__ * (d__[1] + d__[2]) + z__[1] + z__[2];
83 b = c__ * d__[1] * d__[2] + z__[1] * d__[2] + z__[2] * d__[1];
85 r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
86 temp = (r__1>r__2) ? r__1 : r__2;
92 } else if (a <= 0.f) {
93 *tau = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / (
96 *tau = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1))));
99 temp = *rho + z__[1] / (d__[1] - *tau) + z__[2] / (d__[2] - *tau) +
100 z__[3] / (d__[3] - *tau);
101 if (fabs(*finit) <= fabs(temp)) {
108 safemin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
109 i__1 = (int) (log(safemin) / log(base) / 3.f);
110 small1 = pow(base, i__1);
111 sminv1 = 1.f / small1;
112 small2 = small1 * small1;
113 sminv2 = sminv1 * sminv1;
116 r__3 = (r__1 = d__[2] - *tau, fabs(r__1)), r__4 = (r__2 = d__[3] - *
118 temp = (r__3<r__4) ? r__3 : r__4;
120 r__3 = (r__1 = d__[1] - *tau, fabs(r__1)), r__4 = (r__2 = d__[2] - *
122 temp = (r__3<r__4) ? r__3 : r__4;
125 if (temp <= small1) {
127 if (temp <= small2) {
138 for (i__ = 1; i__ <= 3; ++i__) {
139 dscale[i__ - 1] = d__[i__] * sclfac;
140 zscale[i__ - 1] = z__[i__] * sclfac;
145 for (i__ = 1; i__ <= 3; ++i__) {
146 dscale[i__ - 1] = d__[i__];
147 zscale[i__ - 1] = z__[i__];
153 for (i__ = 1; i__ <= 3; ++i__) {
154 temp = 1.f / (dscale[i__ - 1] - *tau);
155 temp1 = zscale[i__ - 1] * temp;
156 temp2 = temp1 * temp;
157 temp3 = temp2 * temp;
158 fc += temp1 / dscale[i__ - 1];
162 f = *finit + *tau * fc;
164 if (fabs(f) <= 0.f) {
168 for (niter = iter; niter <= 20; ++niter) {
170 temp1 = dscale[1] - *tau;
171 temp2 = dscale[2] - *tau;
173 temp1 = dscale[0] - *tau;
174 temp2 = dscale[1] - *tau;
176 a = (temp1 + temp2) * f - temp1 * temp2 * df;
177 b = temp1 * temp2 * f;
178 c__ = f - (temp1 + temp2) * df + temp1 * temp2 * ddf;
179 r__1 = fabs(a), r__2 = fabs(b), r__1 = ((r__1>r__2)? r__1:r__2), r__2 = fabs(c__);
180 temp = (r__1>r__2) ? r__1 : r__2;
186 } else if (a <= 0.f) {
187 eta = (a - sqrt((r__1 = a * a - b * 4.f * c__, fabs(r__1)))) / ( c__ * 2.f);
189 eta = b * 2.f / (a + sqrt((r__1 = a * a - b * 4.f * c__, fabs( r__1))));
191 if (f * eta >= 0.f) {
196 if (eta > 0.f && temp >= dscale[2]) {
197 eta = (dscale[2] - *tau) / 2.f;
200 if (eta < 0.f && temp <= dscale[1]) {
201 eta = (dscale[1] - *tau) / 2.f;
204 if (eta > 0.f && temp >= dscale[1]) {
205 eta = (dscale[1] - *tau) / 2.f;
207 if (eta < 0.f && temp <= dscale[0]) {
208 eta = (dscale[0] - *tau) / 2.f;
216 for (i__ = 1; i__ <= 3; ++i__) {
217 temp = 1.f / (dscale[i__ - 1] - *tau);
218 temp1 = zscale[i__ - 1] * temp;
219 temp2 = temp1 * temp;
220 temp3 = temp2 * temp;
221 temp4 = temp1 / dscale[i__ - 1];
223 erretm += fabs(temp4);
227 f = *finit + *tau * fc;
228 erretm = (fabs(*finit) + fabs(*tau) * erretm) * 8.f + fabs(*tau) * df;
229 if (fabs(f) <= eps * erretm) {