src/gromacs/linearalgebra/gmx_lapack/dbdsqr.c

   1 #include <ctype.h>
   2 #include <math.h>
   3
   4 #include "../gmx_blas.h"
   5 #include "../gmx_lapack.h"
   6
   7 #include "gromacs/utility/real.h"
   8
   9 void
  10 F77_FUNC(dbdsqr,DBDSQR)(const char *uplo,
  11                         int *n,
  12                         int *ncvt,
  13                         int *nru,
  14                         int *ncc,
  15                         double *d__,
  16                         double *e,
  17                         double *vt,
  18                         int *ldvt,
  19                         double *u,
  20                         int *ldu,
  21                         double *c__,
  22                         int *ldc,
  23                         double *work,
  24                         int *info)
  25 {
  26     const char xuplo = toupper(*uplo);
  27     int c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  28             i__2;
  29     double r__1, r__2, r__3, r__4;
  30     double c_b15 = -.125;
  31
  32     int c__1 = 1;
  33     double c_b49 = 1.f;
  34     double c_b72 = -1.f;
  35
  36     double f, g, h__;
  37     int i__, j, m;
  38     double r__, cs;
  39     int ll;
  40     double sn, mu;
  41     int nm1, nm12, nm13, lll;
  42     double eps, sll, tol, abse;
  43     int idir;
  44     double abss;
  45     int oldm;
  46     double cosl;
  47     int isub, iter;
  48     double unfl, sinl, cosr, smin, smax, sinr;
  49     double oldcs;
  50     int oldll;
  51     double shift, sigmn, oldsn;
  52     int maxit;
  53     double sminl;
  54     double sigmx;
  55     int lower;
  56     double sminoa;
  57     double thresh;
  58     int rotate;
  59     double sminlo, tolmul;
  60     int itmp1,itmp2;
  61     double ftmp;
  62
  63     --d__;
  64     --e;
  65     vt_dim1 = *ldvt;
  66     vt_offset = 1 + vt_dim1;
  67     vt -= vt_offset;
  68     u_dim1 = *ldu;
  69     u_offset = 1 + u_dim1;
  70     u -= u_offset;
  71     c_dim1 = *ldc;
  72     c_offset = 1 + c_dim1;
  73     c__ -= c_offset;
  74     --work;
  75
  76     *info = 0;
  77
  78     itmp1 = (*n > 1) ? *n : 1;
  79     itmp2 = (*nru > 1) ? *nru : 1;
  80
  81     lower = (xuplo == 'L');
  82     if ( (xuplo!='U') && !lower) {
  83         *info = -1;
  84     } else if (*n < 0) {
  85         *info = -2;
  86     } else if (*ncvt < 0) {
  87         *info = -3;
  88     } else if (*nru < 0) {
  89         *info = -4;
  90     } else if (*ncc < 0) {
  91         *info = -5;
  92     } else if ( ((*ncvt == 0) && (*ldvt < 1)) || ((*ncvt > 0) && (*ldvt < itmp1)) ) {
  93         *info = -9;
  94     } else if (*ldu < itmp2) {
  95         *info = -11;
  96     } else if ( ((*ncc == 0) && (*ldc < 1)) || ((*ncc > 0) && (*ldc < itmp1))) {
  97         *info = -13;
  98     }
  99     if (*info != 0) {
 100         return;
 101     }
 102     if (*n == 0) {
 103         return;
 104     }
 105     if (*n == 1) {
 106         goto L160;
 107     }
 108
 109     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
 110
 111     if (! rotate) {
 112         F77_FUNC(dlasq1,DLASQ1)(n, &d__[1], &e[1], &work[1], info);
 113         return;
 114     }
 115
 116     nm1 = *n - 1;
 117     nm12 = nm1 + nm1;
 118     nm13 = nm12 + nm1;
 119     idir = 0;
 120
 121     eps = GMX_DOUBLE_EPS;
 122     unfl = GMX_DOUBLE_MIN/GMX_DOUBLE_EPS;
 123
 124     if (lower) {
 125         i__1 = *n - 1;
 126         for (i__ = 1; i__ <= i__1; ++i__) {
 127             F77_FUNC(dlartg,DLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
 128             d__[i__] = r__;
 129             e[i__] = sn * d__[i__ + 1];
 130             d__[i__ + 1] = cs * d__[i__ + 1];
 131             work[i__] = cs;
 132             work[nm1 + i__] = sn;
 133         }
 134
 135         if (*nru > 0) {
 136             F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset],
 137                     ldu);
 138         }
 139         if (*ncc > 0) {
 140             F77_FUNC(dlasr,DLASR)("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset],
 141                      ldc);
 142         }
 143     }
 144
 145     r__3 = 100.f, r__4 = pow(GMX_DOUBLE_EPS,c_b15);
 146     r__1 = 10.f, r__2 = (r__3<r__4) ? r__3 : r__4;
 147     tolmul = (r__1>r__2) ? r__1 : r__2;
 148     tol = tolmul * eps;
 149     smax = 0.f;
 150     i__1 = *n;
 151     for (i__ = 1; i__ <= i__1; ++i__) {
 152         r__2 = smax, r__3 = (r__1 = d__[i__], fabs(r__1));
 153         smax = (r__2>r__3) ? r__2 : r__3;
 154     }
 155     i__1 = *n - 1;
 156     for (i__ = 1; i__ <= i__1; ++i__) {
 157         r__2 = smax, r__3 = (r__1 = e[i__], fabs(r__1));
 158         smax = (r__2>r__3) ? r__2 : r__3;
 159     }
 160     sminl = 0.f;
 161     if (tol >= 0.f) {
 162         sminoa = fabs(d__[1]);
 163         if (sminoa == 0.f) {
 164             goto L50;
 165         }
 166         mu = sminoa;
 167         i__1 = *n;
 168         for (i__ = 2; i__ <= i__1; ++i__) {
 169             mu = (r__2 = d__[i__], fabs(r__2)) * (mu / (mu + (r__1 = e[i__ -
 170                     1], fabs(r__1))));
 171             sminoa = (sminoa<mu) ? sminoa : mu;
 172             if (sminoa == 0.f) {
 173                 goto L50;
 174             }
 175         }
 176 L50:
 177         sminoa /= sqrt((double) (*n));
 178         r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl;
 179         thresh = (r__1>r__2) ? r__1 : r__2;
 180     } else {
 181         r__1 = fabs(tol) * smax, r__2 = *n * 6 * *n * unfl;
 182         thresh = (r__1>r__2) ? r__1 : r__2;
 183     }
 184     maxit = *n * 6 * *n;
 185     iter = 0;
 186     oldll = -1;
 187     oldm = -1;
 188     m = *n;
 189
 190 L60:
 191
 192     if (m <= 1) {
 193         goto L160;
 194     }
 195     if (iter > maxit) {
 196         goto L200;
 197     }
 198
 199     if (tol < 0.f && (r__1 = d__[m], fabs(r__1)) <= thresh) {
 200         d__[m] = 0.f;
 201     }
 202     smax = (r__1 = d__[m], fabs(r__1));
 203     smin = smax;
 204     i__1 = m - 1;
 205     for (lll = 1; lll <= i__1; ++lll) {
 206         ll = m - lll;
 207         abss = (r__1 = d__[ll], fabs(r__1));
 208         abse = (r__1 = e[ll], fabs(r__1));
 209         if (tol < 0.f && abss <= thresh) {
 210             d__[ll] = 0.f;
 211         }
 212         if (abse <= thresh) {
 213             goto L80;
 214         }
 215         smin = (smin<abss) ? smin : abss;
 216         r__1 = (smax>abss) ? smax : abss;
 217         smax = (r__1>abse) ? r__1 : abse;
 218     }
 219     ll = 0;
 220     goto L90;
 221 L80:
 222     e[ll] = 0.f;
 223     if (ll == m - 1) {
 224         --m;
 225         goto L60;
 226     }
 227 L90:
 228     ++ll;
 229     if (ll == m - 1) {
 230         F77_FUNC(dlasv2,DLASV2)(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr,
 231                  &sinl, &cosl);
 232         d__[m - 1] = sigmx;
 233         e[m - 1] = 0.f;
 234         d__[m] = sigmn;
 235         if (*ncvt > 0) {
 236             F77_FUNC(drot,DROT)(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &
 237                     cosr, &sinr);
 238         }
 239         if (*nru > 0) {
 240             F77_FUNC(drot,DROT)(nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &
 241                     c__1, &cosl, &sinl);
 242         }
 243         if (*ncc > 0) {
 244             F77_FUNC(drot,DROT)(ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, &
 245                     cosl, &sinl);
 246         }
 247         m += -2;
 248         goto L60;
 249     }
 250     if (ll > oldm || m < oldll) {
 251         if ((r__1 = d__[ll], fabs(r__1)) >= (r__2 = d__[m], fabs(r__2))) {
 252             idir = 1;
 253         } else {
 254             idir = 2;
 255         }
 256     }
 257     if (idir == 1) {
 258
 259         if( (fabs(e[m-1]) <= fabs(tol) * fabs(d__[m])) ||
 260             (tol<0.0 && fabs(e[m-1])<=thresh)) {
 261             e[m - 1] = 0.f;
 262             goto L60;
 263         }
 264         if (tol >= 0.f) {
 265             mu = (r__1 = d__[ll], fabs(r__1));
 266             sminl = mu;
 267             i__1 = m - 1;
 268             for (lll = ll; lll <= i__1; ++lll) {
 269                 if ((r__1 = e[lll], fabs(r__1)) <= tol * mu) {
 270                     e[lll] = 0.f;
 271                     goto L60;
 272                 }
 273                 sminlo = sminl;
 274                 mu = (r__2 = d__[lll + 1], fabs(r__2)) * (mu / (mu + (r__1 =
 275                         e[lll], fabs(r__1))));
 276                 sminl = (sminl<mu) ? sminl : mu;
 277             }
 278         }
 279     } else {
 280         if( (fabs(e[ll]) <= fabs(tol)*fabs(d__[ll])) ||
 281             (tol<0.0 && fabs(e[ll])<=thresh)) {
 282             e[ll] = 0.f;
 283             goto L60;
 284         }
 285         if (tol >= 0.f) {
 286             mu = (r__1 = d__[m], fabs(r__1));
 287             sminl = mu;
 288             i__1 = ll;
 289             for (lll = m - 1; lll >= i__1; --lll) {
 290                 if ((r__1 = e[lll], fabs(r__1)) <= tol * mu) {
 291                     e[lll] = 0.f;
 292                     goto L60;
 293                 }
 294                 sminlo = sminl;
 295                 mu = (r__2 = d__[lll], fabs(r__2)) * (mu / (mu + (r__1 = e[
 296                         lll], fabs(r__1))));
 297                 sminl = (sminl<mu) ? sminl : mu;
 298             }
 299         }
 300     }
 301     oldll = ll;
 302     oldm = m;
 303
 304     r__1 = eps, r__2 = tol * .01f;
 305     if (tol >= 0.f && *n * tol * (sminl / smax) <= ((r__1>r__2) ? r__1 : r__2)) {
 306         shift = 0.f;
 307     } else {
 308         if (idir == 1) {
 309             sll = (r__1 = d__[ll], fabs(r__1));
 310             F77_FUNC(dlas2,DLAS2)(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
 311         } else {
 312             sll = (r__1 = d__[m], fabs(r__1));
 313             F77_FUNC(dlas2,DLAS2)(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
 314         }
 315         if (sll > 0.f) {
 316             r__1 = shift / sll;
 317             if (r__1 * r__1 < eps) {
 318                 shift = 0.f;
 319             }
 320         }
 321     }
 322     iter = iter + m - ll;
 323     if (shift == 0.f) {
 324         if (idir == 1) {
 325             cs = 1.f;
 326             oldcs = 1.f;
 327             i__1 = m - 1;
 328             for (i__ = ll; i__ <= i__1; ++i__) {
 329                 r__1 = d__[i__] * cs;
 330                 F77_FUNC(dlartg,DLARTG)(&r__1, &e[i__], &cs, &sn, &r__);
 331                 if (i__ > ll) {
 332                     e[i__ - 1] = oldsn * r__;
 333                 }
 334                 r__1 = oldcs * r__;
 335                 r__2 = d__[i__ + 1] * sn;
 336                 F77_FUNC(dlartg,DLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
 337                 work[i__ - ll + 1] = cs;
 338                 work[i__ - ll + 1 + nm1] = sn;
 339                 work[i__ - ll + 1 + nm12] = oldcs;
 340                 work[i__ - ll + 1 + nm13] = oldsn;
 341             }
 342             h__ = d__[m] * cs;
 343             d__[m] = h__ * oldcs;
 344             e[m - 1] = h__ * oldsn;
 345             if (*ncvt > 0) {
 346                 i__1 = m - ll + 1;
 347                 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
 348                         ll + vt_dim1], ldvt);
 349             }
 350             if (*nru > 0) {
 351                 i__1 = m - ll + 1;
 352                 F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
 353                         + 1], &u[ll * u_dim1 + 1], ldu);
 354             }
 355             if (*ncc > 0) {
 356                 i__1 = m - ll + 1;
 357                 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
 358                         + 1], &c__[ll + c_dim1], ldc);
 359             }
 360             if ((r__1 = e[m - 1], fabs(r__1)) <= thresh) {
 361                 e[m - 1] = 0.f;
 362             }
 363         } else {
 364             cs = 1.f;
 365             oldcs = 1.f;
 366             i__1 = ll + 1;
 367             for (i__ = m; i__ >= i__1; --i__) {
 368                 r__1 = d__[i__] * cs;
 369                 F77_FUNC(dlartg,DLARTG)(&r__1, &e[i__ - 1], &cs, &sn, &r__);
 370                 if (i__ < m) {
 371                     e[i__] = oldsn * r__;
 372                 }
 373                 r__1 = oldcs * r__;
 374                 r__2 = d__[i__ - 1] * sn;
 375                 F77_FUNC(dlartg,DLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
 376                 work[i__ - ll] = cs;
 377                 work[i__ - ll + nm1] = -sn;
 378                 work[i__ - ll + nm12] = oldcs;
 379                 work[i__ - ll + nm13] = -oldsn;
 380             }
 381             h__ = d__[ll] * cs;
 382             d__[ll] = h__ * oldcs;
 383             e[ll] = h__ * oldsn;
 384             if (*ncvt > 0) {
 385                 i__1 = m - ll + 1;
 386                 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
 387                         nm13 + 1], &vt[ll + vt_dim1], ldvt);
 388             }
 389             if (*nru > 0) {
 390                 i__1 = m - ll + 1;
 391                 F77_FUNC(dlasr,DLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
 392                          u_dim1 + 1], ldu);
 393             }
 394             if (*ncc > 0) {
 395                 i__1 = m - ll + 1;
 396                 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
 397                         ll + c_dim1], ldc);
 398             }
 399             if ((r__1 = e[ll], fabs(r__1)) <= thresh) {
 400                 e[ll] = 0.f;
 401             }
 402         }
 403     } else {
 404
 405         if (idir == 1) {
 406             f = ((r__1 = d__[ll], fabs(r__1)) - shift) * ( ((d__[ll] > 0) ? c_b49 : -c_b49) + shift / d__[ll]);
 407             g = e[ll];
 408             i__1 = m - 1;
 409             for (i__ = ll; i__ <= i__1; ++i__) {
 410                 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosr, &sinr, &r__);
 411                 if (i__ > ll) {
 412                     e[i__ - 1] = r__;
 413                 }
 414                 f = cosr * d__[i__] + sinr * e[i__];
 415                 e[i__] = cosr * e[i__] - sinr * d__[i__];
 416                 g = sinr * d__[i__ + 1];
 417                 d__[i__ + 1] = cosr * d__[i__ + 1];
 418                 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosl, &sinl, &r__);
 419                 d__[i__] = r__;
 420                 f = cosl * e[i__] + sinl * d__[i__ + 1];
 421                 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
 422                 if (i__ < m - 1) {
 423                     g = sinl * e[i__ + 1];
 424                     e[i__ + 1] = cosl * e[i__ + 1];
 425                 }
 426                 work[i__ - ll + 1] = cosr;
 427                 work[i__ - ll + 1 + nm1] = sinr;
 428                 work[i__ - ll + 1 + nm12] = cosl;
 429                 work[i__ - ll + 1 + nm13] = sinl;
 430             }
 431             e[m - 1] = f;
 432
 433             if (*ncvt > 0) {
 434                 i__1 = m - ll + 1;
 435                 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
 436                         ll + vt_dim1], ldvt);
 437             }
 438             if (*nru > 0) {
 439                 i__1 = m - ll + 1;
 440                 F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
 441                         + 1], &u[ll * u_dim1 + 1], ldu);
 442             }
 443             if (*ncc > 0) {
 444                 i__1 = m - ll + 1;
 445                 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
 446                         + 1], &c__[ll + c_dim1], ldc);
 447             }
 448             if ((r__1 = e[m - 1], fabs(r__1)) <= thresh) {
 449                 e[m - 1] = 0.f;
 450             }
 451         } else {
 452
 453             f = ((r__1 = d__[m], fabs(r__1)) - shift) * ( ((d__[m] > 0) ? c_b49 : -c_b49) + shift / d__[m]);
 454             g = e[m - 1];
 455             i__1 = ll + 1;
 456             for (i__ = m; i__ >= i__1; --i__) {
 457                 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosr, &sinr, &r__);
 458                 if (i__ < m) {
 459                     e[i__] = r__;
 460                 }
 461                 f = cosr * d__[i__] + sinr * e[i__ - 1];
 462                 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
 463                 g = sinr * d__[i__ - 1];
 464                 d__[i__ - 1] = cosr * d__[i__ - 1];
 465                 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosl, &sinl, &r__);
 466                 d__[i__] = r__;
 467                 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
 468                 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
 469                 if (i__ > ll + 1) {
 470                     g = sinl * e[i__ - 2];
 471                     e[i__ - 2] = cosl * e[i__ - 2];
 472                 }
 473                 work[i__ - ll] = cosr;
 474                 work[i__ - ll + nm1] = -sinr;
 475                 work[i__ - ll + nm12] = cosl;
 476                 work[i__ - ll + nm13] = -sinl;
 477             }
 478             e[ll] = f;
 479
 480             if ((r__1 = e[ll], fabs(r__1)) <= thresh) {
 481                 e[ll] = 0.f;
 482             }
 483             if (*ncvt > 0) {
 484                 i__1 = m - ll + 1;
 485                 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
 486                         nm13 + 1], &vt[ll + vt_dim1], ldvt);
 487             }
 488             if (*nru > 0) {
 489                 i__1 = m - ll + 1;
 490                 F77_FUNC(dlasr,DLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
 491                          u_dim1 + 1], ldu);
 492             }
 493             if (*ncc > 0) {
 494                 i__1 = m - ll + 1;
 495                 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
 496                         ll + c_dim1], ldc);
 497             }
 498         }
 499     }
 500
 501     goto L60;
 502
 503 L160:
 504     i__1 = *n;
 505     for (i__ = 1; i__ <= i__1; ++i__) {
 506         if (d__[i__] < 0.f) {
 507             d__[i__] = -d__[i__];
 508
 509             if (*ncvt > 0) {
 510                 F77_FUNC(dscal,DSCAL)(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
 511             }
 512         }
 513     }
 514
 515     i__1 = *n - 1;
 516     for (i__ = 1; i__ <= i__1; ++i__) {
 517
 518         isub = 1;
 519         smin = d__[1];
 520         i__2 = *n + 1 - i__;
 521         for (j = 2; j <= i__2; ++j) {
 522             if (d__[j] <= smin) {
 523                 isub = j;
 524                 smin = d__[j];
 525             }
 526         }
 527         if (isub != *n + 1 - i__) {
 528             d__[isub] = d__[*n + 1 - i__];
 529             d__[*n + 1 - i__] = smin;
 530             if (*ncvt > 0) {
 531                 F77_FUNC(dswap,DSWAP)(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ +
 532                         vt_dim1], ldvt);
 533             }
 534             if (*nru > 0) {
 535                 F77_FUNC(dswap,DSWAP)(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) *
 536                         u_dim1 + 1], &c__1);
 537             }
 538             if (*ncc > 0) {
 539                 F77_FUNC(dswap,DSWAP)(ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ +
 540                         c_dim1], ldc);
 541             }
 542         }
 543     }
 544     goto L220;
 545
 546 L200:
 547     *info = 0;
 548     i__1 = *n - 1;
 549     for (i__ = 1; i__ <= i__1; ++i__) {
 550         if (e[i__] != 0.f) {
 551             ++(*info);
 552         }
 553     }
 554 L220:
 555     return;
 556
 557 }
 558
 559