src/gromacs/linearalgebra/gmx_lapack/sstebz.c

   1 #include <math.h>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "gromacs/utility/real.h"
   6
   7 void
   8 F77_FUNC(sstebz,SSTEBZ)(const char *range,
   9                         const char *order,
  10                         int *n,
  11                         float *vl,
  12                         float *vu,
  13                         int *il,
  14                         int *iu,
  15                         float *abstol,
  16                         float *d__,
  17                         float *e,
  18                         int *m,
  19                         int *nsplit,
  20                         float *w,
  21                         int *iblock,
  22                         int *isplit,
  23                         float *work,
  24                         int *iwork,
  25                         int *info)
  26 {
  27     int i__1, i__2, i__3;
  28     float d__1, d__2, d__3, d__4, d__5;
  29     int c__1 = 1;
  30     int c__3 = 3;
  31     int c__2 = 2;
  32     int c__0 = 0;
  33
  34     int j, ib, jb, ie, je, nb;
  35     float gl;
  36     int im, in;
  37     float gu;
  38     int iw;
  39     float wl, wu;
  40     int nwl;
  41     float ulp, wlu, wul;
  42     int nwu;
  43     float tmp1, tmp2;
  44     int iend, ioff, iout, itmp1, jdisc;
  45     int iinfo;
  46     float atoli;
  47     int iwoff;
  48     float bnorm;
  49     int itmax;
  50     float wkill, rtoli, tnorm;
  51     int ibegin;
  52     int irange, idiscl;
  53     int idumma[1];
  54     int idiscu, iorder;
  55     int ncnvrg;
  56     float pivmin;
  57     int toofew;
  58     const float safemn = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
  59
  60     --iwork;
  61     --work;
  62     --isplit;
  63     --iblock;
  64     --w;
  65     --e;
  66     --d__;
  67
  68     *info = 0;
  69
  70     if (*range=='A' || *range=='a') {
  71         irange = 1;
  72     } else if (*range=='V' || *range=='v') {
  73         irange = 2;
  74     } else if (*range=='I' || *range=='i') {
  75         irange = 3;
  76     } else {
  77         irange = 0;
  78     }
  79
  80     if (*order=='B' || *order=='b') {
  81         iorder = 2;
  82     } else if (*order=='E' || *order=='e') {
  83         iorder = 1;
  84     } else {
  85         iorder = 0;
  86     }
  87
  88     if (irange <= 0) {
  89         *info = -1;
  90     } else if (iorder <= 0) {
  91         *info = -2;
  92     } else if (*n < 0) {
  93         *info = -3;
  94     } else if (irange == 2) {
  95         if (*vl >= *vu) {
  96             *info = -5;
  97         }
  98     } else if (irange == 3 && (*il < 1 || *il > (*n))) {
  99         *info = -6;
 100     } else if (irange == 3 && (*iu < ((*n<*il) ? *n : *il) || *iu > *n)) {
 101         *info = -7;
 102     }
 103
 104     if (*info != 0) {
 105         i__1 = -(*info);
 106         return;
 107     }
 108
 109     *info = 0;
 110     ncnvrg = 0;
 111     toofew = 0;
 112
 113     *m = 0;
 114     if (*n == 0) {
 115         return;
 116     }
 117
 118     if (irange == 3 && *il == 1 && *iu == *n) {
 119         irange = 1;
 120     }
 121
 122     ulp = 2*GMX_FLOAT_EPS;
 123     rtoli = ulp * 2.;
 124     nb = DSTEBZ_BLOCKSIZE;
 125     if (nb <= 1) {
 126         nb = 0;
 127     }
 128
 129     if (*n == 1) {
 130         *nsplit = 1;
 131         isplit[1] = 1;
 132         if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
 133             *m = 0;
 134         } else {
 135             w[1] = d__[1];
 136             iblock[1] = 1;
 137             *m = 1;
 138         }
 139         return;
 140     }
 141
 142     *nsplit = 1;
 143     work[*n] = 0.;
 144     pivmin = 1.;
 145     i__1 = *n;
 146     for (j = 2; j <= i__1; ++j) {
 147         d__1 = e[j - 1];
 148         tmp1 = d__1 * d__1;
 149         d__2 = ulp;
 150         if (fabs(d__[j] * d__[j - 1]) * (d__2 * d__2) + safemn
 151                 > tmp1) {
 152             isplit[*nsplit] = j - 1;
 153             ++(*nsplit);
 154             work[j - 1] = 0.;
 155         } else {
 156             work[j - 1] = tmp1;
 157             pivmin = (pivmin>tmp1) ? pivmin : tmp1;
 158         }
 159     }
 160     isplit[*nsplit] = *n;
 161     pivmin *= safemn;
 162
 163     if (irange == 3) {
 164
 165         gu = d__[1];
 166         gl = d__[1];
 167         tmp1 = 0.;
 168
 169         i__1 = *n - 1;
 170         for (j = 1; j <= i__1; ++j) {
 171             tmp2 = sqrt(work[j]);
 172             d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 173             gu = (d__1>d__2) ? d__1 : d__2;
 174             d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 175             gl = (d__1<d__2) ? d__1 : d__2;
 176             tmp1 = tmp2;
 177         }
 178
 179         d__1 = gu, d__2 = d__[*n] + tmp1;
 180         gu = (d__1>d__2) ? d__1 : d__2;
 181         d__1 = gl, d__2 = d__[*n] - tmp1;
 182         gl = (d__1<d__2) ? d__1 : d__2;
 183         d__1 = fabs(gl);
 184         d__2 = fabs(gu);
 185         tnorm = (d__1>d__2) ? d__1 : d__2;
 186         gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.;
 187         gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.;
 188
 189         itmax = (int) ((log(tnorm + pivmin) - log(pivmin)) / log(2.)) + 2;
 190         if (*abstol <= 0.) {
 191             atoli = ulp * tnorm;
 192         } else {
 193             atoli = *abstol;
 194         }
 195
 196         work[*n + 1] = gl;
 197         work[*n + 2] = gl;
 198         work[*n + 3] = gu;
 199         work[*n + 4] = gu;
 200         work[*n + 5] = gl;
 201         work[*n + 6] = gu;
 202         iwork[1] = -1;
 203         iwork[2] = -1;
 204         iwork[3] = *n + 1;
 205         iwork[4] = *n + 1;
 206         iwork[5] = *il - 1;
 207         iwork[6] = *iu;
 208
 209         F77_FUNC(slaebz,SLAEBZ)(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
 210                 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
 211                 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
 212
 213         if (iwork[6] == *iu) {
 214             wl = work[*n + 1];
 215             wlu = work[*n + 3];
 216             nwl = iwork[1];
 217             wu = work[*n + 4];
 218             wul = work[*n + 2];
 219             nwu = iwork[4];
 220         } else {
 221             wl = work[*n + 2];
 222             wlu = work[*n + 4];
 223             nwl = iwork[2];
 224             wu = work[*n + 3];
 225             wul = work[*n + 1];
 226             nwu = iwork[3];
 227         }
 228
 229         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
 230             *info = 4;
 231             return;
 232         }
 233     } else {
 234
 235
 236       /* avoid warnings for high gcc optimization */
 237       wlu = wul = 1.0;
 238
 239         d__3 = fabs(d__[1]) + fabs(e[1]);
 240         d__4 = fabs(d__[*n]) + fabs(e[*n - 1]);
 241         tnorm = (d__3>d__4) ? d__3 : d__4;
 242
 243         i__1 = *n - 1;
 244         for (j = 2; j <= i__1; ++j) {
 245             d__4 = tnorm;
 246             d__5 = fabs(d__[j]) + fabs(e[j - 1]) + fabs(e[j]);
 247             tnorm = (d__4>d__5) ? d__4 : d__5;
 248         }
 249
 250         if (*abstol <= 0.) {
 251             atoli = ulp * tnorm;
 252         } else {
 253             atoli = *abstol;
 254         }
 255
 256         if (irange == 2) {
 257             wl = *vl;
 258             wu = *vu;
 259         } else {
 260             wl = 0.;
 261             wu = 0.;
 262         }
 263     }
 264
 265     *m = 0;
 266     iend = 0;
 267     *info = 0;
 268     nwl = 0;
 269     nwu = 0;
 270
 271     i__1 = *nsplit;
 272     for (jb = 1; jb <= i__1; ++jb) {
 273         ioff = iend;
 274         ibegin = ioff + 1;
 275         iend = isplit[jb];
 276         in = iend - ioff;
 277
 278         if (in == 1) {
 279
 280             if (irange == 1 || wl >= d__[ibegin] - pivmin) {
 281                 ++nwl;
 282             }
 283             if (irange == 1 || wu >= d__[ibegin] - pivmin) {
 284                 ++nwu;
 285             }
 286             if (irange == 1 || ((wl < d__[ibegin] - pivmin) && (wu >= d__[ibegin] - pivmin))) {
 287                 ++(*m);
 288                 w[*m] = d__[ibegin];
 289                 iblock[*m] = jb;
 290             }
 291         } else {
 292
 293             gu = d__[ibegin];
 294             gl = d__[ibegin];
 295             tmp1 = 0.;
 296
 297             i__2 = iend - 1;
 298             for (j = ibegin; j <= i__2; ++j) {
 299                 tmp2 = fabs(e[j]);
 300                 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 301                 gu = (d__1>d__2) ? d__1 : d__2;
 302                 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 303                 gl = (d__1<d__2) ? d__1 : d__2;
 304                 tmp1 = tmp2;
 305             }
 306
 307             d__1 = gu, d__2 = d__[iend] + tmp1;
 308             gu = (d__1>d__2) ? d__1 : d__2;
 309             d__1 = gl, d__2 = d__[iend] - tmp1;
 310             gl = (d__1<d__2) ? d__1 : d__2;
 311             d__1 = fabs(gl);
 312             d__2 = fabs(gu);
 313             bnorm = (d__1>d__2) ? d__1 : d__2;
 314             gl = gl - bnorm * 2. * ulp * in - pivmin * 2.;
 315             gu = gu + bnorm * 2. * ulp * in + pivmin * 2.;
 316
 317             if (*abstol <= 0.) {
 318                 d__1 = fabs(gl);
 319                 d__2 = fabs(gu);
 320                 atoli = ulp * ((d__1>d__2) ? d__1 : d__2);
 321             } else {
 322                 atoli = *abstol;
 323             }
 324
 325             if (irange > 1) {
 326                 if (gu < wl) {
 327                     nwl += in;
 328                     nwu += in;
 329                 }
 330                 gl = (gl>wl) ? gl : wl;
 331                 gu = (gu<wu) ? gu : wu;
 332                 if (gl >= gu) {
 333                 }
 334                 continue;
 335             }
 336
 337             work[*n + 1] = gl;
 338             work[*n + in + 1] = gu;
 339             F77_FUNC(slaebz,SLAEBZ)(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 340                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 341                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
 342                     w[*m + 1], &iblock[*m + 1], &iinfo);
 343
 344             nwl += iwork[1];
 345             nwu += iwork[in + 1];
 346             iwoff = *m - iwork[1];
 347
 348             itmax = (int) ((log(gu - gl + pivmin) - log(pivmin)) / log(2.)
 349                     ) + 2;
 350             F77_FUNC(slaebz,SLAEBZ)(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 351                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 352                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
 353                      &w[*m + 1], &iblock[*m + 1], &iinfo);
 354
 355             i__2 = iout;
 356             for (j = 1; j <= i__2; ++j) {
 357                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
 358
 359                 if (j > iout - iinfo) {
 360                     ncnvrg = 1;
 361                     ib = -jb;
 362                 } else {
 363                     ib = jb;
 364                 }
 365                 i__3 = iwork[j + in] + iwoff;
 366                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
 367                     w[je] = tmp1;
 368                     iblock[je] = ib;
 369                 }
 370             }
 371
 372             *m += im;
 373         }
 374     }
 375
 376     if (irange == 3) {
 377         im = 0;
 378         idiscl = *il - 1 - nwl;
 379         idiscu = nwu - *iu;
 380
 381         if (idiscl > 0 || idiscu > 0) {
 382             i__1 = *m;
 383             for (je = 1; je <= i__1; ++je) {
 384                 if (w[je] <= wlu && idiscl > 0) {
 385                     --idiscl;
 386                 } else if (w[je] >= wul && idiscu > 0) {
 387                     --idiscu;
 388                 } else {
 389                     ++im;
 390                     w[im] = w[je];
 391                     iblock[im] = iblock[je];
 392                 }
 393             }
 394             *m = im;
 395         }
 396         if (idiscl > 0 || idiscu > 0) {
 397
 398             if (idiscl > 0) {
 399                 wkill = wu;
 400                 i__1 = idiscl;
 401                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 402                     iw = 0;
 403                     i__2 = *m;
 404                     for (je = 1; je <= i__2; ++je) {
 405                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
 406                             iw = je;
 407                             wkill = w[je];
 408                         }
 409                     }
 410                     iblock[iw] = 0;
 411                 }
 412             }
 413             if (idiscu > 0) {
 414
 415                 wkill = wl;
 416                 i__1 = idiscu;
 417                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 418                     iw = 0;
 419                     i__2 = *m;
 420                     for (je = 1; je <= i__2; ++je) {
 421                         if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
 422                             iw = je;
 423                             wkill = w[je];
 424                         }
 425                     }
 426                     iblock[iw] = 0;
 427                 }
 428             }
 429             im = 0;
 430             i__1 = *m;
 431             for (je = 1; je <= i__1; ++je) {
 432                 if (iblock[je] != 0) {
 433                     ++im;
 434                     w[im] = w[je];
 435                     iblock[im] = iblock[je];
 436                 }
 437             }
 438             *m = im;
 439         }
 440         if (idiscl < 0 || idiscu < 0) {
 441             toofew = 1;
 442         }
 443     }
 444
 445     if (iorder == 1 && *nsplit > 1) {
 446         i__1 = *m - 1;
 447         for (je = 1; je <= i__1; ++je) {
 448             ie = 0;
 449             tmp1 = w[je];
 450             i__2 = *m;
 451             for (j = je + 1; j <= i__2; ++j) {
 452                 if (w[j] < tmp1) {
 453                     ie = j;
 454                     tmp1 = w[j];
 455                 }
 456             }
 457
 458             if (ie != 0) {
 459                 itmp1 = iblock[ie];
 460                 w[ie] = w[je];
 461                 iblock[ie] = iblock[je];
 462                 w[je] = tmp1;
 463                 iblock[je] = itmp1;
 464             }
 465         }
 466     }
 467
 468     *info = 0;
 469     if (ncnvrg) {
 470         ++(*info);
 471     }
 472     if (toofew) {
 473         *info += 2;
 474     }
 475     return;
 476
 477 }
 478
 479