src/gromacs/linearalgebra/gmx_lapack/dlar1vx.cpp

   1 #include <cmath>
   2
   3 #include "gromacs/utility/real.h"
   4
   5 #include "../gmx_lapack.h"
   6 #include "lapack_limits.h"
   7
   8 void F77_FUNC(dlar1vx,DLAR1VX)(int *n,
   9               int *b1,
  10               int *bn,
  11               double *sigma,
  12               double *d__,
  13               double *l,
  14               double *ld,
  15               double *lld,
  16               double *eval,
  17               double *gersch,
  18               double *z__,
  19               double *ztz,
  20               double *mingma,
  21               int *r__,
  22               int *isuppz,
  23               double *work)
  24 {
  25     int i__1;
  26
  27     int i__, j;
  28     double s;
  29     int r1, r2;
  30     int to;
  31     double eps, tmp;
  32     int indp, inds, from;
  33     double dplus;
  34     int sawnan;
  35     int indumn;
  36     double dminus;
  37
  38     --work;
  39     --isuppz;
  40     --z__;
  41     --gersch;
  42     --lld;
  43     --ld;
  44     --l;
  45     --d__;
  46
  47     /* Function Body */
  48     eps = GMX_DOUBLE_EPS;
  49     if (*r__ == 0) {
  50
  51         r1 = *b1;
  52         r2 = *bn;
  53         i__1 = *bn;
  54         for (i__ = *b1; i__ <= i__1; ++i__) {
  55             if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
  56                 r1 = i__;
  57                 goto L20;
  58             }
  59         }
  60         goto L40;
  61 L20:
  62         i__1 = *b1;
  63         for (i__ = *bn; i__ >= i__1; --i__) {
  64             if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
  65                 r2 = i__;
  66                 goto L40;
  67             }
  68         }
  69     } else {
  70         r1 = *r__;
  71         r2 = *r__;
  72     }
  73
  74 L40:
  75     indumn = *n;
  76     inds = (*n << 1) + 1;
  77     indp = *n * 3 + 1;
  78     sawnan = 0;
  79
  80     if (*b1 == 1) {
  81         work[inds] = 0.;
  82     } else {
  83         work[inds] = lld[*b1 - 1];
  84     }
  85     s = work[inds] - *sigma;
  86     i__1 = r2 - 1;
  87     for (i__ = *b1; i__ <= i__1; ++i__) {
  88         dplus = d__[i__] + s;
  89         work[i__] = ld[i__] / dplus;
  90         work[inds + i__] = s * work[i__] * l[i__];
  91         s = work[inds + i__] - *sigma;
  92     }
  93
  94     if (std::isnan(s)) {
  95
  96         sawnan = 1;
  97         j = *b1 + 1;
  98 L60:
  99     if (!std::isnan(work[inds + j])) {
 100             ++j;
 101             goto L60;
 102         }
 103         work[inds + j] = lld[j];
 104         s = work[inds + j] - *sigma;
 105         i__1 = r2 - 1;
 106         for (i__ = j + 1; i__ <= i__1; ++i__) {
 107             dplus = d__[i__] + s;
 108             work[i__] = ld[i__] / dplus;
 109             if (std::abs(work[i__])<GMX_DOUBLE_MIN) {
 110                 work[inds + i__] = lld[i__];
 111             } else {
 112                 work[inds + i__] = s * work[i__] * l[i__];
 113             }
 114             s = work[inds + i__] - *sigma;
 115         }
 116     }
 117
 118     work[indp + *bn - 1] = d__[*bn] - *sigma;
 119     i__1 = r1;
 120     for (i__ = *bn - 1; i__ >= i__1; --i__) {
 121         dminus = lld[i__] + work[indp + i__];
 122         tmp = d__[i__] / dminus;
 123         work[indumn + i__] = l[i__] * tmp;
 124         work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
 125     }
 126     tmp = work[indp + r1 - 1];
 127     if (std::isnan(tmp)) {
 128
 129         sawnan = 1;
 130         j = *bn - 3;
 131 L90:
 132     if (!std::isnan(work[indp + j])) {
 133             --j;
 134             goto L90;
 135         }
 136         work[indp + j] = d__[j + 1] - *sigma;
 137         i__1 = r1;
 138         for (i__ = j; i__ >= i__1; --i__) {
 139             dminus = lld[i__] + work[indp + i__];
 140             tmp = d__[i__] / dminus;
 141             work[indumn + i__] = l[i__] * tmp;
 142             if (std::abs(tmp)<GMX_DOUBLE_MIN) {
 143                 work[indp + i__ - 1] = d__[i__] - *sigma;
 144             } else {
 145                 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
 146             }
 147         }
 148     }
 149
 150     *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
 151     if (std::abs(*mingma)<GMX_DOUBLE_MIN) {
 152         *mingma = eps * work[inds + r1 - 1];
 153     }
 154     *r__ = r1;
 155     i__1 = r2 - 1;
 156     for (i__ = r1; i__ <= i__1; ++i__) {
 157         tmp = work[inds + i__] + work[indp + i__];
 158         if (std::abs(tmp)<GMX_DOUBLE_MIN) {
 159             tmp = eps * work[inds + i__];
 160         }
 161         if (std::abs(tmp) < std::abs(*mingma)) {
 162             *mingma = tmp;
 163             *r__ = i__ + 1;
 164         }
 165     }
 166
 167     isuppz[1] = *b1;
 168     isuppz[2] = *bn;
 169     z__[*r__] = 1.;
 170     *ztz = 1.;
 171     if (! sawnan) {
 172         from = *r__ - 1;
 173         i__1 = *r__ - 32;
 174         to = (i__1>(*b1)) ? i__1 : (*b1);
 175 L120:
 176         if (from >= *b1) {
 177             i__1 = to;
 178             for (i__ = from; i__ >= i__1; --i__) {
 179                 z__[i__] = -(work[i__] * z__[i__ + 1]);
 180                 *ztz += z__[i__] * z__[i__];
 181             }
 182             if (std::abs(z__[to]) <= eps && std::abs(z__[to + 1]) <= eps) {
 183                 isuppz[1] = to + 2;
 184             } else {
 185                 from = to - 1;
 186                 i__1 = to - 32;
 187                 to = (i__1>*b1) ? i__1 : *b1;
 188                 goto L120;
 189             }
 190         }
 191         from = *r__ + 1;
 192         i__1 = *r__ + 32;
 193         to = (i__1<*bn) ? i__1 : *bn;
 194 L140:
 195         if (from <= *bn) {
 196             i__1 = to;
 197             for (i__ = from; i__ <= i__1; ++i__) {
 198                 z__[i__] = -(work[indumn + i__ - 1] * z__[i__ - 1]);
 199                 *ztz += z__[i__] * z__[i__];
 200             }
 201             if (std::abs(z__[to]) <= eps && std::abs(z__[to - 1]) <= eps) {
 202                 isuppz[2] = to - 2;
 203             } else {
 204                 from = to + 1;
 205                 i__1 = to + 32;
 206                 to = (i__1<*bn) ? i__1 : *bn;
 207                 goto L140;
 208             }
 209         }
 210     } else {
 211         i__1 = *b1;
 212         for (i__ = *r__ - 1; i__ >= i__1; --i__) {
 213             if (std::abs(z__[i__ + 1])<GMX_DOUBLE_MIN) {
 214                 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
 215             } else {
 216                 z__[i__] = -(work[i__] * z__[i__ + 1]);
 217             }
 218             if (std::abs(z__[i__]) <= eps && std::abs(z__[i__ + 1]) <= eps) {
 219                 isuppz[1] = i__ + 2;
 220                 goto L170;
 221             }
 222             *ztz += z__[i__] * z__[i__];
 223         }
 224 L170:
 225         i__1 = *bn - 1;
 226         for (i__ = *r__; i__ <= i__1; ++i__) {
 227             if (std::abs(z__[i__])<GMX_DOUBLE_MIN) {
 228                 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
 229             } else {
 230                 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
 231             }
 232             if (std::abs(z__[i__]) <= eps && std::abs(z__[i__ + 1]) <= eps) {
 233                 isuppz[2] = i__ - 1;
 234                 break;
 235             }
 236             *ztz += z__[i__ + 1] * z__[i__ + 1];
 237         }
 238     }
 239
 240     return;
 241
 242 }
 243
 244