src/gromacs/linearalgebra/gmx_lapack/dlagtf.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
   9 void
  10 F77_FUNC(dlagtf,DLAGTF)(int *n,
  11         double *a,
  12         double *lambda,
  13         double *b,
  14         double *c__,
  15         double *tol,
  16         double *d__,
  17         int *in,
  18         int *info)
  19 {
  20     int i__1;
  21
  22     int k;
  23     double tl, eps, piv1, piv2, temp, mult, scale1, scale2;
  24
  25     --in;
  26     --d__;
  27     --c__;
  28     --b;
  29     --a;
  30
  31     *info = 0;
  32     if (*n < 0) {
  33         *info = -1;
  34         i__1 = -(*info);
  35         return;
  36     }
  37
  38     if (*n == 0)
  39         return;
  40
  41     a[1] -= *lambda;
  42     in[*n] = 0;
  43     if (*n == 1) {
  44         if (fabs(a[1])<GMX_DOUBLE_MIN) {
  45             in[1] = 1;
  46         }
  47         return;
  48     }
  49
  50     eps = GMX_DOUBLE_EPS;
  51
  52     tl = (*tol>eps) ? *tol : eps;
  53     scale1 = fabs(a[1]) + fabs(b[1]);
  54     i__1 = *n - 1;
  55     for (k = 1; k <= i__1; ++k) {
  56         a[k + 1] -= *lambda;
  57         scale2 = fabs(c__[k]) + fabs(a[k + 1]);
  58         if (k < *n - 1) {
  59             scale2 += fabs(b[k + 1]);
  60         }
  61         if (fabs(a[k])<GMX_DOUBLE_MIN) {
  62             piv1 = 0.;
  63         } else {
  64             piv1 = fabs(a[k]) / scale1;
  65         }
  66         if (fabs(c__[k])<GMX_DOUBLE_MIN) {
  67             in[k] = 0;
  68             piv2 = 0.;
  69             scale1 = scale2;
  70             if (k < *n - 1) {
  71                 d__[k] = 0.;
  72             }
  73         } else {
  74             piv2 = fabs(c__[k]) / scale2;
  75             if (piv2 <= piv1) {
  76                 in[k] = 0;
  77                 scale1 = scale2;
  78                 c__[k] /= a[k];
  79                 a[k + 1] -= c__[k] * b[k];
  80                 if (k < *n - 1) {
  81                     d__[k] = 0.;
  82                 }
  83             } else {
  84                 in[k] = 1;
  85                 mult = a[k] / c__[k];
  86                 a[k] = c__[k];
  87                 temp = a[k + 1];
  88                 a[k + 1] = b[k] - mult * temp;
  89                 if (k < *n - 1) {
  90                     d__[k] = b[k + 1];
  91                     b[k + 1] = -mult * d__[k];
  92                 }
  93                 b[k] = temp;
  94                 c__[k] = mult;
  95             }
  96         }
  97         if (((piv1>piv2) ? piv1 : piv2) <= tl && in[*n] == 0) {
  98             in[*n] = k;
  99         }
 100     }
 101     if (fabs(a[*n]) <= scale1 * tl && in[*n] == 0) {
 102         in[*n] = *n;
 103     }
 104
 105     return;
 106
 107 }
 108
 109