src/gmxlib/gmx_lapack/slasq2.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 #ifdef _MSC_VER
  42 #pragma warning(disable: 4723) /*division by zero - is used on purpose here*/
  43 #endif
  44
  45 void
  46 F77_FUNC(slasq2,SLASQ2)(int *n,
  47                         float *z__,
  48                         int *info)
  49 {
  50     int i__1, i__2, i__3;
  51     float d__1, d__2;
  52
  53     float d__, e;
  54     int k;
  55     float s, t;
  56     int i0, i4, n0, pp;
  57     float dee, eps, tol;
  58     int ipn4;
  59     float tol2;
  60     int ieee;
  61     int nbig;
  62     float dmin__, emin, emax;
  63     int kmin, ndiv, iter;
  64     float qmin, temp, qmax, zmax;
  65     int splt, nfail;
  66     float desig, trace, sigma;
  67     int iinfo;
  68     float deemin;
  69     int iwhila, iwhilb;
  70     float oldemn, safmin, minval;
  71     float posinf,neginf,negzro,newzro;
  72     float zero = 0.0;
  73     float one = 1.0;
  74
  75     --z__;
  76
  77     *info = 0;
  78     eps = GMX_FLOAT_EPS;
  79     minval = GMX_FLOAT_MIN;
  80     safmin = minval*(1.0+eps);
  81
  82     tol = eps * 100.;
  83
  84     d__1 = tol;
  85     tol2 = d__1 * d__1;
  86
  87     if (*n < 0) {
  88         *info = -1;
  89         return;
  90     } else if (*n == 0) {
  91         return;
  92     } else if (*n == 1) {
  93
  94         if (z__[1] < 0.) {
  95             *info = -201;
  96         }
  97         return;
  98     } else if (*n == 2) {
  99
 100         if (z__[2] < 0. || z__[3] < 0.) {
 101             *info = -2;
 102             return;
 103         } else if (z__[3] > z__[1]) {
 104             d__ = z__[3];
 105             z__[3] = z__[1];
 106             z__[1] = d__;
 107         }
 108         z__[5] = z__[1] + z__[2] + z__[3];
 109         if (z__[2] > z__[3] * tol2) {
 110             t = (z__[1] - z__[3] + z__[2]) * .5;
 111             s = z__[3] * (z__[2] / t);
 112             if (s <= t) {
 113                 s = z__[3] * (z__[2] / (t * (sqrt(s / t + 1.) + 1.)));
 114             } else {
 115                 s = z__[3] * (z__[2] / (t + sqrt(t) * sqrt(t + s)));
 116             }
 117             t = z__[1] + (s + z__[2]);
 118             z__[3] *= z__[1] / t;
 119             z__[1] = t;
 120         }
 121         z__[2] = z__[3];
 122         z__[6] = z__[2] + z__[1];
 123         return;
 124     }
 125     z__[*n * 2] = 0.;
 126     emin = z__[2];
 127     qmax = 0.;
 128     zmax = 0.;
 129     d__ = 0.;
 130     e = 0.;
 131
 132     i__1 = 2*(*n - 1);
 133     for (k = 1; k <= i__1; k += 2) {
 134         if (z__[k] < 0.) {
 135             *info = -(k + 200);
 136             return;
 137         } else if (z__[k + 1] < 0.) {
 138             *info = -(k + 201);
 139             return;
 140         }
 141         d__ += z__[k];
 142         e += z__[k + 1];
 143         d__1 = qmax, d__2 = z__[k];
 144         qmax = (d__1>d__2) ? d__1 : d__2;
 145         d__1 = emin, d__2 = z__[k + 1];
 146         emin = (d__1<d__2) ? d__1 : d__2;
 147         d__1 = (qmax>zmax) ? qmax : zmax;
 148         d__2 = z__[k + 1];
 149         zmax = (d__1>d__2) ? d__1 : d__2;
 150     }
 151     if (z__[(*n << 1) - 1] < 0.) {
 152         *info = -((*n << 1) + 199);
 153         return;
 154     }
 155     d__ += z__[(*n << 1) - 1];
 156     d__1 = qmax, d__2 = z__[(*n << 1) - 1];
 157     qmax = (d__1>d__2) ? d__1 : d__2;
 158     zmax = (qmax>zmax) ? qmax : zmax;
 159
 160     if (fabs(e)<GMX_FLOAT_MIN) {
 161         i__1 = *n;
 162         for (k = 2; k <= i__1; ++k) {
 163             z__[k] = z__[(k << 1) - 1];
 164         }
 165         F77_FUNC(slasrt,SLASRT)("D", n, &z__[1], &iinfo);
 166         z__[(*n << 1) - 1] = d__;
 167         return;
 168     }
 169
 170     trace = d__ + e;
 171
 172     if (fabs(trace)<GMX_FLOAT_MIN) {
 173         z__[(*n << 1) - 1] = 0.;
 174         return;
 175     }
 176
 177     ieee = 1;
 178     posinf = one/zero;
 179     if(posinf<=1.0)
 180       ieee = 0;
 181     neginf = -one/zero;
 182     if(neginf>=0.0)
 183       ieee = 0;
 184     negzro = one/(neginf+one);
 185     if(fabs(negzro)>GMX_FLOAT_MIN)
 186       ieee = 0;
 187     neginf = one/negzro;
 188     if(neginf>=0)
 189       ieee = 0;
 190     newzro = negzro + zero;
 191     if(fabs(newzro-zero)>GMX_FLOAT_MIN)
 192       ieee = 0;
 193     posinf = one /newzro;
 194     if(posinf<=one)
 195       ieee = 0;
 196     neginf = neginf*posinf;
 197     if(neginf>=zero)
 198       ieee = 0;
 199     posinf = posinf*posinf;
 200     if(posinf<=1.0)
 201       ieee = 0;
 202
 203     for (k = *n << 1; k >= 2; k += -2) {
 204         z__[k * 2] = 0.;
 205         z__[(k << 1) - 1] = z__[k];
 206         z__[(k << 1) - 2] = 0.;
 207         z__[(k << 1) - 3] = z__[k - 1];
 208     }
 209
 210     i0 = 1;
 211     n0 = *n;
 212
 213     if (z__[(i0 << 2) - 3] * 1.5 < z__[(n0 << 2) - 3]) {
 214         ipn4 = 4*(i0 + n0);
 215         i__1 = 2*(i0 + n0 - 1);
 216         for (i4 = i0 << 2; i4 <= i__1; i4 += 4) {
 217             temp = z__[i4 - 3];
 218             z__[i4 - 3] = z__[ipn4 - i4 - 3];
 219             z__[ipn4 - i4 - 3] = temp;
 220             temp = z__[i4 - 1];
 221             z__[i4 - 1] = z__[ipn4 - i4 - 5];
 222             z__[ipn4 - i4 - 5] = temp;
 223         }
 224     }
 225
 226     pp = 0;
 227
 228     for (k = 1; k <= 2; ++k) {
 229
 230         d__ = z__[(n0 << 2) + pp - 3];
 231         i__1 = (i0 << 2) + pp;
 232         for (i4 = 4*(n0 - 1) + pp; i4 >= i__1; i4 += -4) {
 233             if (z__[i4 - 1] <= tol2 * d__) {
 234                 z__[i4 - 1] = -0.;
 235                 d__ = z__[i4 - 3];
 236             } else {
 237                 d__ = z__[i4 - 3] * (d__ / (d__ + z__[i4 - 1]));
 238             }
 239         }
 240
 241         emin = z__[(i0 << 2) + pp + 1];
 242         d__ = z__[(i0 << 2) + pp - 3];
 243         i__1 = 4*(n0 - 1) + pp;
 244         for (i4 = (i0 << 2) + pp; i4 <= i__1; i4 += 4) {
 245             z__[i4 - (pp << 1) - 2] = d__ + z__[i4 - 1];
 246             if (z__[i4 - 1] <= tol2 * d__) {
 247                 z__[i4 - 1] = -0.;
 248                 z__[i4 - (pp << 1) - 2] = d__;
 249                 z__[i4 - (pp << 1)] = 0.;
 250                 d__ = z__[i4 + 1];
 251             } else if (safmin * z__[i4 + 1] < z__[i4 - (pp << 1) - 2] &&
 252                     safmin * z__[i4 - (pp << 1) - 2] < z__[i4 + 1]) {
 253                 temp = z__[i4 + 1] / z__[i4 - (pp << 1) - 2];
 254                 z__[i4 - (pp << 1)] = z__[i4 - 1] * temp;
 255                 d__ *= temp;
 256             } else {
 257                 z__[i4 - (pp << 1)] = z__[i4 + 1] * (z__[i4 - 1] / z__[i4 - (
 258                         pp << 1) - 2]);
 259                 d__ = z__[i4 + 1] * (d__ / z__[i4 - (pp << 1) - 2]);
 260             }
 261             d__1 = emin, d__2 = z__[i4 - (pp << 1)];
 262             emin = (d__1<d__2) ? d__1 : d__2;
 263         }
 264         z__[(n0 << 2) - pp - 2] = d__;
 265
 266
 267         qmax = z__[(i0 << 2) - pp - 2];
 268         i__1 = (n0 << 2) - pp - 2;
 269         for (i4 = (i0 << 2) - pp + 2; i4 <= i__1; i4 += 4) {
 270             d__1 = qmax, d__2 = z__[i4];
 271             qmax = (d__1>d__2) ? d__1 : d__2;
 272         }
 273
 274         pp = 1 - pp;
 275     }
 276
 277     iter = 2;
 278     nfail = 0;
 279     ndiv = 2*(n0 - i0);
 280
 281     i__1 = *n + 1;
 282     for (iwhila = 1; iwhila <= i__1; ++iwhila) {
 283         if (n0 < 1) {
 284             goto L170;
 285         }
 286
 287         desig = 0.;
 288         if (n0 == *n) {
 289             sigma = 0.;
 290         } else {
 291             sigma = -z__[(n0 << 2) - 1];
 292         }
 293         if (sigma < 0.) {
 294             *info = 1;
 295             return;
 296         }
 297
 298         emax = 0.;
 299         if (n0 > i0) {
 300             emin = fabs(z__[(n0 << 2) - 5]);
 301         } else {
 302             emin = 0.;
 303         }
 304         qmin = z__[(n0 << 2) - 3];
 305         qmax = qmin;
 306         for (i4 = n0 << 2; i4 >= 8; i4 += -4) {
 307             if (z__[i4 - 5] <= 0.) {
 308                 goto L100;
 309             }
 310             if (qmin >= emax * 4.) {
 311                 d__1 = qmin, d__2 = z__[i4 - 3];
 312                 qmin = (d__1<d__2) ? d__1 : d__2;
 313                 d__1 = emax, d__2 = z__[i4 - 5];
 314                 emax = (d__1>d__2) ? d__1 : d__2;
 315             }
 316             d__1 = qmax, d__2 = z__[i4 - 7] + z__[i4 - 5];
 317             qmax = (d__1>d__2) ? d__1 : d__2;
 318             d__1 = emin, d__2 = z__[i4 - 5];
 319             emin = (d__1<d__2) ? d__1 : d__2;
 320         }
 321         i4 = 4;
 322
 323 L100:
 324         i0 = i4 / 4;
 325         pp = 0;
 326
 327         if (n0 - i0 > 1) {
 328             dee = z__[(i0 << 2) - 3];
 329             deemin = dee;
 330             kmin = i0;
 331             i__2 = (n0 << 2) - 3;
 332             for (i4 = (i0 << 2) - 3; i4 <= i__2; i4 += 4) {
 333                 dee = z__[i4] * (dee / (dee + z__[i4 - 2]));
 334                 if (dee <= deemin) {
 335                     deemin = dee;
 336                     kmin = (i4 + 3) / 4;
 337                 }
 338             }
 339             if (2*(kmin - i0) < n0 - kmin && deemin <= z__[(n0 << 2) - 3] *
 340                     .5) {
 341                 ipn4 = 4*(i0 + n0);
 342                 pp = 2;
 343                 i__2 = 2*(i0 + n0 - 1);
 344                 for (i4 = i0 << 2; i4 <= i__2; i4 += 4) {
 345                     temp = z__[i4 - 3];
 346                     z__[i4 - 3] = z__[ipn4 - i4 - 3];
 347                     z__[ipn4 - i4 - 3] = temp;
 348                     temp = z__[i4 - 2];
 349                     z__[i4 - 2] = z__[ipn4 - i4 - 2];
 350                     z__[ipn4 - i4 - 2] = temp;
 351                     temp = z__[i4 - 1];
 352                     z__[i4 - 1] = z__[ipn4 - i4 - 5];
 353                     z__[ipn4 - i4 - 5] = temp;
 354                     temp = z__[i4];
 355                     z__[i4] = z__[ipn4 - i4 - 4];
 356                     z__[ipn4 - i4 - 4] = temp;
 357                 }
 358             }
 359         }
 360
 361
 362         d__1 = 0., d__2 = qmin - sqrt(qmin) * 2. * sqrt(emax);
 363         dmin__ = -((d__1>d__2) ? d__1 : d__2);
 364
 365         nbig = (n0 - i0 + 1) * 30;
 366         i__2 = nbig;
 367         for (iwhilb = 1; iwhilb <= i__2; ++iwhilb) {
 368             if (i0 > n0) {
 369                 goto L150;
 370             }
 371
 372             F77_FUNC(slasq3,SLASQ3)(&i0, &n0, &z__[1], &pp, &dmin__, &sigma, &desig, &qmax, &
 373                     nfail, &iter, &ndiv, &ieee);
 374
 375             pp = 1 - pp;
 376
 377             if (pp == 0 && n0 - i0 >= 3) {
 378                 if (z__[n0 * 4] <= tol2 * qmax || z__[(n0 << 2) - 1] <= tol2 *
 379                          sigma) {
 380                     splt = i0 - 1;
 381                     qmax = z__[(i0 << 2) - 3];
 382                     emin = z__[(i0 << 2) - 1];
 383                     oldemn = z__[i0 * 4];
 384                     i__3 = 4*(n0 - 3);
 385                     for (i4 = i0 << 2; i4 <= i__3; i4 += 4) {
 386                         if (z__[i4] <= tol2 * z__[i4 - 3] || z__[i4 - 1] <=
 387                                 tol2 * sigma) {
 388                             z__[i4 - 1] = -sigma;
 389                             splt = i4 / 4;
 390                             qmax = 0.;
 391                             emin = z__[i4 + 3];
 392                             oldemn = z__[i4 + 4];
 393                         } else {
 394                             d__1 = qmax, d__2 = z__[i4 + 1];
 395                             qmax = (d__1>d__2) ? d__1 : d__2;
 396                             d__1 = emin, d__2 = z__[i4 - 1];
 397                             emin = (d__1<d__2) ? d__1 : d__2;
 398                             d__1 = oldemn, d__2 = z__[i4];
 399                             oldemn = (d__1<d__2) ? d__1 : d__2;
 400                         }
 401                     }
 402                     z__[(n0 << 2) - 1] = emin;
 403                     z__[n0 * 4] = oldemn;
 404                     i0 = splt + 1;
 405                 }
 406             }
 407         }
 408
 409         *info = 2;
 410         return;
 411
 412 L150:
 413         ;
 414     }
 415
 416     *info = 3;
 417     return;
 418
 419
 420 L170:
 421
 422     i__1 = *n;
 423     for (k = 2; k <= i__1; ++k) {
 424         z__[k] = z__[(k << 2) - 3];
 425     }
 426
 427     F77_FUNC(slasrt,SLASRT)("D", n, &z__[1], &iinfo);
 428
 429     e = 0.;
 430     for (k = *n; k >= 1; --k) {
 431         e += z__[k];
 432     }
 433
 434
 435     z__[(*n << 1) + 1] = trace;
 436     z__[(*n << 1) + 2] = e;
 437     z__[(*n << 1) + 3] = (float) iter;
 438     i__1 = *n;
 439     z__[(*n << 1) + 4] = (float) ndiv / (float) (i__1 * i__1);
 440     z__[(*n << 1) + 5] = nfail * 100. / (float) iter;
 441
 442     return;
 443
 444 }
 445
 446
 447