src/gmxlib/gmx_lapack/slar1vx.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
  37 #include <types/simple.h>
  38
  39 #include "gmx_lapack.h"
  40 #include "lapack_limits.h"
  41
  42 void F77_FUNC(slar1vx,SLAR1VX)(int *n,
  43               int *b1,
  44               int *bn,
  45               float *sigma,
  46               float *d__,
  47               float *l,
  48               float *ld,
  49               float *lld,
  50               float *eval,
  51               float *gersch,
  52               float *z__,
  53               float *ztz,
  54               float *mingma,
  55               int *r__,
  56               int *isuppz,
  57               float *work)
  58 {
  59     int i__1;
  60
  61     int i__, j;
  62     float s;
  63     int r1, r2;
  64     int to;
  65     float eps, tmp;
  66     int indp, inds, from;
  67     float dplus;
  68     int sawnan;
  69     int indumn;
  70     float dminus;
  71
  72     --work;
  73     --isuppz;
  74     --z__;
  75     --gersch;
  76     --lld;
  77     --ld;
  78     --l;
  79     --d__;
  80
  81     /* Function Body */
  82     eps = GMX_FLOAT_EPS;
  83     if (*r__ == 0) {
  84
  85         r1 = *b1;
  86         r2 = *bn;
  87         i__1 = *bn;
  88         for (i__ = *b1; i__ <= i__1; ++i__) {
  89             if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
  90                 r1 = i__;
  91                 goto L20;
  92             }
  93         }
  94         goto L40;
  95 L20:
  96         i__1 = *b1;
  97         for (i__ = *bn; i__ >= i__1; --i__) {
  98             if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
  99                 r2 = i__;
 100                 goto L40;
 101             }
 102         }
 103     } else {
 104         r1 = *r__;
 105         r2 = *r__;
 106     }
 107
 108 L40:
 109     indumn = *n;
 110     inds = (*n << 1) + 1;
 111     indp = *n * 3 + 1;
 112     sawnan = 0;
 113
 114     if (*b1 == 1) {
 115         work[inds] = 0.;
 116     } else {
 117         work[inds] = lld[*b1 - 1];
 118     }
 119     s = work[inds] - *sigma;
 120     i__1 = r2 - 1;
 121     for (i__ = *b1; i__ <= i__1; ++i__) {
 122         dplus = d__[i__] + s;
 123         work[i__] = ld[i__] / dplus;
 124         work[inds + i__] = s * work[i__] * l[i__];
 125         s = work[inds + i__] - *sigma;
 126     }
 127
 128     if (! (s > 0. || s < 1.)) {
 129
 130         sawnan = 1;
 131         j = *b1 + 1;
 132 L60:
 133         if (work[inds + j] > 0. || work[inds + j] < 1.) {
 134             ++j;
 135             goto L60;
 136         }
 137         work[inds + j] = lld[j];
 138         s = work[inds + j] - *sigma;
 139         i__1 = r2 - 1;
 140         for (i__ = j + 1; i__ <= i__1; ++i__) {
 141             dplus = d__[i__] + s;
 142             work[i__] = ld[i__] / dplus;
 143             if (fabs(work[i__])<GMX_FLOAT_MIN) {
 144                 work[inds + i__] = lld[i__];
 145             } else {
 146                 work[inds + i__] = s * work[i__] * l[i__];
 147             }
 148             s = work[inds + i__] - *sigma;
 149         }
 150     }
 151
 152     work[indp + *bn - 1] = d__[*bn] - *sigma;
 153     i__1 = r1;
 154     for (i__ = *bn - 1; i__ >= i__1; --i__) {
 155         dminus = lld[i__] + work[indp + i__];
 156         tmp = d__[i__] / dminus;
 157         work[indumn + i__] = l[i__] * tmp;
 158         work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
 159     }
 160     tmp = work[indp + r1 - 1];
 161     if (! (tmp > 0. || tmp < 1.)) {
 162
 163         sawnan = 1;
 164         j = *bn - 3;
 165 L90:
 166         if (work[indp + j] > 0. || work[indp + j] < 1.) {
 167             --j;
 168             goto L90;
 169         }
 170         work[indp + j] = d__[j + 1] - *sigma;
 171         i__1 = r1;
 172         for (i__ = j; i__ >= i__1; --i__) {
 173             dminus = lld[i__] + work[indp + i__];
 174             tmp = d__[i__] / dminus;
 175             work[indumn + i__] = l[i__] * tmp;
 176             if (fabs(tmp)<GMX_FLOAT_MIN) {
 177                 work[indp + i__ - 1] = d__[i__] - *sigma;
 178             } else {
 179                 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
 180             }
 181         }
 182     }
 183
 184     *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
 185     if (fabs(*mingma)<GMX_FLOAT_MIN) {
 186         *mingma = eps * work[inds + r1 - 1];
 187     }
 188     *r__ = r1;
 189     i__1 = r2 - 1;
 190     for (i__ = r1; i__ <= i__1; ++i__) {
 191         tmp = work[inds + i__] + work[indp + i__];
 192         if (fabs(tmp)<GMX_FLOAT_MIN) {
 193             tmp = eps * work[inds + i__];
 194         }
 195         if (fabs(tmp) < fabs(*mingma)) {
 196             *mingma = tmp;
 197             *r__ = i__ + 1;
 198         }
 199     }
 200
 201     isuppz[1] = *b1;
 202     isuppz[2] = *bn;
 203     z__[*r__] = 1.;
 204     *ztz = 1.;
 205     if (! sawnan) {
 206         from = *r__ - 1;
 207         i__1 = *r__ - 32;
 208         to = (i__1>(*b1)) ? i__1 : (*b1);
 209 L120:
 210         if (from >= *b1) {
 211             i__1 = to;
 212             for (i__ = from; i__ >= i__1; --i__) {
 213                 z__[i__] = -(work[i__] * z__[i__ + 1]);
 214                 *ztz += z__[i__] * z__[i__];
 215             }
 216             if (fabs(z__[to]) <= eps && fabs(z__[to + 1]) <= eps) {
 217                 isuppz[1] = to + 2;
 218             } else {
 219                 from = to - 1;
 220                 i__1 = to - 32;
 221                 to = (i__1>*b1) ? i__1 : *b1;
 222                 goto L120;
 223             }
 224         }
 225         from = *r__ + 1;
 226         i__1 = *r__ + 32;
 227         to = (i__1<*bn) ? i__1 : *bn;
 228 L140:
 229         if (from <= *bn) {
 230             i__1 = to;
 231             for (i__ = from; i__ <= i__1; ++i__) {
 232                 z__[i__] = -(work[indumn + i__ - 1] * z__[i__ - 1]);
 233                 *ztz += z__[i__] * z__[i__];
 234             }
 235             if (fabs(z__[to]) <= eps && fabs(z__[to - 1]) <= eps) {
 236                 isuppz[2] = to - 2;
 237             } else {
 238                 from = to + 1;
 239                 i__1 = to + 32;
 240                 to = (i__1<*bn) ? i__1 : *bn;
 241                 goto L140;
 242             }
 243         }
 244     } else {
 245         i__1 = *b1;
 246         for (i__ = *r__ - 1; i__ >= i__1; --i__) {
 247             if (fabs(z__[i__ + 1])<GMX_FLOAT_MIN) {
 248                 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
 249             } else {
 250                 z__[i__] = -(work[i__] * z__[i__ + 1]);
 251             }
 252             if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
 253                 isuppz[1] = i__ + 2;
 254                 goto L170;
 255             }
 256             *ztz += z__[i__] * z__[i__];
 257         }
 258 L170:
 259         i__1 = *bn - 1;
 260         for (i__ = *r__; i__ <= i__1; ++i__) {
 261             if (fabs(z__[i__])<GMX_FLOAT_MIN) {
 262                 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
 263             } else {
 264                 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
 265             }
 266             if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
 267                 isuppz[2] = i__ - 1;
 268                 break;
 269             }
 270             *ztz += z__[i__ + 1] * z__[i__ + 1];
 271         }
 272     }
 273
 274     return;
 275
 276 }
 277
 278