src/gromacs/linearalgebra/gmx_lapack/dlasq5.c

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