src/gromacs/linearalgebra/gmx_lapack/dlasq2.c

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