src/gmxlib/gmx_lapack/dlarrex.c

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