2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012, 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(slasd4,SLASD4)(int *n,
61 float eta, phi, eps, tau, psi;
65 float temp, prew, sg2lb, sg2ub, temp1, temp2, dtiim, delsq,
74 float erretm, dtipsq, rhoinv;
84 *sigma = sqrt(d__[1] * d__[1] + *rho * z__[1] * z__[1]);
90 F77_FUNC(slasd5,SLASD5)(i__, &d__[1], &z__[1], &delta[1], rho, sigma, &work[1]);
104 temp1 = temp / (d__[*n] + sqrt(d__[*n] * d__[*n] + temp));
106 for (j = 1; j <= i__1; ++j) {
107 work[j] = d__[j] + d__[*n] + temp1;
108 delta[j] = d__[j] - d__[*n] - temp1;
113 for (j = 1; j <= i__1; ++j) {
114 psi += z__[j] * z__[j] / (delta[j] * work[j]);
118 w = c__ + z__[ii] * z__[ii] / (delta[ii] * work[ii]) + z__[*n] * z__[*
119 n] / (delta[*n] * work[*n]);
122 temp1 = sqrt(d__[*n] * d__[*n] + *rho);
123 temp = z__[*n - 1] * z__[*n - 1] / ((d__[*n - 1] + temp1) * (d__[*
124 n] - d__[*n - 1] + *rho / (d__[*n] + temp1))) + z__[*n] *
130 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
131 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*
133 b = z__[*n] * z__[*n] * delsq;
135 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
137 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
142 delsq = (d__[*n] - d__[*n - 1]) * (d__[*n] + d__[*n - 1]);
143 a = -c__ * delsq + z__[*n - 1] * z__[*n - 1] + z__[*n] * z__[*n];
144 b = z__[*n] * z__[*n] * delsq;
147 tau = b * 2. / (sqrt(a * a + b * 4. * c__) - a);
149 tau = (a + sqrt(a * a + b * 4. * c__)) / (c__ * 2.);
154 eta = tau / (d__[*n] + sqrt(d__[*n] * d__[*n] + tau));
156 *sigma = d__[*n] + eta;
158 for (j = 1; j <= i__1; ++j) {
159 delta[j] = d__[j] - d__[*i__] - eta;
160 work[j] = d__[j] + d__[*i__] + eta;
167 for (j = 1; j <= i__1; ++j) {
168 temp = z__[j] / (delta[j] * work[j]);
169 psi += z__[j] * temp;
173 erretm = fabs(erretm);
175 temp = z__[*n] / (delta[*n] * work[*n]);
176 phi = z__[*n] * temp;
178 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
181 w = rhoinv + phi + psi;
183 if (fabs(w) <= eps * erretm) {
188 dtnsq1 = work[*n - 1] * delta[*n - 1];
189 dtnsq = work[*n] * delta[*n];
190 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
191 a = (dtnsq + dtnsq1) * w - dtnsq * dtnsq1 * (dpsi + dphi);
192 b = dtnsq * dtnsq1 * w;
196 if ( fabs(c__)<GMX_FLOAT_MIN) {
197 eta = *rho - *sigma * *sigma;
198 } else if (a >= 0.) {
199 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
201 eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
205 eta = -w / (dpsi + dphi);
213 eta /= *sigma + sqrt(eta + *sigma * *sigma);
215 for (j = 1; j <= i__1; ++j) {
226 for (j = 1; j <= i__1; ++j) {
227 temp = z__[j] / (work[j] * delta[j]);
228 psi += z__[j] * temp;
232 erretm = fabs(erretm);
234 temp = z__[*n] / (work[*n] * delta[*n]);
235 phi = z__[*n] * temp;
237 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (dpsi
240 w = rhoinv + phi + psi;
244 for (niter = iter; niter <= 20; ++niter) {
246 if (fabs(w) <= eps * erretm) {
249 dtnsq1 = work[*n - 1] * delta[*n - 1];
250 dtnsq = work[*n] * delta[*n];
251 c__ = w - dtnsq1 * dpsi - dtnsq * dphi;
252 a = (dtnsq + dtnsq1) * w - dtnsq1 * dtnsq * (dpsi + dphi);
253 b = dtnsq1 * dtnsq * w;
255 eta = (a + sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
257 eta = b * 2. / (a - sqrt(fabs(a * a - b * 4. * c__)));
261 eta = -w / (dpsi + dphi);
269 eta /= *sigma + sqrt(eta + *sigma * *sigma);
271 for (j = 1; j <= i__1; ++j) {
282 for (j = 1; j <= i__1; ++j) {
283 temp = z__[j] / (work[j] * delta[j]);
284 psi += z__[j] * temp;
288 erretm = fabs(erretm);
290 temp = z__[*n] / (work[*n] * delta[*n]);
291 phi = z__[*n] * temp;
293 erretm = (-phi - psi) * 8. + erretm - phi + rhoinv + fabs(tau) * (
296 w = rhoinv + phi + psi;
307 delsq = (d__[ip1] - d__[*i__]) * (d__[ip1] + d__[*i__]);
309 temp = delsq2 / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + delsq2));
311 for (j = 1; j <= i__1; ++j) {
312 work[j] = d__[j] + d__[*i__] + temp;
313 delta[j] = d__[j] - d__[*i__] - temp;
318 for (j = 1; j <= i__1; ++j) {
319 psi += z__[j] * z__[j] / (work[j] * delta[j]);
324 for (j = *n; j >= i__1; --j) {
325 phi += z__[j] * z__[j] / (work[j] * delta[j]);
327 c__ = rhoinv + psi + phi;
328 w = c__ + z__[*i__] * z__[*i__] / (work[*i__] * delta[*i__]) + z__[
329 ip1] * z__[ip1] / (work[ip1] * delta[ip1]);
336 a = c__ * delsq + z__[*i__] * z__[*i__] + z__[ip1] * z__[ip1];
337 b = z__[*i__] * z__[*i__] * delsq;
339 tau = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
341 tau = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
343 eta = tau / (d__[*i__] + sqrt(d__[*i__] * d__[*i__] + tau));
349 a = c__ * delsq - z__[*i__] * z__[*i__] - z__[ip1] * z__[ip1];
350 b = z__[ip1] * z__[ip1] * delsq;
352 tau = b * 2. / (a - sqrt(fabs(a * a + b * 4. * c__)));
354 tau = -(a + sqrt(fabs(a * a + b * 4. * c__))) / (c__ * 2.);
356 eta = tau / (d__[ip1] + sqrt(fabs(d__[ip1] * d__[ip1] + tau)));
361 *sigma = d__[*i__] + eta;
363 for (j = 1; j <= i__1; ++j) {
364 work[j] = d__[j] + d__[*i__] + eta;
365 delta[j] = d__[j] - d__[*i__] - eta;
369 *sigma = d__[ip1] + eta;
371 for (j = 1; j <= i__1; ++j) {
372 work[j] = d__[j] + d__[ip1] + eta;
373 delta[j] = d__[j] - d__[ip1] - eta;
383 for (j = 1; j <= i__1; ++j) {
384 temp = z__[j] / (work[j] * delta[j]);
385 psi += z__[j] * temp;
389 erretm = fabs(erretm);
394 for (j = *n; j >= i__1; --j) {
395 temp = z__[j] / (work[j] * delta[j]);
396 phi += z__[j] * temp;
401 w = rhoinv + phi + psi;
413 if (ii == 1 || ii == *n) {
417 temp = z__[ii] / (work[ii] * delta[ii]);
418 dw = dpsi + dphi + temp * temp;
419 temp = z__[ii] * temp;
421 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
424 if (fabs(w) <= eps * erretm) {
429 sg2lb = (sg2lb > tau) ? sg2lb : tau;
431 sg2ub = (sg2ub < tau) ? sg2ub : tau;
436 dtipsq = work[ip1] * delta[ip1];
437 dtisq = work[*i__] * delta[*i__];
439 d__1 = z__[*i__] / dtisq;
440 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
442 d__1 = z__[ip1] / dtipsq;
443 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
445 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
446 b = dtipsq * dtisq * w;
447 if ( fabs(c__)<GMX_FLOAT_MIN) {
448 if ( fabs(a)<GMX_FLOAT_MIN) {
450 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq * (dpsi +
453 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (dpsi +
458 } else if (a <= 0.) {
459 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
461 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
465 dtiim = work[iim1] * delta[iim1];
466 dtiip = work[iip1] * delta[iip1];
467 temp = rhoinv + psi + phi;
469 temp1 = z__[iim1] / dtiim;
471 c__ = temp - dtiip * (dpsi + dphi) - (d__[iim1] - d__[iip1]) *
472 (d__[iim1] + d__[iip1]) * temp1;
473 zz[0] = z__[iim1] * z__[iim1];
475 zz[2] = dtiip * dtiip * dphi;
477 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
480 temp1 = z__[iip1] / dtiip;
482 c__ = temp - dtiim * (dpsi + dphi) - (d__[iip1] - d__[iim1]) *
483 (d__[iim1] + d__[iip1]) * temp1;
485 zz[0] = dtiim * dtiim * dpsi;
487 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
489 zz[2] = z__[iip1] * z__[iip1];
491 zz[1] = z__[ii] * z__[ii];
493 dd[1] = delta[ii] * work[ii];
495 F77_FUNC(slaed6,SLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
505 temp1 = work[*i__] * delta[*i__];
508 temp1 = work[ip1] * delta[ip1];
511 if (temp > sg2ub || temp < sg2lb) {
513 eta = (sg2ub - tau) / 2.;
515 eta = (sg2lb - tau) / 2.;
520 eta /= *sigma + sqrt(*sigma * *sigma + eta);
526 for (j = 1; j <= i__1; ++j) {
535 for (j = 1; j <= i__1; ++j) {
536 temp = z__[j] / (work[j] * delta[j]);
537 psi += z__[j] * temp;
541 erretm = fabs(erretm);
546 for (j = *n; j >= i__1; --j) {
547 temp = z__[j] / (work[j] * delta[j]);
548 phi += z__[j] * temp;
553 temp = z__[ii] / (work[ii] * delta[ii]);
554 dw = dpsi + dphi + temp * temp;
555 temp = z__[ii] * temp;
556 w = rhoinv + phi + psi + temp;
557 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3. +
561 sg2lb = (sg2lb > tau) ? sg2lb : tau;
563 sg2ub = (sg2ub < tau) ? sg2ub : tau;
568 if (-w > fabs(prew) / 10.) {
572 if (w > fabs(prew) / 10.) {
579 for (niter = iter; niter <= 20; ++niter) {
581 if (fabs(w) <= eps * erretm) {
586 dtipsq = work[ip1] * delta[ip1];
587 dtisq = work[*i__] * delta[*i__];
590 d__1 = z__[*i__] / dtisq;
591 c__ = w - dtipsq * dw + delsq * (d__1 * d__1);
593 d__1 = z__[ip1] / dtipsq;
594 c__ = w - dtisq * dw - delsq * (d__1 * d__1);
597 temp = z__[ii] / (work[ii] * delta[ii]);
603 c__ = w - dtisq * dpsi - dtipsq * dphi;
605 a = (dtipsq + dtisq) * w - dtipsq * dtisq * dw;
606 b = dtipsq * dtisq * w;
607 if (fabs(c__)<GMX_FLOAT_MIN) {
608 if (fabs(a)<GMX_FLOAT_MIN) {
611 a = z__[*i__] * z__[*i__] + dtipsq * dtipsq *
614 a = z__[ip1] * z__[ip1] + dtisq * dtisq * (
618 a = dtisq * dtisq * dpsi + dtipsq * dtipsq * dphi;
622 } else if (a <= 0.) {
623 eta = (a - sqrt(fabs(a * a - b * 4. * c__))) / (c__ * 2.);
625 eta = b * 2. / (a + sqrt(fabs(a * a - b * 4. * c__)));
629 dtiim = work[iim1] * delta[iim1];
630 dtiip = work[iip1] * delta[iip1];
631 temp = rhoinv + psi + phi;
633 c__ = temp - dtiim * dpsi - dtiip * dphi;
634 zz[0] = dtiim * dtiim * dpsi;
635 zz[2] = dtiip * dtiip * dphi;
638 temp1 = z__[iim1] / dtiim;
640 temp2 = (d__[iim1] - d__[iip1]) * (d__[iim1] + d__[
642 c__ = temp - dtiip * (dpsi + dphi) - temp2;
643 zz[0] = z__[iim1] * z__[iim1];
645 zz[2] = dtiip * dtiip * dphi;
647 zz[2] = dtiip * dtiip * (dpsi - temp1 + dphi);
650 temp1 = z__[iip1] / dtiip;
652 temp2 = (d__[iip1] - d__[iim1]) * (d__[iim1] + d__[
654 c__ = temp - dtiim * (dpsi + dphi) - temp2;
656 zz[0] = dtiim * dtiim * dpsi;
658 zz[0] = dtiim * dtiim * (dpsi + (dphi - temp1));
660 zz[2] = z__[iip1] * z__[iip1];
664 dd[1] = delta[ii] * work[ii];
666 F77_FUNC(slaed6,SLAED6)(&niter, &orgati, &c__, dd, zz, &w, &eta, info);
676 temp1 = work[*i__] * delta[*i__];
679 temp1 = work[ip1] * delta[ip1];
682 if (temp > sg2ub || temp < sg2lb) {
684 eta = (sg2ub - tau) / 2.;
686 eta = (sg2lb - tau) / 2.;
691 eta /= *sigma + sqrt(*sigma * *sigma + eta);
695 for (j = 1; j <= i__1; ++j) {
706 for (j = 1; j <= i__1; ++j) {
707 temp = z__[j] / (work[j] * delta[j]);
708 psi += z__[j] * temp;
712 erretm = fabs(erretm);
717 for (j = *n; j >= i__1; --j) {
718 temp = z__[j] / (work[j] * delta[j]);
719 phi += z__[j] * temp;
724 temp = z__[ii] / (work[ii] * delta[ii]);
725 dw = dpsi + dphi + temp * temp;
726 temp = z__[ii] * temp;
727 w = rhoinv + phi + psi + temp;
728 erretm = (phi - psi) * 8. + erretm + rhoinv * 2. + fabs(temp) * 3.
730 if (w * prew > 0. && fabs(w) > fabs(prew) / 10.) {
735 sg2lb = (sg2lb > tau) ? sg2lb : tau;
737 sg2ub = (sg2ub < tau) ? sg2ub : tau;