src/gmxlib/gmx_lapack/dsteqr.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_blas.h"
  39 #include "gmx_lapack.h"
  40 #include "lapack_limits.h"
  41
  42 void
  43 F77_FUNC(dsteqr,DSTEQR)(const char *    compz,
  44         int *     n,
  45         double *  d__,
  46         double *  e,
  47         double *  z__,
  48         int *     ldz,
  49         double *  work,
  50         int *     info)
  51 {
  52     double c_b9 = 0.;
  53     double c_b10 = 1.;
  54     int c__0 = 0;
  55     int c__1 = 1;
  56     int c__2 = 2;
  57     int z_dim1, z_offset, i__1, i__2;
  58     double d__1, d__2;
  59
  60     double b, c__, f, g;
  61     int i__, j, k, l, m;
  62     double p, r__, s;
  63     int l1, ii, mm, lm1, mm1, nm1;
  64     double rt1, rt2, eps;
  65     int lsv;
  66     double tst, eps2;
  67     int lend, jtot;
  68     double anorm;
  69     int lendm1, lendp1;
  70     int iscale;
  71     double safmin,minval;
  72     double safmax;
  73     int lendsv;
  74     double ssfmin;
  75     int nmaxit, icompz;
  76     double ssfmax;
  77
  78
  79     --d__;
  80     --e;
  81     z_dim1 = *ldz;
  82     z_offset = 1 + z_dim1;
  83     z__ -= z_offset;
  84     --work;
  85
  86     *info = 0;
  87
  88     if (*compz=='N' || *compz=='n') {
  89         icompz = 0;
  90     } else if (*compz=='V' || *compz=='v') {
  91         icompz = 1;
  92     } else if (*compz=='I' || *compz=='i') {
  93         icompz = 2;
  94     } else {
  95         icompz = -1;
  96     }
  97     if (icompz < 0) {
  98         *info = -1;
  99     } else if (*n < 0) {
 100         *info = -2;
 101     } else if (*ldz < 1 || (icompz > 0 && *ldz < ((*n>1) ? *n : 1))) {
 102         *info = -6;
 103     }
 104     if (*info != 0) {
 105         return;
 106     }
 107
 108
 109     if (*n == 0) {
 110         return;
 111     }
 112
 113     if (*n == 1) {
 114         if (icompz == 2) {
 115             z__[z_dim1 + 1] = 1.;
 116         }
 117         return;
 118     }
 119
 120     eps = GMX_DOUBLE_EPS;
 121     d__1 = eps;
 122     eps2 = d__1 * d__1;
 123     minval = GMX_DOUBLE_MIN;
 124     safmin = minval*(1.0+GMX_DOUBLE_EPS);
 125
 126     safmax = 1. / safmin;
 127     ssfmax = sqrt(safmax) / 3.;
 128     ssfmin = sqrt(safmin) / eps2;
 129
 130     if (icompz == 2) {
 131         F77_FUNC(dlaset,DLASET)("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
 132     }
 133
 134     nmaxit = *n * 30;
 135     jtot = 0;
 136
 137     l1 = 1;
 138     nm1 = *n - 1;
 139
 140 L10:
 141     if (l1 > *n) {
 142         goto L160;
 143     }
 144     if (l1 > 1) {
 145         e[l1 - 1] = 0.;
 146     }
 147     if (l1 <= nm1) {
 148         i__1 = nm1;
 149         for (m = l1; m <= i__1; ++m) {
 150             tst = fabs(e[m]);
 151             if (fabs(tst)<GMX_DOUBLE_MIN) {
 152                 goto L30;
 153             }
 154             if (tst <= sqrt(fabs(d__[m])) * sqrt(fabs(d__[m + 1])) * eps) {
 155                 e[m] = 0.;
 156                 goto L30;
 157             }
 158         }
 159     }
 160     m = *n;
 161
 162 L30:
 163     l = l1;
 164     lsv = l;
 165     lend = m;
 166     lendsv = lend;
 167     l1 = m + 1;
 168     if (lend == l) {
 169         goto L10;
 170     }
 171
 172     i__1 = lend - l + 1;
 173     anorm = F77_FUNC(dlanst,DLANST)("I", &i__1, &d__[l], &e[l]);
 174     iscale = 0;
 175     if (fabs(anorm)<GMX_DOUBLE_MIN) {
 176         goto L10;
 177     }
 178     if (anorm > ssfmax) {
 179         iscale = 1;
 180         i__1 = lend - l + 1;
 181         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
 182                 info);
 183         i__1 = lend - l;
 184         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
 185                 info);
 186     } else if (anorm < ssfmin) {
 187         iscale = 2;
 188         i__1 = lend - l + 1;
 189         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
 190                 info);
 191         i__1 = lend - l;
 192         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
 193                 info);
 194     }
 195
 196     if (fabs(d__[lend]) < fabs(d__[l])) {
 197         lend = lsv;
 198         l = lendsv;
 199     }
 200
 201     if (lend > l) {
 202
 203 L40:
 204         if (l != lend) {
 205             lendm1 = lend - 1;
 206             i__1 = lendm1;
 207             for (m = l; m <= i__1; ++m) {
 208                 d__2 = fabs(e[m]);
 209                 tst = d__2 * d__2;
 210                 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m+ 1]) + safmin) {
 211                     goto L60;
 212                 }
 213             }
 214         }
 215
 216         m = lend;
 217
 218 L60:
 219         if (m < lend) {
 220             e[m] = 0.;
 221         }
 222         p = d__[l];
 223         if (m == l) {
 224             goto L80;
 225         }
 226
 227         if (m == l + 1) {
 228             if (icompz > 0) {
 229                 F77_FUNC(dlaev2,DLAEV2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
 230                 work[l] = c__;
 231                 work[*n - 1 + l] = s;
 232                 F77_FUNC(dlasr,DLASR)("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
 233                         z__[l * z_dim1 + 1], ldz);
 234             } else {
 235                 F77_FUNC(dlae2,DLAE2)(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
 236             }
 237             d__[l] = rt1;
 238             d__[l + 1] = rt2;
 239             e[l] = 0.;
 240             l += 2;
 241             if (l <= lend) {
 242                 goto L40;
 243             }
 244             goto L140;
 245         }
 246
 247         if (jtot == nmaxit) {
 248             goto L140;
 249         }
 250         ++jtot;
 251
 252         g = (d__[l + 1] - p) / (e[l] * 2.);
 253         r__ = F77_FUNC(dlapy2,DLAPY2)(&g, &c_b10);
 254         g = d__[m] - p + e[l] / (g + ( (g>0) ? r__ : -r__ ) );
 255
 256         s = 1.;
 257         c__ = 1.;
 258         p = 0.;
 259
 260         mm1 = m - 1;
 261         i__1 = l;
 262         for (i__ = mm1; i__ >= i__1; --i__) {
 263             f = s * e[i__];
 264             b = c__ * e[i__];
 265             F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
 266             if (i__ != m - 1) {
 267                 e[i__ + 1] = r__;
 268             }
 269             g = d__[i__ + 1] - p;
 270             r__ = (d__[i__] - g) * s + c__ * 2. * b;
 271             p = s * r__;
 272             d__[i__ + 1] = g + p;
 273             g = c__ * r__ - b;
 274
 275             if (icompz > 0) {
 276                 work[i__] = c__;
 277                 work[*n - 1 + i__] = -s;
 278             }
 279         }
 280
 281         if (icompz > 0) {
 282             mm = m - l + 1;
 283             F77_FUNC(dlasr,DLASR)("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
 284                     * z_dim1 + 1], ldz);
 285         }
 286
 287         d__[l] -= p;
 288         e[l] = g;
 289         goto L40;
 290
 291 L80:
 292         d__[l] = p;
 293
 294         ++l;
 295         if (l <= lend) {
 296             goto L40;
 297         }
 298         goto L140;
 299
 300     } else {
 301
 302 L90:
 303         if (l != lend) {
 304             lendp1 = lend + 1;
 305             i__1 = lendp1;
 306             for (m = l; m >= i__1; --m) {
 307                 d__2 = fabs(e[m - 1]);
 308                 tst = d__2 * d__2;
 309                 if (tst <= eps2 * fabs(d__[m]) * fabs(d__[m- 1]) + safmin) {
 310                     goto L110;
 311                 }
 312             }
 313         }
 314
 315         m = lend;
 316
 317 L110:
 318         if (m > lend) {
 319             e[m - 1] = 0.;
 320         }
 321         p = d__[l];
 322         if (m == l) {
 323             goto L130;
 324         }
 325         if (m == l - 1) {
 326             if (icompz > 0) {
 327                 F77_FUNC(dlaev2,DLAEV2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
 328                         ;
 329                 work[m] = c__;
 330                 work[*n - 1 + m] = s;
 331                 F77_FUNC(dlasr,DLASR)("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
 332                         z__[(l - 1) * z_dim1 + 1], ldz);
 333             } else {
 334                 F77_FUNC(dlae2,DLAE2)(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
 335             }
 336             d__[l - 1] = rt1;
 337             d__[l] = rt2;
 338             e[l - 1] = 0.;
 339             l += -2;
 340             if (l >= lend) {
 341                 goto L90;
 342             }
 343             goto L140;
 344         }
 345
 346         if (jtot == nmaxit) {
 347             goto L140;
 348         }
 349         ++jtot;
 350
 351         g = (d__[l - 1] - p) / (e[l - 1] * 2.);
 352         r__ = F77_FUNC(dlapy2,DLAPY2)(&g, &c_b10);
 353         g = d__[m] - p + e[l - 1] / (g + ( (g>0) ? r__ : -r__ ));
 354
 355         s = 1.;
 356         c__ = 1.;
 357         p = 0.;
 358
 359         lm1 = l - 1;
 360         i__1 = lm1;
 361         for (i__ = m; i__ <= i__1; ++i__) {
 362             f = s * e[i__];
 363             b = c__ * e[i__];
 364             F77_FUNC(dlartg,DLARTG)(&g, &f, &c__, &s, &r__);
 365             if (i__ != m) {
 366                 e[i__ - 1] = r__;
 367             }
 368             g = d__[i__] - p;
 369             r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
 370             p = s * r__;
 371             d__[i__] = g + p;
 372             g = c__ * r__ - b;
 373
 374             if (icompz > 0) {
 375                 work[i__] = c__;
 376                 work[*n - 1 + i__] = s;
 377             }
 378         }
 379
 380         if (icompz > 0) {
 381             mm = l - m + 1;
 382             F77_FUNC(dlasr,DLASR)("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
 383                     * z_dim1 + 1], ldz);
 384         }
 385
 386         d__[l] -= p;
 387         e[lm1] = g;
 388         goto L90;
 389
 390 L130:
 391         d__[l] = p;
 392
 393         --l;
 394         if (l >= lend) {
 395             goto L90;
 396         }
 397         goto L140;
 398
 399     }
 400
 401 L140:
 402     if (iscale == 1) {
 403         i__1 = lendsv - lsv + 1;
 404         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
 405                 n, info);
 406         i__1 = lendsv - lsv;
 407         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
 408                 info);
 409     } else if (iscale == 2) {
 410         i__1 = lendsv - lsv + 1;
 411         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
 412                 n, info);
 413         i__1 = lendsv - lsv;
 414         F77_FUNC(dlascl,DLASCL)("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
 415                 info);
 416     }
 417
 418     if (jtot < nmaxit) {
 419         goto L10;
 420     }
 421     i__1 = *n - 1;
 422     for (i__ = 1; i__ <= i__1; ++i__) {
 423         if (fabs(e[i__])>GMX_DOUBLE_MIN) {
 424             ++(*info);
 425         }
 426     }
 427     goto L190;
 428
 429 L160:
 430     if (icompz == 0) {
 431
 432         F77_FUNC(dlasrt,DLASRT)("I", n, &d__[1], info);
 433
 434     } else {
 435
 436         i__1 = *n;
 437         for (ii = 2; ii <= i__1; ++ii) {
 438             i__ = ii - 1;
 439             k = i__;
 440             p = d__[i__];
 441             i__2 = *n;
 442             for (j = ii; j <= i__2; ++j) {
 443                 if (d__[j] < p) {
 444                     k = j;
 445                     p = d__[j];
 446                 }
 447             }
 448             if (k != i__) {
 449                 d__[k] = d__[i__];
 450                 d__[i__] = p;
 451                 F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],
 452                          &c__1);
 453             }
 454         }
 455     }
 456
 457 L190:
 458     return;
 459 }
 460
 461