src/gmxlib/gmx_lapack/dstebz.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include <math.h>
  36 #include "gmx_lapack.h"
  37 #include "lapack_limits.h"
  38
  39 #include <types/simple.h>
  40
  41 void
  42 F77_FUNC(dstebz,DSTEBZ)(const char *range,
  43                         const char *order,
  44                         int *n,
  45                         double *vl,
  46                         double *vu,
  47                         int *il,
  48                         int *iu,
  49                         double *abstol,
  50                         double *d__,
  51                         double *e,
  52                         int *m,
  53                         int *nsplit,
  54                         double *w,
  55                         int *iblock,
  56                         int *isplit,
  57                         double *work,
  58                         int *iwork,
  59                         int *info)
  60 {
  61     int i__1, i__2, i__3;
  62     double d__1, d__2, d__3, d__4, d__5;
  63     int c__1 = 1;
  64     int c__3 = 3;
  65     int c__2 = 2;
  66     int c__0 = 0;
  67
  68     int j, ib, jb, ie, je, nb;
  69     double gl;
  70     int im, in;
  71     double gu;
  72     int iw;
  73     double wl, wu;
  74     int nwl;
  75     double ulp, wlu, wul;
  76     int nwu;
  77     double tmp1, tmp2;
  78     int iend, ioff, iout, itmp1, jdisc;
  79     int iinfo;
  80     double atoli;
  81     int iwoff;
  82     double bnorm;
  83     int itmax;
  84     double wkill, rtoli, tnorm;
  85     int ibegin;
  86     int irange, idiscl;
  87     int idumma[1];
  88     int idiscu, iorder;
  89     int ncnvrg;
  90     double pivmin;
  91     int toofew;
  92     const double safemn = GMX_DOUBLE_MIN*(1.0+GMX_DOUBLE_EPS);
  93
  94     --iwork;
  95     --work;
  96     --isplit;
  97     --iblock;
  98     --w;
  99     --e;
 100     --d__;
 101
 102     *info = 0;
 103
 104     if (*range=='A' || *range=='a') {
 105         irange = 1;
 106     } else if (*range=='V' || *range=='v') {
 107         irange = 2;
 108     } else if (*range=='I' || *range=='i') {
 109         irange = 3;
 110     } else {
 111         irange = 0;
 112     }
 113
 114     if (*order=='B' || *order=='b') {
 115         iorder = 2;
 116     } else if (*order=='E' || *order=='e') {
 117         iorder = 1;
 118     } else {
 119         iorder = 0;
 120     }
 121
 122     if (irange <= 0) {
 123         *info = -1;
 124     } else if (iorder <= 0) {
 125         *info = -2;
 126     } else if (*n < 0) {
 127         *info = -3;
 128     } else if (irange == 2) {
 129         if (*vl >= *vu) {
 130             *info = -5;
 131         }
 132     } else if (irange == 3 && (*il < 1 || *il > (*n))) {
 133         *info = -6;
 134     } else if (irange == 3 && (*iu < ((*n<*il) ? *n : *il) || *iu > *n)) {
 135         *info = -7;
 136     }
 137
 138     if (*info != 0) {
 139         i__1 = -(*info);
 140         return;
 141     }
 142
 143     *info = 0;
 144     ncnvrg = 0;
 145     toofew = 0;
 146
 147     *m = 0;
 148     if (*n == 0) {
 149         return;
 150     }
 151
 152     if (irange == 3 && *il == 1 && *iu == *n) {
 153         irange = 1;
 154     }
 155
 156     ulp = 2*GMX_DOUBLE_EPS;
 157     rtoli = ulp * 2.;
 158     nb = DSTEBZ_BLOCKSIZE;
 159     if (nb <= 1) {
 160         nb = 0;
 161     }
 162
 163     if (*n == 1) {
 164         *nsplit = 1;
 165         isplit[1] = 1;
 166         if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
 167             *m = 0;
 168         } else {
 169             w[1] = d__[1];
 170             iblock[1] = 1;
 171             *m = 1;
 172         }
 173         return;
 174     }
 175
 176     *nsplit = 1;
 177     work[*n] = 0.;
 178     pivmin = 1.;
 179     i__1 = *n;
 180     for (j = 2; j <= i__1; ++j) {
 181         d__1 = e[j - 1];
 182         tmp1 = d__1 * d__1;
 183         d__2 = ulp;
 184         if (fabs(d__[j] * d__[j - 1]) * (d__2 * d__2) + safemn
 185                 > tmp1) {
 186             isplit[*nsplit] = j - 1;
 187             ++(*nsplit);
 188             work[j - 1] = 0.;
 189         } else {
 190             work[j - 1] = tmp1;
 191             pivmin = (pivmin>tmp1) ? pivmin : tmp1;
 192         }
 193     }
 194     isplit[*nsplit] = *n;
 195     pivmin *= safemn;
 196
 197     if (irange == 3) {
 198
 199         gu = d__[1];
 200         gl = d__[1];
 201         tmp1 = 0.;
 202
 203         i__1 = *n - 1;
 204         for (j = 1; j <= i__1; ++j) {
 205             tmp2 = sqrt(work[j]);
 206             d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 207             gu = (d__1>d__2) ? d__1 : d__2;
 208             d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 209             gl = (d__1<d__2) ? d__1 : d__2;
 210             tmp1 = tmp2;
 211         }
 212
 213         d__1 = gu, d__2 = d__[*n] + tmp1;
 214         gu = (d__1>d__2) ? d__1 : d__2;
 215         d__1 = gl, d__2 = d__[*n] - tmp1;
 216         gl = (d__1<d__2) ? d__1 : d__2;
 217         d__1 = fabs(gl);
 218         d__2 = fabs(gu);
 219         tnorm = (d__1>d__2) ? d__1 : d__2;
 220         gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.;
 221         gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.;
 222
 223         itmax = (int) ((log(tnorm + pivmin) - log(pivmin)) / log(2.)) + 2;
 224         if (*abstol <= 0.) {
 225             atoli = ulp * tnorm;
 226         } else {
 227             atoli = *abstol;
 228         }
 229
 230         work[*n + 1] = gl;
 231         work[*n + 2] = gl;
 232         work[*n + 3] = gu;
 233         work[*n + 4] = gu;
 234         work[*n + 5] = gl;
 235         work[*n + 6] = gu;
 236         iwork[1] = -1;
 237         iwork[2] = -1;
 238         iwork[3] = *n + 1;
 239         iwork[4] = *n + 1;
 240         iwork[5] = *il - 1;
 241         iwork[6] = *iu;
 242
 243         F77_FUNC(dlaebz,DLAEBZ)(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
 244                 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
 245                 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
 246
 247         if (iwork[6] == *iu) {
 248             wl = work[*n + 1];
 249             wlu = work[*n + 3];
 250             nwl = iwork[1];
 251             wu = work[*n + 4];
 252             wul = work[*n + 2];
 253             nwu = iwork[4];
 254         } else {
 255             wl = work[*n + 2];
 256             wlu = work[*n + 4];
 257             nwl = iwork[2];
 258             wu = work[*n + 3];
 259             wul = work[*n + 1];
 260             nwu = iwork[3];
 261         }
 262
 263         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
 264             *info = 4;
 265             return;
 266         }
 267     } else {
 268
 269
 270       /* avoid warnings for high gcc optimization */
 271       wlu = wul = 1.0;
 272
 273         d__3 = fabs(d__[1]) + fabs(e[1]);
 274         d__4 = fabs(d__[*n]) + fabs(e[*n - 1]);
 275         tnorm = (d__3>d__4) ? d__3 : d__4;
 276
 277         i__1 = *n - 1;
 278         for (j = 2; j <= i__1; ++j) {
 279             d__4 = tnorm;
 280             d__5 = fabs(d__[j]) + fabs(e[j - 1]) + fabs(e[j]);
 281             tnorm = (d__4>d__5) ? d__4 : d__5;
 282         }
 283
 284         if (*abstol <= 0.) {
 285             atoli = ulp * tnorm;
 286         } else {
 287             atoli = *abstol;
 288         }
 289
 290         if (irange == 2) {
 291             wl = *vl;
 292             wu = *vu;
 293         } else {
 294             wl = 0.;
 295             wu = 0.;
 296         }
 297     }
 298
 299     *m = 0;
 300     iend = 0;
 301     *info = 0;
 302     nwl = 0;
 303     nwu = 0;
 304
 305     i__1 = *nsplit;
 306     for (jb = 1; jb <= i__1; ++jb) {
 307         ioff = iend;
 308         ibegin = ioff + 1;
 309         iend = isplit[jb];
 310         in = iend - ioff;
 311
 312         if (in == 1) {
 313
 314             if (irange == 1 || wl >= d__[ibegin] - pivmin) {
 315                 ++nwl;
 316             }
 317             if (irange == 1 || wu >= d__[ibegin] - pivmin) {
 318                 ++nwu;
 319             }
 320             if (irange == 1 || ((wl < d__[ibegin] - pivmin) && (wu >= d__[ibegin] - pivmin))) {
 321                 ++(*m);
 322                 w[*m] = d__[ibegin];
 323                 iblock[*m] = jb;
 324             }
 325         } else {
 326
 327             gu = d__[ibegin];
 328             gl = d__[ibegin];
 329             tmp1 = 0.;
 330
 331             i__2 = iend - 1;
 332             for (j = ibegin; j <= i__2; ++j) {
 333                 tmp2 = fabs(e[j]);
 334                 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 335                 gu = (d__1>d__2) ? d__1 : d__2;
 336                 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 337                 gl = (d__1<d__2) ? d__1 : d__2;
 338                 tmp1 = tmp2;
 339             }
 340
 341             d__1 = gu, d__2 = d__[iend] + tmp1;
 342             gu = (d__1>d__2) ? d__1 : d__2;
 343             d__1 = gl, d__2 = d__[iend] - tmp1;
 344             gl = (d__1<d__2) ? d__1 : d__2;
 345             d__1 = fabs(gl);
 346             d__2 = fabs(gu);
 347             bnorm = (d__1>d__2) ? d__1 : d__2;
 348             gl = gl - bnorm * 2. * ulp * in - pivmin * 2.;
 349             gu = gu + bnorm * 2. * ulp * in + pivmin * 2.;
 350
 351             if (*abstol <= 0.) {
 352                 d__1 = fabs(gl);
 353                 d__2 = fabs(gu);
 354                 atoli = ulp * ((d__1>d__2) ? d__1 : d__2);
 355             } else {
 356                 atoli = *abstol;
 357             }
 358
 359             if (irange > 1) {
 360                 if (gu < wl) {
 361                     nwl += in;
 362                     nwu += in;
 363                 }
 364                 gl = (gl>wl) ? gl : wl;
 365                 gu = (gu<wu) ? gu : wu;
 366                 if (gl >= gu) {
 367                 }
 368                 continue;
 369             }
 370
 371             work[*n + 1] = gl;
 372             work[*n + in + 1] = gu;
 373             F77_FUNC(dlaebz,DLAEBZ)(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 374                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 375                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
 376                     w[*m + 1], &iblock[*m + 1], &iinfo);
 377
 378             nwl += iwork[1];
 379             nwu += iwork[in + 1];
 380             iwoff = *m - iwork[1];
 381
 382             itmax = (int) ((log(gu - gl + pivmin) - log(pivmin)) / log(2.)
 383                     ) + 2;
 384             F77_FUNC(dlaebz,DLAEBZ)(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 385                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 386                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
 387                      &w[*m + 1], &iblock[*m + 1], &iinfo);
 388
 389             i__2 = iout;
 390             for (j = 1; j <= i__2; ++j) {
 391                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
 392
 393                 if (j > iout - iinfo) {
 394                     ncnvrg = 1;
 395                     ib = -jb;
 396                 } else {
 397                     ib = jb;
 398                 }
 399                 i__3 = iwork[j + in] + iwoff;
 400                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
 401                     w[je] = tmp1;
 402                     iblock[je] = ib;
 403                 }
 404             }
 405
 406             *m += im;
 407         }
 408     }
 409
 410     if (irange == 3) {
 411         im = 0;
 412         idiscl = *il - 1 - nwl;
 413         idiscu = nwu - *iu;
 414
 415         if (idiscl > 0 || idiscu > 0) {
 416             i__1 = *m;
 417             for (je = 1; je <= i__1; ++je) {
 418                 if (w[je] <= wlu && idiscl > 0) {
 419                     --idiscl;
 420                 } else if (w[je] >= wul && idiscu > 0) {
 421                     --idiscu;
 422                 } else {
 423                     ++im;
 424                     w[im] = w[je];
 425                     iblock[im] = iblock[je];
 426                 }
 427             }
 428             *m = im;
 429         }
 430         if (idiscl > 0 || idiscu > 0) {
 431
 432             if (idiscl > 0) {
 433                 wkill = wu;
 434                 i__1 = idiscl;
 435                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 436                     iw = 0;
 437                     i__2 = *m;
 438                     for (je = 1; je <= i__2; ++je) {
 439                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
 440                             iw = je;
 441                             wkill = w[je];
 442                         }
 443                     }
 444                     iblock[iw] = 0;
 445                 }
 446             }
 447             if (idiscu > 0) {
 448
 449                 wkill = wl;
 450                 i__1 = idiscu;
 451                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 452                     iw = 0;
 453                     i__2 = *m;
 454                     for (je = 1; je <= i__2; ++je) {
 455                         if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
 456                             iw = je;
 457                             wkill = w[je];
 458                         }
 459                     }
 460                     iblock[iw] = 0;
 461                 }
 462             }
 463             im = 0;
 464             i__1 = *m;
 465             for (je = 1; je <= i__1; ++je) {
 466                 if (iblock[je] != 0) {
 467                     ++im;
 468                     w[im] = w[je];
 469                     iblock[im] = iblock[je];
 470                 }
 471             }
 472             *m = im;
 473         }
 474         if (idiscl < 0 || idiscu < 0) {
 475             toofew = 1;
 476         }
 477     }
 478
 479     if (iorder == 1 && *nsplit > 1) {
 480         i__1 = *m - 1;
 481         for (je = 1; je <= i__1; ++je) {
 482             ie = 0;
 483             tmp1 = w[je];
 484             i__2 = *m;
 485             for (j = je + 1; j <= i__2; ++j) {
 486                 if (w[j] < tmp1) {
 487                     ie = j;
 488                     tmp1 = w[j];
 489                 }
 490             }
 491
 492             if (ie != 0) {
 493                 itmp1 = iblock[ie];
 494                 w[ie] = w[je];
 495                 iblock[ie] = iblock[je];
 496                 w[je] = tmp1;
 497                 iblock[je] = itmp1;
 498             }
 499         }
 500     }
 501
 502     *info = 0;
 503     if (ncnvrg) {
 504         ++(*info);
 505     }
 506     if (toofew) {
 507         *info += 2;
 508     }
 509     return;
 510
 511 }
 512
 513