src/gromacs/linearalgebra/gmx_lapack/slasq6.cpp

   1 #include <cmath>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "gromacs/utility/real.h"
   6
   7 void
   8 F77_FUNC(slasq6,SLASQ6)(int *i0,
   9         int *n0,
  10         float *z__,
  11         int *pp,
  12         float *dmin__,
  13         float *dmin1,
  14         float *dmin2,
  15         float *dn,
  16         float *dnm1,
  17         float *dnm2)
  18 {
  19     int i__1;
  20     float d__1, d__2;
  21
  22     /* Local variables */
  23     float d__;
  24     int j4, j4p2;
  25     float emin, temp;
  26     const float safemin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
  27
  28     --z__;
  29
  30     if (*n0 - *i0 - 1 <= 0) {
  31         return;
  32     }
  33
  34     j4 = (*i0 << 2) + *pp - 3;
  35     emin = z__[j4 + 4];
  36     d__ = z__[j4];
  37     *dmin__ = d__;
  38
  39     if (*pp == 0) {
  40         i__1 = 4*(*n0 - 3);
  41         for (j4 = *i0*4; j4 <= i__1; j4 += 4) {
  42             z__[j4 - 2] = d__ + z__[j4 - 1];
  43             if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
  44                 z__[j4] = 0.;
  45                 d__ = z__[j4 + 1];
  46                 *dmin__ = d__;
  47                 emin = 0.;
  48             } else if (safemin * z__[j4 + 1] < z__[j4 - 2] && safemin * z__[j4
  49                     - 2] < z__[j4 + 1]) {
  50                 temp = z__[j4 + 1] / z__[j4 - 2];
  51                 z__[j4] = z__[j4 - 1] * temp;
  52                 d__ *= temp;
  53             } else {
  54                 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
  55                 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
  56             }
  57             if(d__<*dmin__)
  58               *dmin__ = d__;
  59
  60             d__1 = emin, d__2 = z__[j4];
  61             emin = (d__1<d__2) ? d__1 : d__2;
  62         }
  63     } else {
  64         i__1 = 4*(*n0 - 3);
  65         for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
  66             z__[j4 - 3] = d__ + z__[j4];
  67             if (std::abs(z__[j4 - 3])<GMX_FLOAT_MIN) {
  68                 z__[j4 - 1] = 0.;
  69                 d__ = z__[j4 + 2];
  70                 *dmin__ = d__;
  71                 emin = 0.;
  72             } else if (safemin * z__[j4 + 2] < z__[j4 - 3] && safemin * z__[j4
  73                     - 3] < z__[j4 + 2]) {
  74                 temp = z__[j4 + 2] / z__[j4 - 3];
  75                 z__[j4 - 1] = z__[j4] * temp;
  76                 d__ *= temp;
  77             } else {
  78                 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
  79                 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
  80             }
  81             if(d__<*dmin__)
  82               *dmin__ = d__;
  83             d__1 = emin, d__2 = z__[j4 - 1];
  84             emin = (d__1<d__2) ? d__1 : d__2;
  85         }
  86     }
  87
  88     *dnm2 = d__;
  89     *dmin2 = *dmin__;
  90     j4 = 4*(*n0 - 2) - *pp;
  91     j4p2 = j4 + (*pp << 1) - 1;
  92     z__[j4 - 2] = *dnm2 + z__[j4p2];
  93     if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
  94         z__[j4] = 0.;
  95         *dnm1 = z__[j4p2 + 2];
  96         *dmin__ = *dnm1;
  97         emin = 0.;
  98     } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
  99             z__[j4p2 + 2]) {
 100         temp = z__[j4p2 + 2] / z__[j4 - 2];
 101         z__[j4] = z__[j4p2] * temp;
 102         *dnm1 = *dnm2 * temp;
 103     } else {
 104         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 105         *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
 106     }
 107     if(*dnm1<*dmin__)
 108       *dmin__ = *dnm1;
 109
 110     *dmin1 = *dmin__;
 111     j4 += 4;
 112     j4p2 = j4 + (*pp << 1) - 1;
 113     z__[j4 - 2] = *dnm1 + z__[j4p2];
 114     if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
 115         z__[j4] = 0.;
 116         *dn = z__[j4p2 + 2];
 117         *dmin__ = *dn;
 118         emin = 0.;
 119     } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
 120             z__[j4p2 + 2]) {
 121         temp = z__[j4p2 + 2] / z__[j4 - 2];
 122         z__[j4] = z__[j4p2] * temp;
 123         *dn = *dnm1 * temp;
 124     } else {
 125         z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
 126         *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
 127     }
 128     if(*dn<*dmin__)
 129       *dmin__ = *dn;
 130
 131     z__[j4 + 2] = *dn;
 132     z__[(*n0 << 2) - *pp] = emin;
 133     return;
 134
 135
 136 }