src/gmxlib/gmx_lapack/slasd5.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 #include <math.h>
  36 #include "gmx_lapack.h"
  37
  38 void
  39 F77_FUNC(slasd5,SLASD5)(int *i__,
  40         float *d__,
  41         float *z__,
  42         float *delta,
  43         float *rho,
  44         float *dsigma,
  45         float *work)
  46 {
  47     float b, c__, w, del, tau, delsq;
  48
  49     --work;
  50     --delta;
  51     --z__;
  52     --d__;
  53
  54     del = d__[2] - d__[1];
  55     delsq = del * (d__[2] + d__[1]);
  56     if (*i__ == 1) {
  57         w = *rho * 4. * (z__[2] * z__[2] / (d__[1] + d__[2] * 3.) - z__[1] *
  58                 z__[1] / (d__[1] * 3. + d__[2])) / del + 1.;
  59         if (w > 0.) {
  60             b = delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  61             c__ = *rho * z__[1] * z__[1] * delsq;
  62
  63             tau = c__ * 2. / (b + sqrt(fabs(b * b - c__ * 4.)));
  64
  65             tau /= d__[1] + sqrt(d__[1] * d__[1] + tau);
  66             *dsigma = d__[1] + tau;
  67             delta[1] = -tau;
  68             delta[2] = del - tau;
  69             work[1] = d__[1] * 2. + tau;
  70             work[2] = d__[1] + tau + d__[2];
  71         } else {
  72             b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  73             c__ = *rho * z__[2] * z__[2] * delsq;
  74
  75             if (b > 0.) {
  76                 tau = c__ * -2. / (b + sqrt(b * b + c__ * 4.));
  77             } else {
  78                 tau = (b - sqrt(b * b + c__ * 4.)) / 2.;
  79             }
  80
  81             tau /= d__[2] + sqrt(fabs(d__[2] * d__[2] + tau));
  82             *dsigma = d__[2] + tau;
  83             delta[1] = -(del + tau);
  84             delta[2] = -tau;
  85             work[1] = d__[1] + tau + d__[2];
  86             work[2] = d__[2] * 2. + tau;
  87         }
  88     } else {
  89
  90         b = -delsq + *rho * (z__[1] * z__[1] + z__[2] * z__[2]);
  91         c__ = *rho * z__[2] * z__[2] * delsq;
  92
  93         if (b > 0.) {
  94             tau = (b + sqrt(b * b + c__ * 4.)) / 2.;
  95         } else {
  96             tau = c__ * 2. / (-b + sqrt(b * b + c__ * 4.));
  97         }
  98         tau /= d__[2] + sqrt(d__[2] * d__[2] + tau);
  99         *dsigma = d__[2] + tau;
 100         delta[1] = -(del + tau);
 101         delta[2] = -tau;
 102         work[1] = d__[1] + tau + d__[2];
 103         work[2] = d__[2] * 2. + tau;
 104     }
 105     return;
 106
 107 }