src/gromacs/linearalgebra/gmx_lapack/slarrbx.c

   1 #include <math.h>
   2
   3 #include "gromacs/utility/real.h"
   4
   5 #include "../gmx_lapack.h"
   6 #include "lapack_limits.h"
   7
   8 void
   9 F77_FUNC(slarrbx,SLARRBX)(int *n,
  10          float *d__,
  11          float *l,
  12          float *ld,
  13          float *lld,
  14          int *ifirst,
  15          int *ilast,
  16          float *rtol1,
  17          float *rtol2,
  18          int *offset,
  19          float *w,
  20          float *wgap,
  21          float *werr,
  22          float *work,
  23          int *iwork,
  24          int *info)
  25 {
  26     int i__1, i__2, i__3;
  27     float d__1, d__2;
  28
  29     int i__, j, k, p;
  30     float s;
  31     int i1, i2, ii, kk;
  32     float fac, gap, mid;
  33     int cnt;
  34     float eps, tmp, left;
  35     int nint, prev, next, nleft;
  36     float right, width, dplus, error;
  37     int nright, olnint;
  38     k = 0;
  39     right = 0.0;
  40
  41     --iwork;
  42     --work;
  43     --werr;
  44     --wgap;
  45     --w;
  46     --lld;
  47     --ld;
  48     --l;
  49     --d__;
  50
  51     *info = 0;
  52     eps = GMX_FLOAT_EPS;
  53     i__1 = *n << 1;
  54     for (i__ = 1; i__ <= i__1; ++i__) {
  55         iwork[i__] = 0;
  56     }
  57     i1 = *ifirst;
  58     i2 = *ifirst;
  59     prev = 0;
  60     i__1 = *ilast;
  61     for (i__ = *ifirst; i__ <= i__1; ++i__) {
  62         ii = i__ - *offset;
  63         if (i__ == *ifirst) {
  64             gap = wgap[ii];
  65         } else if (i__ == *ilast) {
  66             gap = wgap[ii - 1];
  67         } else {
  68             d__1 = wgap[ii - 1], d__2 = wgap[ii];
  69             gap = (d__1<d__2) ? d__1 : d__2;
  70         }
  71         error = werr[ii];
  72         k = i__ << 1;
  73         iwork[k - 1] = 1;
  74         i2 = i__;
  75     }
  76
  77     i__ = i1;
  78     nint = 0;
  79 L30:
  80     if (i__ <= i2) {
  81         ii = i__ - *offset;
  82         if (iwork[(i__ << 1) - 1] == 1) {
  83             fac = 1.;
  84             left = w[ii] - werr[ii];
  85
  86
  87 L40:
  88             if (i__ > i1 && left <= right) {
  89                 left = right;
  90                 cnt = i__ - 1;
  91             } else {
  92                 s = -left;
  93                 cnt = 0;
  94                 i__1 = *n - 1;
  95                 for (j = 1; j <= i__1; ++j) {
  96                     dplus = d__[j] + s;
  97                     s = s * lld[j] / dplus - left;
  98                     if (dplus < 0.) {
  99                         ++cnt;
 100                     }
 101                 }
 102                 dplus = d__[*n] + s;
 103                 if (dplus < 0.) {
 104                     ++cnt;
 105                 }
 106                 if (! (s > 0. || s < 1.)) {
 107
 108                     cnt = 0;
 109                     s = -left;
 110                     i__1 = *n - 1;
 111                     for (j = 1; j <= i__1; ++j) {
 112                         dplus = d__[j] + s;
 113                         if (dplus < 0.) {
 114                             ++cnt;
 115                         }
 116                         tmp = lld[j] / dplus;
 117                         if (fabs(tmp)<GMX_FLOAT_MIN) {
 118                             s = lld[j] - left;
 119                         } else {
 120                             s = s * tmp - left;
 121                         }
 122                     }
 123                     dplus = d__[*n] + s;
 124                     if (dplus < 0.) {
 125                         ++cnt;
 126                     }
 127                 }
 128                 if (cnt > i__ - 1) {
 129                     left -= werr[ii] * fac;
 130                     fac *= 2.;
 131                     goto L40;
 132                 }
 133             }
 134             nleft = cnt + 1;
 135             i1 = (i1<nleft) ? i1 : nleft;
 136             fac = 1.;
 137             right = w[ii] + werr[ii];
 138 L60:
 139             s = -right;
 140             cnt = 0;
 141             i__1 = *n - 1;
 142             for (j = 1; j <= i__1; ++j) {
 143                 dplus = d__[j] + s;
 144                 s = s * lld[j] / dplus - right;
 145                 if (dplus < 0.) {
 146                     ++cnt;
 147                 }
 148             }
 149             dplus = d__[*n] + s;
 150             if (dplus < 0.) {
 151                 ++cnt;
 152             }
 153             if (! (s > 0. || s < 1.)) {
 154
 155                 cnt = 0;
 156                 s = -right;
 157                 i__1 = *n - 1;
 158                 for (j = 1; j <= i__1; ++j) {
 159                     dplus = d__[j] + s;
 160                     if (dplus < 0.) {
 161                         ++cnt;
 162                     }
 163                     tmp = lld[j] / dplus;
 164                     if (fabs(tmp)<GMX_FLOAT_MIN) {
 165                         s = lld[j] - right;
 166                     } else {
 167                         s = s * tmp - right;
 168                     }
 169                 }
 170                 dplus = d__[*n] + s;
 171                 if (dplus < 0.) {
 172                     ++cnt;
 173                 }
 174             }
 175             if (cnt < i__) {
 176                 right += werr[ii] * fac;
 177                 fac *= 2.;
 178                 goto L60;
 179             }
 180             cnt = (cnt<i2) ? cnt : i2;
 181             ++nint;
 182             k = nleft << 1;
 183             work[k - 1] = left;
 184             work[k] = right;
 185             i__ = cnt + 1;
 186             iwork[k - 1] = i__;
 187             iwork[k] = cnt;
 188             if (prev != nleft - 1) {
 189                 work[k - 2] = left;
 190             }
 191             prev = nleft;
 192         } else {
 193             right = work[i__ * 2];
 194
 195             ++iwork[k - 1];
 196             prev = i__;
 197             ++i__;
 198         }
 199         goto L30;
 200     }
 201     if (i__ <= *n && iwork[(i__ << 1) - 1] != -1) {
 202         work[(i__ << 1) - 1] = work[prev * 2];
 203     }
 204
 205 L80:
 206     prev = i1 - 1;
 207     olnint = nint;
 208     i__ = i1;
 209     i__1 = olnint;
 210     for (p = 1; p <= i__1; ++p) {
 211         k = i__ << 1;
 212         left = work[k - 1];
 213         right = work[k];
 214         next = iwork[k - 1];
 215         nright = iwork[k];
 216         mid = (left + right) * .5;
 217         width = right - mid;
 218         d__1 = fabs(left);
 219         d__2 = fabs(right);
 220         tmp = (d__1>d__2) ? d__1 : d__2;
 221
 222         gap = 0.;
 223         if (i__ == nright) {
 224             if (prev > 0 && next <= *n) {
 225                 d__1 = left - work[k - 2], d__2 = work[k + 1] - right;
 226                 gap = (d__1<d__2) ? d__1 : d__2;
 227             } else if (prev > 0) {
 228                 gap = left - work[k - 2];
 229             } else if (next <= *n) {
 230                 gap = work[k + 1] - right;
 231             }
 232         }
 233         d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
 234         if (width < ((d__1>d__2) ? d__1 : d__2)) {
 235             --nint;
 236             iwork[k - 1] = 0;
 237             kk = k;
 238             i__2 = nright;
 239             for (j = i__ + 1; j <= i__2; ++j) {
 240                 kk += 2;
 241                 iwork[kk - 1] = 0;
 242                 work[kk - 1] = left;
 243                 work[kk] = right;
 244                 wgap[j - 1 - *offset] = 0.;
 245             }
 246             if (i1 == i__) {
 247                 i1 = next;
 248             } else {
 249                 iwork[(prev << 1) - 1] = next;
 250             }
 251             i__ = next;
 252             continue;
 253         }
 254         prev = i__;
 255
 256         s = -mid;
 257         cnt = 0;
 258         i__2 = *n - 1;
 259         for (j = 1; j <= i__2; ++j) {
 260             dplus = d__[j] + s;
 261             s = s * lld[j] / dplus - mid;
 262             if (dplus < 0.) {
 263                 ++cnt;
 264             }
 265         }
 266         dplus = d__[*n] + s;
 267         if (dplus < 0.) {
 268             ++cnt;
 269         }
 270         if (! (s > 0. || s < 1.)) {
 271             cnt = 0;
 272             s = -mid;
 273             i__2 = *n - 1;
 274             for (j = 1; j <= i__2; ++j) {
 275                 dplus = d__[j] + s;
 276                 if (dplus < 0.) {
 277                     ++cnt;
 278                 }
 279                 tmp = lld[j] / dplus;
 280                 if (fabs(tmp)<GMX_FLOAT_MIN) {
 281                     s = lld[j] - mid;
 282                 } else {
 283                     s = s * tmp - mid;
 284                 }
 285             }
 286             dplus = d__[*n] + s;
 287             if (dplus < 0.) {
 288                 ++cnt;
 289             }
 290         }
 291         i__2 = i__ - 1, i__3 = (nright<cnt) ? nright : cnt;
 292         cnt = (i__2>i__3) ? i__2 : i__3;
 293         if (cnt == i__ - 1) {
 294             work[k - 1] = mid;
 295         } else if (cnt == nright) {
 296             work[k] = mid;
 297         } else {
 298             iwork[k] = cnt;
 299             ++cnt;
 300             iwork[k - 1] = cnt;
 301             kk = cnt << 1;
 302             iwork[kk - 1] = next;
 303             iwork[kk] = nright;
 304             work[k] = mid;
 305             work[kk - 1] = mid;
 306             work[kk] = right;
 307             prev = cnt;
 308             if (cnt - 1 > i__) {
 309                 work[kk - 2] = mid;
 310             }
 311             if (cnt > *ifirst && cnt <= *ilast) {
 312                 ++nint;
 313             } else if (cnt <= *ifirst) {
 314                 i1 = cnt;
 315             }
 316         }
 317         i__ = next;
 318     }
 319     if (nint > 0) {
 320         goto L80;
 321     }
 322     i__1 = *ilast;
 323     for (i__ = *ifirst; i__ <= i__1; ++i__) {
 324         k = i__ << 1;
 325         ii = i__ - *offset;
 326         if (iwork[k - 1] != -1) {
 327             w[ii] = (work[k - 1] + work[k]) * .5;
 328             werr[ii] = work[k] - w[ii];
 329             if (i__ != *ilast) {
 330                 wgap[ii] = work[k + 1] - work[k];
 331             }
 332         }
 333     }
 334
 335     return;
 336
 337 }