src/gromacs/linearalgebra/gmx_lapack/ssterf.c

   1 #include <math.h>
   2 #include "../gmx_lapack.h"
   3 #include "lapack_limits.h"
   4
   5 #include "types/simple.h"
   6
   7 void
   8 F77_FUNC(ssterf,SSTERF)(int *n,
   9         float *d__,
  10         float *e,
  11         int *info)
  12 {
  13     int i__1;
  14     float d__1;
  15
  16     float c__;
  17     int i__, l, m;
  18     float p, r__, s;
  19     int l1;
  20     float bb, rt1, rt2, eps, rte;
  21     int lsv;
  22     float eps2, oldc;
  23     int lend, jtot;
  24     float gamma, alpha, sigma, anorm;
  25       int iscale;
  26     float oldgam;
  27     float safmax;
  28     int lendsv;
  29     float ssfmin;
  30     int nmaxit;
  31     float ssfmax;
  32     int c__0 = 0;
  33     int c__1 = 1;
  34     float c_b32 = 1.;
  35     const float safmin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
  36
  37     --e;
  38     --d__;
  39
  40     *info = 0;
  41
  42     if (*n < 0) {
  43         *info = -1;
  44         i__1 = -(*info);
  45         return;
  46     }
  47     if (*n <= 1) {
  48         return;
  49     }
  50
  51     eps = GMX_FLOAT_EPS;
  52     d__1 = eps;
  53     eps2 = d__1 * d__1;
  54     safmax = 1. / safmin;
  55     ssfmax = sqrt(safmax) / 3.;
  56     ssfmin = sqrt(safmin) / eps2;
  57
  58     nmaxit = *n * 30;
  59     sigma = 0.;
  60     jtot = 0;
  61
  62     l1 = 1;
  63
  64 L10:
  65     if (l1 > *n) {
  66       F77_FUNC(slasrt,SLASRT)("I", n, &d__[1], info);
  67       return;
  68     }
  69     if (l1 > 1) {
  70         e[l1 - 1] = 0.;
  71     }
  72     i__1 = *n - 1;
  73     for (m = l1; m <= i__1; ++m) {
  74         if (fabs(e[m]) <= sqrt(fabs(d__[m])) *
  75                 sqrt(fabs(d__[m + 1])) * eps) {
  76             e[m] = 0.;
  77             goto L30;
  78         }
  79     }
  80     m = *n;
  81
  82 L30:
  83     l = l1;
  84     lsv = l;
  85     lend = m;
  86     lendsv = lend;
  87     l1 = m + 1;
  88     if (lend == l) {
  89         goto L10;
  90     }
  91
  92     i__1 = lend - l + 1;
  93     anorm = F77_FUNC(slanst,SLANST)("I", &i__1, &d__[l], &e[l]);
  94     iscale = 0;
  95     if (anorm > ssfmax) {
  96         iscale = 1;
  97         i__1 = lend - l + 1;
  98         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
  99                 info);
 100         i__1 = lend - l;
 101         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
 102                 info);
 103     } else if (anorm < ssfmin) {
 104         iscale = 2;
 105         i__1 = lend - l + 1;
 106         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
 107                 info);
 108         i__1 = lend - l;
 109         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
 110                 info);
 111     }
 112
 113     i__1 = lend - 1;
 114     for (i__ = l; i__ <= i__1; ++i__) {
 115         d__1 = e[i__];
 116         e[i__] = d__1 * d__1;
 117     }
 118
 119     if (fabs(d__[lend]) < fabs(d__[l])) {
 120         lend = lsv;
 121         l = lendsv;
 122     }
 123
 124     if (lend >= l) {
 125
 126 L50:
 127         if (l != lend) {
 128             i__1 = lend - 1;
 129             for (m = l; m <= i__1; ++m) {
 130                 if (fabs(e[m]) <= eps2 * fabs(d__[m] * d__[m + 1])) {
 131                     goto L70;
 132                 }
 133             }
 134         }
 135         m = lend;
 136
 137 L70:
 138         if (m < lend) {
 139             e[m] = 0.;
 140         }
 141         p = d__[l];
 142         if (m == l) {
 143             goto L90;
 144         }
 145         if (m == l + 1) {
 146             rte = sqrt(e[l]);
 147             F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
 148             d__[l] = rt1;
 149             d__[l + 1] = rt2;
 150             e[l] = 0.;
 151             l += 2;
 152             if (l <= lend) {
 153                 goto L50;
 154             }
 155             goto L150;
 156         }
 157
 158         if (jtot == nmaxit) {
 159             goto L150;
 160         }
 161         ++jtot;
 162
 163         rte = sqrt(e[l]);
 164         sigma = (d__[l + 1] - p) / (rte * 2.);
 165         r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
 166         sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
 167
 168         c__ = 1.;
 169         s = 0.;
 170         gamma = d__[m] - sigma;
 171         p = gamma * gamma;
 172
 173         i__1 = l;
 174         for (i__ = m - 1; i__ >= i__1; --i__) {
 175             bb = e[i__];
 176             r__ = p + bb;
 177             if (i__ != m - 1) {
 178                 e[i__ + 1] = s * r__;
 179             }
 180             oldc = c__;
 181             c__ = p / r__;
 182             s = bb / r__;
 183             oldgam = gamma;
 184             alpha = d__[i__];
 185             gamma = c__ * (alpha - sigma) - s * oldgam;
 186             d__[i__ + 1] = oldgam + (alpha - gamma);
 187             if (fabs(c__)>GMX_FLOAT_MIN) {
 188                 p = gamma * gamma / c__;
 189             } else {
 190                 p = oldc * bb;
 191             }
 192         }
 193
 194         e[l] = s * p;
 195         d__[l] = sigma + gamma;
 196         goto L50;
 197
 198 L90:
 199         d__[l] = p;
 200
 201         ++l;
 202         if (l <= lend) {
 203             goto L50;
 204         }
 205         goto L150;
 206
 207     } else {
 208
 209 L100:
 210         i__1 = lend + 1;
 211         for (m = l; m >= i__1; --m) {
 212             if (fabs(e[m - 1]) <= eps2 * fabs(d__[m] * d__[m - 1])) {
 213                 goto L120;
 214             }
 215         }
 216         m = lend;
 217
 218 L120:
 219         if (m > lend) {
 220             e[m - 1] = 0.;
 221         }
 222         p = d__[l];
 223         if (m == l) {
 224             goto L140;
 225         }
 226
 227         if (m == l - 1) {
 228             rte = sqrt(e[l - 1]);
 229             F77_FUNC(slae2,SLAE2)(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
 230             d__[l] = rt1;
 231             d__[l - 1] = rt2;
 232             e[l - 1] = 0.;
 233             l += -2;
 234             if (l >= lend) {
 235                 goto L100;
 236             }
 237             goto L150;
 238         }
 239
 240         if (jtot == nmaxit) {
 241             goto L150;
 242         }
 243         ++jtot;
 244
 245         rte = sqrt(e[l - 1]);
 246         sigma = (d__[l - 1] - p) / (rte * 2.);
 247         r__ = F77_FUNC(slapy2,SLAPY2)(&sigma, &c_b32);
 248         sigma = p - rte / (sigma + ( (sigma>0) ? r__ : -r__));
 249
 250         c__ = 1.;
 251         s = 0.;
 252         gamma = d__[m] - sigma;
 253         p = gamma * gamma;
 254
 255         i__1 = l - 1;
 256         for (i__ = m; i__ <= i__1; ++i__) {
 257             bb = e[i__];
 258             r__ = p + bb;
 259             if (i__ != m) {
 260                 e[i__ - 1] = s * r__;
 261             }
 262             oldc = c__;
 263             c__ = p / r__;
 264             s = bb / r__;
 265             oldgam = gamma;
 266             alpha = d__[i__ + 1];
 267             gamma = c__ * (alpha - sigma) - s * oldgam;
 268             d__[i__] = oldgam + (alpha - gamma);
 269             if (fabs(c__)>GMX_FLOAT_MIN) {
 270                 p = gamma * gamma / c__;
 271             } else {
 272                 p = oldc * bb;
 273             }
 274         }
 275
 276         e[l - 1] = s * p;
 277         d__[l] = sigma + gamma;
 278         goto L100;
 279
 280 L140:
 281         d__[l] = p;
 282
 283         --l;
 284         if (l >= lend) {
 285             goto L100;
 286         }
 287         goto L150;
 288
 289     }
 290
 291 L150:
 292     if (iscale == 1) {
 293         i__1 = lendsv - lsv + 1;
 294         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
 295                 n, info);
 296     }
 297     if (iscale == 2) {
 298         i__1 = lendsv - lsv + 1;
 299         F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
 300                 n, info);
 301     }
 302
 303     if (jtot < nmaxit) {
 304         goto L10;
 305     }
 306     i__1 = *n - 1;
 307     for (i__ = 1; i__ <= i__1; ++i__) {
 308         if (fabs(e[i__])>GMX_FLOAT_MIN) {
 309             ++(*info);
 310         }
 311     }
 312     return;
 313 }
 314
 315