src/gmxlib/gmx_lapack/slasv2.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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 #include "lapack_limits.h"
  38
  39 #include <types/simple.h>
  40
  41 void
  42 F77_FUNC(slasv2,SLASV2)(float *f,
  43                         float *g,
  44                         float *h__,
  45                         float *ssmin,
  46                         float *ssmax,
  47                         float *snr,
  48                         float *csr,
  49                         float *snl,
  50                         float *csl)
  51 {
  52     float d__1;
  53
  54     float a, d__, l, m, r__, s, t, fa, ga, ha, ft, gt, ht, mm, tt,
  55              clt, crt, slt, srt;
  56     int pmax;
  57     float temp;
  58     int swap;
  59     float tsign=1.0;
  60     int gasmal;
  61
  62     ft = *f;
  63     fa = fabs(ft);
  64     ht = *h__;
  65     ha = fabs(*h__);
  66
  67     pmax = 1;
  68     swap = ha > fa;
  69     if (swap) {
  70         pmax = 3;
  71         temp = ft;
  72         ft = ht;
  73         ht = temp;
  74         temp = fa;
  75         fa = ha;
  76         ha = temp;
  77
  78     }
  79     gt = *g;
  80     ga = fabs(gt);
  81     if (fabs(ga)<GMX_FLOAT_MIN) {
  82
  83         *ssmin = ha;
  84         *ssmax = fa;
  85         clt = 1.;
  86         crt = 1.;
  87         slt = 0.;
  88         srt = 0.;
  89     } else {
  90         gasmal = 1;
  91         if (ga > fa) {
  92             pmax = 2;
  93             if (fa / ga < GMX_FLOAT_EPS) {
  94
  95                 gasmal = 0;
  96                 *ssmax = ga;
  97                 if (ha > 1.) {
  98                     *ssmin = fa / (ga / ha);
  99                 } else {
 100                     *ssmin = fa / ga * ha;
 101                 }
 102                 clt = 1.;
 103                 slt = ht / gt;
 104                 srt = 1.;
 105                 crt = ft / gt;
 106             }
 107         }
 108         if (gasmal) {
 109
 110             d__ = fa - ha;
 111             if ( fabs( fa - d__ )<GMX_FLOAT_EPS*fabs( fa + d__ )) {
 112                 l = 1.;
 113             } else {
 114                 l = d__ / fa;
 115             }
 116
 117             m = gt / ft;
 118             t = 2. - l;
 119
 120             mm = m * m;
 121             tt = t * t;
 122             s = sqrt(tt + mm);
 123
 124             if ( fabs(l)<GMX_FLOAT_MIN) {
 125                 r__ = fabs(m);
 126             } else {
 127                 r__ = sqrt(l * l + mm);
 128             }
 129             a = (s + r__) * .5;
 130
 131             *ssmin = ha / a;
 132             *ssmax = fa * a;
 133             if ( fabs(mm)<GMX_FLOAT_MIN) {
 134
 135                 if (fabs(l)<GMX_FLOAT_MIN) {
 136                     t = ( (ft>0) ? 2.0 : -2.0) * ( (gt>0) ? 1.0 : -1.0);
 137                 } else {
 138                     t = gt / ( (ft>0) ? d__ : d__) + m / t;
 139                 }
 140             } else {
 141                 t = (m / (s + t) + m / (r__ + l)) * (a + 1.);
 142             }
 143             l = sqrt(t * t + 4.);
 144             crt = 2. / l;
 145             srt = t / l;
 146             clt = (crt + srt * m) / a;
 147             slt = ht / ft * srt / a;
 148         }
 149     }
 150     if (swap) {
 151         *csl = srt;
 152         *snl = crt;
 153         *csr = slt;
 154         *snr = clt;
 155     } else {
 156         *csl = clt;
 157         *snl = slt;
 158         *csr = crt;
 159         *snr = srt;
 160     }
 161
 162     if (pmax == 1) {
 163         tsign = ( (*csr>0) ? 1.0 : -1.0) * ( (*csl>0) ? 1.0 : -1.0) * ( (*f>0) ? 1.0 : -1.0);
 164     }
 165     if (pmax == 2) {
 166         tsign = ( (*snr>0) ? 1.0 : -1.0) * ( (*csl>0) ? 1.0 : -1.0) * ( (*g>0) ? 1.0 : -1.0);
 167     }
 168     if (pmax == 3) {
 169         tsign = ( (*snr>0) ? 1.0 : -1.0) * ( (*snl>0) ? 1.0 : -1.0) * ( (*h__>0) ? 1.0 : -1.0);
 170     }
 171     if(tsign<0)
 172       *ssmax *= -1.0;
 173     d__1 = tsign * ( (*f>0) ? 1.0 : -1.0) * ( (*h__>0) ? 1.0 : -1.0);
 174     if(d__1<0)
 175       *ssmin *= -1.0;
 176     return;
 177
 178 }