src/gmxlib/gmx_lapack/slasq3.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 <types/simple.h>
  37
  38 #include "gmx_lapack.h"
  39 #include "lapack_limits.h"
  40
  41 void
  42 F77_FUNC(slasq3,SLASQ3)(int *i0,
  43                         int *n0,
  44                         float *z__,
  45                         int *pp,
  46                         float *dmin__,
  47                         float *sigma,
  48                         float *desig,
  49                         float *qmax,
  50                         int *nfail,
  51                         int *iter,
  52                         int *ndiv,
  53         int *ieee)
  54 {
  55
  56     int ttype = 0;
  57     float dmin1 = 0.;
  58     float dmin2 = 0.;
  59     float dn = 0.;
  60     float dn1 = 0.;
  61     float dn2 = 0.;
  62     float tau = 0.;
  63
  64     int i__1;
  65     float d__1, d__2;
  66     float s, t;
  67     int j4, nn;
  68     float eps, tol;
  69     int n0in, ipn4;
  70     float tol2, temp;
  71     --z__;
  72
  73     n0in = *n0;
  74     eps = GMX_FLOAT_EPS;
  75     tol = eps * 100.;
  76     d__1 = tol;
  77     tol2 = d__1 * d__1;
  78
  79
  80 L10:
  81
  82     if (*n0 < *i0) {
  83         return;
  84     }
  85     if (*n0 == *i0) {
  86         goto L20;
  87     }
  88     nn = (*n0 << 2) + *pp;
  89     if (*n0 == *i0 + 1) {
  90         goto L40;
  91     }
  92
  93     if (z__[nn - 5] > tol2 * (*sigma + z__[nn - 3]) && z__[nn - (*pp << 1) -
  94             4] > tol2 * z__[nn - 7]) {
  95         goto L30;
  96     }
  97
  98 L20:
  99
 100     z__[(*n0 << 2) - 3] = z__[(*n0 << 2) + *pp - 3] + *sigma;
 101     --(*n0);
 102     goto L10;
 103
 104 L30:
 105
 106     if (z__[nn - 9] > tol2 * *sigma && z__[nn - (*pp << 1) - 8] > tol2 * z__[
 107             nn - 11]) {
 108         goto L50;
 109     }
 110
 111 L40:
 112
 113     if (z__[nn - 3] > z__[nn - 7]) {
 114         s = z__[nn - 3];
 115         z__[nn - 3] = z__[nn - 7];
 116         z__[nn - 7] = s;
 117     }
 118     if (z__[nn - 5] > z__[nn - 3] * tol2) {
 119         t = (z__[nn - 7] - z__[nn - 3] + z__[nn - 5]) * .5;
 120         s = z__[nn - 3] * (z__[nn - 5] / t);
 121         if (s <= t) {
 122             s = z__[nn - 3] * (z__[nn - 5] / (t * (sqrt(s / t + 1.) + 1.)));
 123         } else {
 124             s = z__[nn - 3] * (z__[nn - 5] / (t + sqrt(t) * sqrt(t + s)));
 125         }
 126         t = z__[nn - 7] + (s + z__[nn - 5]);
 127         z__[nn - 3] *= z__[nn - 7] / t;
 128         z__[nn - 7] = t;
 129     }
 130     z__[(*n0 << 2) - 7] = z__[nn - 7] + *sigma;
 131     z__[(*n0 << 2) - 3] = z__[nn - 3] + *sigma;
 132     *n0 += -2;
 133     goto L10;
 134
 135 L50:
 136     if (*pp == 2) {
 137         *pp = 0;
 138     }
 139
 140     if (*dmin__ <= 0. || *n0 < n0in) {
 141         if (z__[(*i0 << 2) + *pp - 3] * 1.5 < z__[(*n0 << 2) + *pp - 3]) {
 142             ipn4 = 4*(*i0 + *n0);
 143             i__1 = 2*(*i0 + *n0 - 1);
 144             for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
 145                 temp = z__[j4 - 3];
 146                 z__[j4 - 3] = z__[ipn4 - j4 - 3];
 147                 z__[ipn4 - j4 - 3] = temp;
 148                 temp = z__[j4 - 2];
 149                 z__[j4 - 2] = z__[ipn4 - j4 - 2];
 150                 z__[ipn4 - j4 - 2] = temp;
 151                 temp = z__[j4 - 1];
 152                 z__[j4 - 1] = z__[ipn4 - j4 - 5];
 153                 z__[ipn4 - j4 - 5] = temp;
 154                 temp = z__[j4];
 155                 z__[j4] = z__[ipn4 - j4 - 4];
 156                 z__[ipn4 - j4 - 4] = temp;
 157             }
 158             if (*n0 - *i0 <= 4) {
 159                 z__[(*n0 << 2) + *pp - 1] = z__[(*i0 << 2) + *pp - 1];
 160                 z__[(*n0 << 2) - *pp] = z__[(*i0 << 2) - *pp];
 161             }
 162             d__1 = dmin2, d__2 = z__[(*n0 << 2) + *pp - 1];
 163             dmin2 = ((d__1<d__2) ? d__1 : d__2);
 164             d__1 = z__[(*n0 << 2) + *pp - 1], d__2 = z__[(*i0 << 2) + *pp - 1]
 165                     , d__1 = ((d__1<d__2) ? d__1 : d__2), d__2 = z__[(*i0 << 2) + *pp + 3];
 166             z__[(*n0 << 2) + *pp - 1] = ((d__1<d__2) ? d__1 : d__2);
 167             d__1 = z__[(*n0 << 2) - *pp], d__2 = z__[(*i0 << 2) - *pp], d__1 =
 168                      ((d__1<d__2) ? d__1 : d__2), d__2 = z__[(*i0 << 2) - *pp + 4];
 169             z__[(*n0 << 2) - *pp] = ((d__1<d__2) ? d__1 : d__2);
 170             d__1 = *qmax;
 171             d__2 = z__[(*i0 << 2) + *pp - 3];
 172             d__1 = (d__1>d__2) ? d__1 : d__2;
 173             d__2 = z__[(*i0 << 2) + *pp + 1];
 174             *qmax = ((d__1>d__2) ? d__1 : d__2);
 175             *dmin__ = -0.;
 176         }
 177     }
 178
 179
 180     F77_FUNC(slasq4,SLASQ4)(i0, n0, &z__[1], pp, &n0in, dmin__, &dmin1, &dmin2, &dn, &dn1, &
 181             dn2, &tau, &ttype);
 182
 183 L70:
 184
 185     F77_FUNC(slasq5,SLASQ5)(i0, n0, &z__[1], pp, &tau, dmin__, &dmin1, &dmin2, &dn, &dn1, &
 186             dn2, ieee);
 187
 188     *ndiv += *n0 - *i0 + 2;
 189     ++(*iter);
 190
 191     if (*dmin__ >= 0. && dmin1 > 0.) {
 192
 193         goto L90;
 194
 195     } else if (*dmin__ < 0. && dmin1 > 0. && z__[4*(*n0 - 1) - *pp] < tol *
 196              (*sigma + dn1) && fabs(dn) < tol * *sigma) {
 197
 198         z__[4*(*n0 - 1) - *pp + 2] = 0.;
 199         *dmin__ = 0.;
 200         goto L90;
 201     } else if (*dmin__ < 0.) {
 202
 203         ++(*nfail);
 204         if (ttype < -22) {
 205
 206             tau = 0.;
 207         } else if (dmin1 > 0.) {
 208
 209             tau = (tau + *dmin__) * (1. - eps * 2.);
 210             ttype += -11;
 211         } else {
 212
 213             tau *= .25;
 214             ttype += -12;
 215         }
 216         goto L70;
 217     }
 218     else {
 219
 220         goto L80;
 221     }
 222
 223 L80:
 224     F77_FUNC(slasq6,SLASQ6)(i0, n0, &z__[1], pp, dmin__, &dmin1, &dmin2, &dn, &dn1, &dn2);
 225     *ndiv += *n0 - *i0 + 2;
 226     ++(*iter);
 227     tau = 0.;
 228
 229 L90:
 230     if (tau < *sigma) {
 231         *desig += tau;
 232         t = *sigma + *desig;
 233         *desig -= t - *sigma;
 234     } else {
 235         t = *sigma + tau;
 236         *desig = *sigma - (t - tau) + *desig;
 237     }
 238     *sigma = t;
 239
 240     return;
 241 }