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