src/gromacs/linearalgebra/gmx_lapack/slarrex.c

   1 #include <math.h>
   2 #include <ctype.h>
   3
   4 #include "gromacs/utility/real.h"
   5
   6 #include "../gmx_blas.h"
   7 #include "../gmx_lapack.h"
   8 #include "lapack_limits.h"
   9
  10
  11
  12 void
  13 F77_FUNC(slarrex,SLARREX)(const char *range,
  14          int *n,
  15          float *vl,
  16          float *vu,
  17          int *il,
  18          int *iu,
  19          float *d__,
  20          float *e,
  21          float *tol,
  22          int *nsplit,
  23          int *isplit,
  24          int *m,
  25          float *w,
  26          int *iblock,
  27          int *indexw,
  28          float *gersch,
  29          float *work,
  30          int *iwork,
  31          int *info)
  32 {
  33     int i__1, i__2, i__3;
  34     float d__1, d__2;
  35     int c__1 = 1;
  36     int c__0 = 0;
  37
  38     int i__, j, k;
  39     float s, gl;
  40     int in;
  41     float gu;
  42     int cnt;
  43     float eps, tau, nrm, tmp, vvl, vvu, offd;
  44     int iend, jblk, till, itmp;
  45     float rtol, delta, sigma;
  46     int iinfo;
  47     float width;
  48     int ibegin;
  49     int irange;
  50     float sgndef;
  51     int maxcnt;
  52     --iwork;
  53     --work;
  54     --gersch;
  55     --indexw;
  56     --iblock;
  57     --w;
  58     --isplit;
  59     --e;
  60     --d__;
  61
  62     sigma = 0;
  63     irange = 0;
  64     sgndef = 0;
  65     maxcnt = 0;
  66
  67     *info = 0;
  68
  69     if (*range=='A' || *range=='a')
  70         irange = 1;
  71     else if (*range=='V' || *range=='v')
  72         irange = 2;
  73     else if (*range=='I' || *range=='i')
  74         irange = 3;
  75
  76
  77     *m = 0;
  78     eps = GMX_FLOAT_EPS;
  79
  80     *nsplit = 1;
  81     i__1 = *n - 1;
  82     for (i__ = 1; i__ <= i__1; ++i__) {
  83         if (fabs(e[i__]) <= *tol) {
  84             isplit[*nsplit] = i__;
  85             ++(*nsplit);
  86         }
  87     }
  88     isplit[*nsplit] = *n;
  89
  90     ibegin = 1;
  91     i__1 = *nsplit;
  92     for (jblk = 1; jblk <= i__1; ++jblk) {
  93         iend = isplit[jblk];
  94         if (ibegin == iend) {
  95             ++(*m);
  96             w[*m] = d__[ibegin];
  97             iblock[*m] = jblk;
  98             indexw[*m] = 1;
  99             e[iend] = 0.;
 100             ibegin = iend + 1;
 101             goto L170;
 102         }
 103         in = iend - ibegin + 1;
 104
 105         gl = d__[ibegin] - fabs(e[ibegin]);
 106         gu = d__[ibegin] + fabs(e[ibegin]);
 107         gersch[(ibegin << 1) - 1] = gl;
 108         gersch[ibegin * 2] = gu;
 109         gersch[(iend << 1) - 1] = d__[iend] - fabs(e[iend - 1]);
 110         gersch[iend * 2] = d__[iend] + fabs(e[iend - 1]);
 111         d__1 = gersch[(iend << 1) - 1];
 112         gl = (d__1<gl) ? d__1 : gl;
 113         d__1 = gersch[iend * 2];
 114         gu = (d__1>gu) ? d__1 : gu;
 115         i__2 = iend - 1;
 116         for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
 117             offd = fabs(e[i__ - 1]) + fabs(e[i__]);
 118             gersch[(i__ << 1) - 1] = d__[i__] - offd;
 119             d__1 = gersch[(i__ << 1) - 1];
 120             gl = (d__1<gl) ? d__1 : gl;
 121             gersch[i__ * 2] = d__[i__] + offd;
 122             d__1 = gersch[i__ * 2];
 123             gu = (d__1>gu) ? d__1 : gu;
 124         }
 125         d__1 = fabs(gl), d__2 = fabs(gu);
 126         nrm = (d__1>d__2) ? d__1 : d__2;
 127
 128         width = gu - gl;
 129         i__2 = iend - 1;
 130         for (i__ = ibegin; i__ <= i__2; ++i__) {
 131             work[i__] = e[i__] * e[i__];
 132         }
 133         for (j = 1; j <= 2; ++j) {
 134             if (j == 1) {
 135                 tau = gl + width * .25;
 136             } else {
 137                 tau = gu - width * .25;
 138             }
 139             tmp = d__[ibegin] - tau;
 140             if (tmp < 0.) {
 141                 cnt = 1;
 142             } else {
 143                 cnt = 0;
 144             }
 145             i__2 = iend;
 146             for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
 147                 tmp = d__[i__] - tau - work[i__ - 1] / tmp;
 148                 if (tmp < 0.) {
 149                     ++cnt;
 150                 }
 151             }
 152             if (cnt == 0) {
 153                 gl = tau;
 154             } else if (cnt == in) {
 155                 gu = tau;
 156             }
 157             if (j == 1) {
 158                 maxcnt = cnt;
 159                 sigma = gl;
 160                 sgndef = 1.;
 161             } else {
 162                 if (in - cnt > maxcnt) {
 163                     sigma = gu;
 164                     sgndef = -1.;
 165                 }
 166             }
 167         }
 168
 169         work[in * 3] = 1.;
 170         delta = eps;
 171         tau = sgndef * nrm;
 172 L60:
 173         sigma -= delta * tau;
 174         work[1] = d__[ibegin] - sigma;
 175         j = ibegin;
 176         i__2 = in - 1;
 177         for (i__ = 1; i__ <= i__2; ++i__) {
 178             work[(in << 1) + i__] = 1. / work[i__];
 179             tmp = e[j] * work[(in << 1) + i__];
 180             work[i__ + 1] = d__[j + 1] - sigma - tmp * e[j];
 181             work[in + i__] = tmp;
 182             ++j;
 183         }
 184         for (i__ = in; i__ >= 1; --i__) {
 185             tmp = sgndef * work[i__];
 186             if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_FLOAT_MIN || ! (tmp > 0. || tmp < 1.)) {
 187                 delta *= 2.;
 188                 goto L60;
 189             }
 190         }
 191
 192         F77_FUNC(scopy,SCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
 193         i__2 = in - 1;
 194         F77_FUNC(scopy,SCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
 195         i__2 = in - 1;
 196         for (i__ = 1; i__ <= i__2; ++i__) {
 197             work[in * 3 + i__] = work[i__] * work[in + i__];
 198             work[(in << 2) + i__] = work[in * 3 + i__] * work[in + i__];
 199         }
 200         if (sgndef > 0.) {
 201             cnt = 1;
 202             work[1] = (gl + gu) / 2. - sigma;
 203             work[in + 1] = 0.;
 204             work[(in << 1) + 1] = (gu - gl) / 2.;
 205         } else {
 206             cnt = in;
 207             work[in] = (gl + gu) / 2. - sigma;
 208             work[in * 2] = 0.;
 209             work[in * 3] = (gu - gl) / 2.;
 210         }
 211         rtol = eps * 4.;
 212         F77_FUNC(slarrbx,SLARRBX)(&in, &d__[ibegin], &e[ibegin], &work[in * 3 + 1], &work[(in <<
 213                  2) + 1], &cnt, &cnt, &rtol, &rtol, &c__0, &work[1], &work[in
 214                 + 1], &work[(in << 1) + 1], &work[in * 5 + 1], &iwork[1], &
 215                 iinfo);
 216         if (sgndef > 0.) {
 217             tau = work[1] - work[(in << 1) + 1];
 218         } else {
 219             tau = work[in] + work[in * 3];
 220         }
 221
 222         work[in * 3] = 1.;
 223         delta = eps * 2.;
 224 L100:
 225         tau *= 1. - delta;
 226
 227         s = -tau;
 228         j = ibegin;
 229         i__2 = in - 1;
 230         for (i__ = 1; i__ <= i__2; ++i__) {
 231             work[i__] = d__[j] + s;
 232             work[(in << 1) + i__] = 1. / work[i__];
 233             work[in + i__] = e[j] * d__[j] * work[(in << 1) + i__];
 234             s = s * work[in + i__] * e[j] - tau;
 235             ++j;
 236         }
 237         work[in] = d__[iend] + s;
 238
 239         for (i__ = in; i__ >= 1; --i__) {
 240             tmp = sgndef * work[i__];
 241             if (tmp < 0. || fabs(work[(in << 1) + i__])<GMX_FLOAT_MIN || ! (tmp > 0. || tmp < 1.)) {
 242                 delta *= 2.;
 243                 goto L100;
 244             }
 245         }
 246
 247         sigma += tau;
 248         F77_FUNC(scopy,SCOPY)(&in, &work[1], &c__1, &d__[ibegin], &c__1);
 249         i__2 = in - 1;
 250         F77_FUNC(scopy,SCOPY)(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
 251         e[iend] = sigma;
 252         tmp = (float) in * 4. * eps * (fabs(sigma) + fabs(tau));
 253         i__2 = iend;
 254         for (i__ = ibegin; i__ <= i__2; ++i__) {
 255             gersch[(i__ << 1) - 1] = gersch[(i__ << 1) - 1] - sigma - tmp;
 256             gersch[i__ * 2] = gersch[i__ * 2] - sigma + tmp;
 257         }
 258
 259         j = ibegin;
 260         i__2 = in - 1;
 261         for (i__ = 1; i__ <= i__2; ++i__) {
 262             work[(i__ << 1) - 1] = fabs(d__[j]);
 263             work[i__ * 2] = e[j] * e[j] * work[(i__ << 1) - 1];
 264             ++j;
 265         }
 266         work[(in << 1) - 1] = fabs(d__[iend]);
 267
 268         F77_FUNC(slasq2,SLASQ2)(&in, &work[1], info);
 269         if (*info != 0) {
 270             return;
 271         }
 272
 273         if (sgndef > 0.) {
 274             i__2 = in;
 275             for (i__ = 1; i__ <= i__2; ++i__) {
 276                 ++(*m);
 277                 w[*m] = work[in - i__ + 1];
 278                 iblock[*m] = jblk;
 279                 indexw[*m] = i__;
 280             }
 281         } else {
 282             i__2 = in;
 283             for (i__ = 1; i__ <= i__2; ++i__) {
 284                 ++(*m);
 285                 w[*m] = -work[i__];
 286                 iblock[*m] = jblk;
 287                 indexw[*m] = i__;
 288             }
 289         }
 290         ibegin = iend + 1;
 291 L170:
 292         ;
 293     }
 294     if (irange == 2) {
 295         *m = 0;
 296         ibegin = 1;
 297         i__1 = *nsplit;
 298         for (i__ = 1; i__ <= i__1; ++i__) {
 299             iend = isplit[i__];
 300             vvl = *vl - e[iend];
 301             vvu = *vu - e[iend];
 302             i__2 = iend;
 303             for (j = ibegin; j <= i__2; ++j) {
 304                 if (vvl <= w[j] && w[j] <= vvu) {
 305                     ++(*m);
 306                     w[*m] = w[j];
 307                     iblock[*m] = i__;
 308                     indexw[*m] = j - ibegin + 1;
 309                 }
 310             }
 311             ibegin = iend + 1;
 312         }
 313     } else if (irange == 3) {
 314         *m = *iu - *il + 1;
 315         if (*nsplit == 1) {
 316             i__1 = *m;
 317             for (i__ = 1; i__ <= i__1; ++i__) {
 318                 w[i__] = w[*il + i__ - 1];
 319                 indexw[i__] = *il + i__ - 1;
 320             }
 321         } else {
 322             ibegin = 1;
 323             i__1 = *nsplit;
 324             for (i__ = 1; i__ <= i__1; ++i__) {
 325                 iend = isplit[i__];
 326                 i__2 = iend;
 327                 for (j = ibegin; j <= i__2; ++j) {
 328                     work[j] = w[j] + e[iend];
 329                 }
 330                 ibegin = iend + 1;
 331             }
 332             i__1 = *n;
 333             for (i__ = 1; i__ <= i__1; ++i__) {
 334                 iwork[i__] = i__;
 335                 iwork[*n + i__] = iblock[i__];
 336             }
 337             F77_FUNC(slasrt2,SLASRT2)("I", n, &work[1], &iwork[1], &iinfo);
 338             i__1 = *m;
 339             for (i__ = 1; i__ <= i__1; ++i__) {
 340                 itmp = iwork[*il + i__ - 1];
 341                 work[i__] = w[itmp];
 342                 iblock[i__] = iwork[*n + itmp];
 343             }
 344             i__1 = *m;
 345             for (i__ = 1; i__ <= i__1; ++i__) {
 346                 iwork[*n + i__] = iwork[*il + i__ - 1];
 347                 iwork[i__] = i__;
 348             }
 349             F77_FUNC(ilasrt2,ILASRT2)("I", m, &iblock[1], &iwork[1], &iinfo);
 350             j = 1;
 351             itmp = iblock[j];
 352             cnt = iwork[*n + iwork[j]];
 353             if (itmp == 1) {
 354                 ibegin = 1;
 355             } else {
 356                 ibegin = isplit[itmp - 1] + 1;
 357             }
 358             i__1 = *m;
 359             for (i__ = 1; i__ <= i__1; ++i__) {
 360                 w[i__] = work[iwork[i__]];
 361                 if (iblock[i__] != itmp || i__ == *m) {
 362                     if (iblock[i__] == itmp) {
 363                         till = *m;
 364                     } else {
 365                         till = i__ - 1;
 366                     }
 367                     i__2 = till - j + 1;
 368                     F77_FUNC(slasrt,SLASRT)("I", &i__2, &w[j], &iinfo);
 369                     cnt = cnt - ibegin + 1;
 370                     i__2 = till;
 371                     for (k = j; k <= i__2; ++k) {
 372                         indexw[k] = cnt + k - j;
 373                     }
 374                     j = i__;
 375                     itmp = iblock[j];
 376                     cnt = iwork[*n + iwork[j]];
 377                     ibegin = isplit[itmp - 1] + 1;
 378                     if (i__ == *m && till < *m) {
 379                         indexw[*m] = cnt - ibegin + 1;
 380                     }
 381                 } else {
 382                     i__2 = cnt, i__3 = iwork[*n + iwork[i__]];
 383                     cnt = (i__2<i__3) ? i__2 : i__3;
 384                 }
 385             }
 386         }
 387     }
 388
 389     return;
 390
 391 }
 392
 393