src/gromacs/linearalgebra/gmx_lapack/slasv2.c

   1 #include <math.h>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "gromacs/utility/real.h"
   6
   7 void
   8 F77_FUNC(slasv2,SLASV2)(float *f,
   9                         float *g,
  10                         float *h__,
  11                         float *ssmin,
  12                         float *ssmax,
  13                         float *snr,
  14                         float *csr,
  15                         float *snl,
  16                         float *csl)
  17 {
  18     float d__1;
  19
  20     float a, d__, l, m, r__, s, t, fa, ga, ha, ft, gt, ht, mm, tt,
  21              clt, crt, slt, srt;
  22     int pmax;
  23     float temp;
  24     int swap;
  25     float tsign=1.0;
  26     int gasmal;
  27
  28     ft = *f;
  29     fa = fabs(ft);
  30     ht = *h__;
  31     ha = fabs(*h__);
  32
  33     pmax = 1;
  34     swap = ha > fa;
  35     if (swap) {
  36         pmax = 3;
  37         temp = ft;
  38         ft = ht;
  39         ht = temp;
  40         temp = fa;
  41         fa = ha;
  42         ha = temp;
  43
  44     }
  45     gt = *g;
  46     ga = fabs(gt);
  47     if (fabs(ga)<GMX_FLOAT_MIN) {
  48
  49         *ssmin = ha;
  50         *ssmax = fa;
  51         clt = 1.;
  52         crt = 1.;
  53         slt = 0.;
  54         srt = 0.;
  55     } else {
  56         gasmal = 1;
  57         if (ga > fa) {
  58             pmax = 2;
  59             if (fa / ga < GMX_FLOAT_EPS) {
  60
  61                 gasmal = 0;
  62                 *ssmax = ga;
  63                 if (ha > 1.) {
  64                     *ssmin = fa / (ga / ha);
  65                 } else {
  66                     *ssmin = fa / ga * ha;
  67                 }
  68                 clt = 1.;
  69                 slt = ht / gt;
  70                 srt = 1.;
  71                 crt = ft / gt;
  72             }
  73         }
  74         if (gasmal) {
  75
  76             d__ = fa - ha;
  77             if ( fabs( fa - d__ )<GMX_FLOAT_EPS*fabs( fa + d__ )) {
  78                 l = 1.;
  79             } else {
  80                 l = d__ / fa;
  81             }
  82
  83             m = gt / ft;
  84             t = 2. - l;
  85
  86             mm = m * m;
  87             tt = t * t;
  88             s = sqrt(tt + mm);
  89
  90             if ( fabs(l)<GMX_FLOAT_MIN) {
  91                 r__ = fabs(m);
  92             } else {
  93                 r__ = sqrt(l * l + mm);
  94             }
  95             a = (s + r__) * .5;
  96
  97             *ssmin = ha / a;
  98             *ssmax = fa * a;
  99             if ( fabs(mm)<GMX_FLOAT_MIN) {
 100
 101                 if (fabs(l)<GMX_FLOAT_MIN) {
 102                     t = ( (ft>0) ? 2.0 : -2.0) * ( (gt>0) ? 1.0 : -1.0);
 103                 } else {
 104                     t = gt / ( (ft>0) ? d__ : d__) + m / t;
 105                 }
 106             } else {
 107                 t = (m / (s + t) + m / (r__ + l)) * (a + 1.);
 108             }
 109             l = sqrt(t * t + 4.);
 110             crt = 2. / l;
 111             srt = t / l;
 112             clt = (crt + srt * m) / a;
 113             slt = ht / ft * srt / a;
 114         }
 115     }
 116     if (swap) {
 117         *csl = srt;
 118         *snl = crt;
 119         *csr = slt;
 120         *snr = clt;
 121     } else {
 122         *csl = clt;
 123         *snl = slt;
 124         *csr = crt;
 125         *snr = srt;
 126     }
 127
 128     if (pmax == 1) {
 129         tsign = ( (*csr>0) ? 1.0 : -1.0) * ( (*csl>0) ? 1.0 : -1.0) * ( (*f>0) ? 1.0 : -1.0);
 130     }
 131     if (pmax == 2) {
 132         tsign = ( (*snr>0) ? 1.0 : -1.0) * ( (*csl>0) ? 1.0 : -1.0) * ( (*g>0) ? 1.0 : -1.0);
 133     }
 134     if (pmax == 3) {
 135         tsign = ( (*snr>0) ? 1.0 : -1.0) * ( (*snl>0) ? 1.0 : -1.0) * ( (*h__>0) ? 1.0 : -1.0);
 136     }
 137     if(tsign<0)
 138       *ssmax *= -1.0;
 139     d__1 = tsign * ( (*f>0) ? 1.0 : -1.0) * ( (*h__>0) ? 1.0 : -1.0);
 140     if(d__1<0)
 141       *ssmin *= -1.0;
 142     return;
 143
 144 }