src/gromacs/linearalgebra/gmx_lapack/slagts.c

   1 #include <math.h>
   2 #include "gromacs/utility/real.h"
   3
   4 #include "../gmx_lapack.h"
   5 #include "lapack_limits.h"
   6
   7
   8 void
   9 F77_FUNC(slagts,SLAGTS)(int *job,
  10         int *n,
  11         float *a,
  12         float *b,
  13         float *c__,
  14         float *d__,
  15         int *in,
  16         float *y,
  17         float *tol,
  18         int *info)
  19 {
  20     int i__1;
  21     float d__1, d__2, d__4, d__5;
  22
  23     int k;
  24     float ak, eps, temp, pert, absak, sfmin;
  25     float bignum,minval;
  26     --y;
  27     --in;
  28     --d__;
  29     --c__;
  30     --b;
  31     --a;
  32
  33     *info = 0;
  34     if (fabs(*job) > 2 || *job == 0) {
  35         *info = -1;
  36     } else if (*n < 0) {
  37         *info = -2;
  38     }
  39     if (*info != 0) {
  40         i__1 = -(*info);
  41         return;
  42     }
  43
  44     if (*n == 0) {
  45         return;
  46     }
  47     eps = GMX_FLOAT_EPS;
  48     minval = GMX_FLOAT_MIN;
  49     sfmin = minval / eps;
  50
  51     bignum = 1. / sfmin;
  52
  53     if (*job < 0) {
  54         if (*tol <= 0.) {
  55             *tol = fabs(a[1]);
  56             if (*n > 1) {
  57                 d__1 = *tol;
  58                 d__2 = fabs(a[2]);
  59                 d__1 = (d__1>d__2) ? d__1 : d__2;
  60                 d__2 = fabs(b[1]);
  61                 *tol = (d__1>d__2) ? d__1 : d__2;
  62             }
  63             i__1 = *n;
  64             for (k = 3; k <= i__1; ++k) {
  65               d__4 = *tol;
  66               d__5 = fabs(a[k]);
  67               d__4 = (d__4>d__5) ? d__4 : d__5;
  68               d__5 = fabs(b[k - 1]);
  69               d__4 = (d__4>d__5) ? d__4 : d__5;
  70               d__5 = fabs(d__[k - 2]);
  71               *tol = (d__4>d__5) ? d__4 : d__5;
  72             }
  73             *tol *= eps;
  74             if (fabs(*tol)<GMX_FLOAT_MIN) {
  75                 *tol = eps;
  76             }
  77         }
  78     }
  79
  80     if (fabs(fabs(*job)-1.0)<GMX_FLOAT_MIN) {
  81         i__1 = *n;
  82         for (k = 2; k <= i__1; ++k) {
  83             if (in[k - 1] == 0) {
  84                 y[k] -= c__[k - 1] * y[k - 1];
  85             } else {
  86                 temp = y[k - 1];
  87                 y[k - 1] = y[k];
  88                 y[k] = temp - c__[k - 1] * y[k];
  89             }
  90         }
  91         if (*job == 1) {
  92             for (k = *n; k >= 1; --k) {
  93                 if (k <= *n - 2) {
  94                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
  95                 } else if (k == *n - 1) {
  96                     temp = y[k] - b[k] * y[k + 1];
  97                 } else {
  98                     temp = y[k];
  99                 }
 100                 ak = a[k];
 101                 absak = fabs(ak);
 102                 if (absak < 1.) {
 103                     if (absak < sfmin) {
 104                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 105                             *info = k;
 106                             return;
 107                         } else {
 108                             temp *= bignum;
 109                             ak *= bignum;
 110                         }
 111                     } else if (fabs(temp) > absak * bignum) {
 112                         *info = k;
 113                         return;
 114                     }
 115                 }
 116                 y[k] = temp / ak;
 117             }
 118         } else {
 119             for (k = *n; k >= 1; --k) {
 120                 if (k <= *n - 2) {
 121                     temp = y[k] - b[k] * y[k + 1] - d__[k] * y[k + 2];
 122                 } else if (k == *n - 1) {
 123                     temp = y[k] - b[k] * y[k + 1];
 124                 } else {
 125                     temp = y[k];
 126                 }
 127                 ak = a[k];
 128
 129                 pert = *tol;
 130                 if(ak<0)
 131                   pert *= -1.0;
 132 L40:
 133                 absak = fabs(ak);
 134                 if (absak < 1.) {
 135                     if (absak < sfmin) {
 136                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 137                             ak += pert;
 138                             pert *= 2;
 139                             goto L40;
 140                         } else {
 141                             temp *= bignum;
 142                             ak *= bignum;
 143                         }
 144                     } else if (fabs(temp) > absak * bignum) {
 145                         ak += pert;
 146                         pert *= 2;
 147                         goto L40;
 148                     }
 149                 }
 150                 y[k] = temp / ak;
 151             }
 152         }
 153     } else {
 154
 155         if (*job == 2) {
 156             i__1 = *n;
 157             for (k = 1; k <= i__1; ++k) {
 158                 if (k >= 3) {
 159                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 160                 } else if (k == 2) {
 161                     temp = y[k] - b[k - 1] * y[k - 1];
 162                 } else {
 163                     temp = y[k];
 164                 }
 165                 ak = a[k];
 166                 absak = fabs(ak);
 167                 if (absak < 1.) {
 168                     if (absak < sfmin) {
 169                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 170                             *info = k;
 171                             return;
 172                         } else {
 173                             temp *= bignum;
 174                             ak *= bignum;
 175                         }
 176                     } else if (fabs(temp) > absak * bignum) {
 177                         *info = k;
 178                         return;
 179                     }
 180                 }
 181                 y[k] = temp / ak;
 182             }
 183         } else {
 184             i__1 = *n;
 185             for (k = 1; k <= i__1; ++k) {
 186                 if (k >= 3) {
 187                     temp = y[k] - b[k - 1] * y[k - 1] - d__[k - 2] * y[k - 2];
 188                 } else if (k == 2) {
 189                     temp = y[k] - b[k - 1] * y[k - 1];
 190                 } else {
 191                     temp = y[k];
 192                 }
 193                 ak = a[k];
 194
 195                 pert = *tol;
 196                 if(ak<0)
 197                   pert *= -1.0;
 198
 199 L70:
 200                 absak = fabs(ak);
 201                 if (absak < 1.) {
 202                     if (absak < sfmin) {
 203                         if (fabs(absak)<GMX_FLOAT_MIN || fabs(temp) * sfmin > absak) {
 204                             ak += pert;
 205                             pert *= 2;
 206                             goto L70;
 207                         } else {
 208                             temp *= bignum;
 209                             ak *= bignum;
 210                         }
 211                     } else if (fabs(temp) > absak * bignum) {
 212                         ak += pert;
 213                         pert *= 2;
 214                         goto L70;
 215                     }
 216                 }
 217                 y[k] = temp / ak;
 218             }
 219         }
 220
 221         for (k = *n; k >= 2; --k) {
 222             if (in[k - 1] == 0) {
 223                 y[k - 1] -= c__[k - 1] * y[k];
 224             } else {
 225                 temp = y[k - 1];
 226                 y[k - 1] = y[k];
 227                 y[k] = temp - c__[k - 1] * y[k];
 228             }
 229         }
 230     }
 231
 232     return;
 233 }
 234
 235