src/gromacs/linearalgebra/gmx_lapack/slarrfx.c

   1 #include <math.h>
   2
   3 #include "gromacs/utility/real.h"
   4
   5 #include "../gmx_blas.h"
   6 #include "../gmx_lapack.h"
   7 #include "lapack_limits.h"
   8
   9
  10 void
  11 F77_FUNC(slarrfx,SLARRFX)(int *n,
  12         float *d__,
  13         float *l,
  14         float *ld,
  15         float *lld,
  16         int *ifirst,
  17         int *ilast,
  18         float *w,
  19         float *sigma,
  20         float *dplus,
  21         float *lplus,
  22         float *work,
  23         int *info)
  24 {
  25     int i1 = 1;
  26     int i__1;
  27     float d__2, d__3;
  28
  29     int i__;
  30     float s, eps, tmp, dmax1, dmax2, delta;
  31     --work;
  32     --lplus;
  33     --dplus;
  34     --w;
  35     --lld;
  36     --ld;
  37     --l;
  38     --d__;
  39     *info = 0;
  40     eps = GMX_FLOAT_EPS;
  41     *sigma = w[*ifirst];
  42     delta = eps * 2.;
  43
  44 L10:
  45     s = -(*sigma);
  46     dplus[1] = d__[1] + s;
  47     dmax1 = fabs(dplus[1]);
  48     i__1 = *n - 1;
  49     for (i__ = 1; i__ <= i__1; ++i__) {
  50         lplus[i__] = ld[i__] / dplus[i__];
  51         s = s * lplus[i__] * l[i__] - *sigma;
  52         dplus[i__ + 1] = d__[i__ + 1] + s;
  53         d__2 = dmax1, d__3 = fabs(dplus[i__ + 1]);
  54         dmax1 = (d__2>d__3) ? d__2 : d__3;
  55     }
  56     if (! (dmax1 > 0. || dmax1 < 1.)) {
  57         *sigma -= fabs(*sigma) * delta;
  58         delta *= 2.;
  59         goto L10;
  60     }
  61
  62     tmp = w[*ilast];
  63     delta = eps * 2.;
  64 L30:
  65     s = -tmp;
  66     work[1] = d__[1] + s;
  67     dmax2 = fabs(work[1]);
  68     i__1 = *n - 1;
  69     for (i__ = 1; i__ <= i__1; ++i__) {
  70         work[*n + i__] = ld[i__] / work[i__];
  71         s = s * work[*n + i__] * l[i__] - tmp;
  72         work[i__ + 1] = d__[i__ + 1] + s;
  73         d__2 = dmax2, d__3 = fabs(work[i__ + 1]);
  74         dmax2 = (d__2>d__3) ? d__2 : d__3;
  75     }
  76     if (! (dmax2 > 0. || dmax2 < 1.)) {
  77         tmp += fabs(tmp) * delta;
  78         delta *= 2.;
  79         goto L30;
  80     }
  81     if (dmax2 < dmax1) {
  82         *sigma = tmp;
  83         F77_FUNC(scopy,SCOPY)(n, &work[1], &i1, &dplus[1], &i1);
  84         i__1 = *n - 1;
  85         F77_FUNC(scopy,SCOPY)(&i__1, &work[*n + 1], &i1, &lplus[1], &i1);
  86     }
  87
  88     return;
  89 }