git.maemo.org Git - opencv/blob - 3rdparty/lapack/dbdsqr.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static doublereal c_b15 = -.125;
   6 static integer c__1 = 1;
   7 static doublereal c_b49 = 1.;
   8 static doublereal c_b72 = -1.;
   9
  10 /* Subroutine */ int dbdsqr_(char *uplo, integer *n, integer *ncvt, integer *
  11         nru, integer *ncc, doublereal *d__, doublereal *e, doublereal *vt,
  12         integer *ldvt, doublereal *u, integer *ldu, doublereal *c__, integer *
  13         ldc, doublereal *work, integer *info)
  14 {
  15     /* System generated locals */
  16     integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  17             i__2;
  18     doublereal d__1, d__2, d__3, d__4;
  19
  20     /* Builtin functions */
  21     double pow_dd(doublereal *, doublereal *), sqrt(doublereal), d_sign(
  22             doublereal *, doublereal *);
  23
  24     /* Local variables */
  25     doublereal f, g, h__;
  26     integer i__, j, m;
  27     doublereal r__, cs;
  28     integer ll;
  29     doublereal sn, mu;
  30     integer nm1, nm12, nm13, lll;
  31     doublereal eps, sll, tol, abse;
  32     integer idir;
  33     doublereal abss;
  34     integer oldm;
  35     doublereal cosl;
  36     integer isub, iter;
  37     doublereal unfl, sinl, cosr, smin, smax, sinr;
  38     extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
  39             doublereal *, integer *, doublereal *, doublereal *), dlas2_(
  40             doublereal *, doublereal *, doublereal *, doublereal *,
  41             doublereal *), dscal_(integer *, doublereal *, doublereal *,
  42             integer *);
  43     extern logical lsame_(char *, char *);
  44     doublereal oldcs;
  45     extern /* Subroutine */ int dlasr_(char *, char *, char *, integer *,
  46             integer *, doublereal *, doublereal *, doublereal *, integer *);
  47     integer oldll;
  48     doublereal shift, sigmn, oldsn;
  49     extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
  50             doublereal *, integer *);
  51     integer maxit;
  52     doublereal sminl, sigmx;
  53     logical lower;
  54     extern /* Subroutine */ int dlasq1_(integer *, doublereal *, doublereal *,
  55              doublereal *, integer *), dlasv2_(doublereal *, doublereal *,
  56             doublereal *, doublereal *, doublereal *, doublereal *,
  57             doublereal *, doublereal *, doublereal *);
  58     extern doublereal dlamch_(char *);
  59     extern /* Subroutine */ int dlartg_(doublereal *, doublereal *,
  60             doublereal *, doublereal *, doublereal *), xerbla_(char *,
  61             integer *);
  62     doublereal sminoa, thresh;
  63     logical rotate;
  64     doublereal tolmul;
  65
  66
  67 /*  -- LAPACK routine (version 3.1.1) -- */
  68 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  69 /*     January 2007 */
  70
  71 /*     .. Scalar Arguments .. */
  72 /*     .. */
  73 /*     .. Array Arguments .. */
  74 /*     .. */
  75
  76 /*  Purpose */
  77 /*  ======= */
  78
  79 /*  DBDSQR computes the singular values and, optionally, the right and/or */
  80 /*  left singular vectors from the singular value decomposition (SVD) of */
  81 /*  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit */
  82 /*  zero-shift QR algorithm.  The SVD of B has the form */
  83
  84 /*     B = Q * S * P**T */
  85
  86 /*  where S is the diagonal matrix of singular values, Q is an orthogonal */
  87 /*  matrix of left singular vectors, and P is an orthogonal matrix of */
  88 /*  right singular vectors.  If left singular vectors are requested, this */
  89 /*  subroutine actually returns U*Q instead of Q, and, if right singular */
  90 /*  vectors are requested, this subroutine returns P**T*VT instead of */
  91 /*  P**T, for given real input matrices U and VT.  When U and VT are the */
  92 /*  orthogonal matrices that reduce a general matrix A to bidiagonal */
  93 /*  form:  A = U*B*VT, as computed by DGEBRD, then */
  94
  95 /*     A = (U*Q) * S * (P**T*VT) */
  96
  97 /*  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C */
  98 /*  for a given real input matrix C. */
  99
 100 /*  See "Computing  Small Singular Values of Bidiagonal Matrices With */
 101 /*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
 102 /*  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, */
 103 /*  no. 5, pp. 873-912, Sept 1990) and */
 104 /*  "Accurate singular values and differential qd algorithms," by */
 105 /*  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics */
 106 /*  Department, University of California at Berkeley, July 1992 */
 107 /*  for a detailed description of the algorithm. */
 108
 109 /*  Arguments */
 110 /*  ========= */
 111
 112 /*  UPLO    (input) CHARACTER*1 */
 113 /*          = 'U':  B is upper bidiagonal; */
 114 /*          = 'L':  B is lower bidiagonal. */
 115
 116 /*  N       (input) INTEGER */
 117 /*          The order of the matrix B.  N >= 0. */
 118
 119 /*  NCVT    (input) INTEGER */
 120 /*          The number of columns of the matrix VT. NCVT >= 0. */
 121
 122 /*  NRU     (input) INTEGER */
 123 /*          The number of rows of the matrix U. NRU >= 0. */
 124
 125 /*  NCC     (input) INTEGER */
 126 /*          The number of columns of the matrix C. NCC >= 0. */
 127
 128 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
 129 /*          On entry, the n diagonal elements of the bidiagonal matrix B. */
 130 /*          On exit, if INFO=0, the singular values of B in decreasing */
 131 /*          order. */
 132
 133 /*  E       (input/output) DOUBLE PRECISION array, dimension (N-1) */
 134 /*          On entry, the N-1 offdiagonal elements of the bidiagonal */
 135 /*          matrix B. */
 136 /*          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E */
 137 /*          will contain the diagonal and superdiagonal elements of a */
 138 /*          bidiagonal matrix orthogonally equivalent to the one given */
 139 /*          as input. */
 140
 141 /*  VT      (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) */
 142 /*          On entry, an N-by-NCVT matrix VT. */
 143 /*          On exit, VT is overwritten by P**T * VT. */
 144 /*          Not referenced if NCVT = 0. */
 145
 146 /*  LDVT    (input) INTEGER */
 147 /*          The leading dimension of the array VT. */
 148 /*          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. */
 149
 150 /*  U       (input/output) DOUBLE PRECISION array, dimension (LDU, N) */
 151 /*          On entry, an NRU-by-N matrix U. */
 152 /*          On exit, U is overwritten by U * Q. */
 153 /*          Not referenced if NRU = 0. */
 154
 155 /*  LDU     (input) INTEGER */
 156 /*          The leading dimension of the array U.  LDU >= max(1,NRU). */
 157
 158 /*  C       (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) */
 159 /*          On entry, an N-by-NCC matrix C. */
 160 /*          On exit, C is overwritten by Q**T * C. */
 161 /*          Not referenced if NCC = 0. */
 162
 163 /*  LDC     (input) INTEGER */
 164 /*          The leading dimension of the array C. */
 165 /*          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. */
 166
 167 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N) */
 168 /*          if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise */
 169
 170 /*  INFO    (output) INTEGER */
 171 /*          = 0:  successful exit */
 172 /*          < 0:  If INFO = -i, the i-th argument had an illegal value */
 173 /*          > 0:  the algorithm did not converge; D and E contain the */
 174 /*                elements of a bidiagonal matrix which is orthogonally */
 175 /*                similar to the input matrix B;  if INFO = i, i */
 176 /*                elements of E have not converged to zero. */
 177
 178 /*  Internal Parameters */
 179 /*  =================== */
 180
 181 /*  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) */
 182 /*          TOLMUL controls the convergence criterion of the QR loop. */
 183 /*          If it is positive, TOLMUL*EPS is the desired relative */
 184 /*             precision in the computed singular values. */
 185 /*          If it is negative, abs(TOLMUL*EPS*sigma_max) is the */
 186 /*             desired absolute accuracy in the computed singular */
 187 /*             values (corresponds to relative accuracy */
 188 /*             abs(TOLMUL*EPS) in the largest singular value. */
 189 /*          abs(TOLMUL) should be between 1 and 1/EPS, and preferably */
 190 /*             between 10 (for fast convergence) and .1/EPS */
 191 /*             (for there to be some accuracy in the results). */
 192 /*          Default is to lose at either one eighth or 2 of the */
 193 /*             available decimal digits in each computed singular value */
 194 /*             (whichever is smaller). */
 195
 196 /*  MAXITR  INTEGER, default = 6 */
 197 /*          MAXITR controls the maximum number of passes of the */
 198 /*          algorithm through its inner loop. The algorithms stops */
 199 /*          (and so fails to converge) if the number of passes */
 200 /*          through the inner loop exceeds MAXITR*N**2. */
 201
 202 /*  ===================================================================== */
 203
 204 /*     .. Parameters .. */
 205 /*     .. */
 206 /*     .. Local Scalars .. */
 207 /*     .. */
 208 /*     .. External Functions .. */
 209 /*     .. */
 210 /*     .. External Subroutines .. */
 211 /*     .. */
 212 /*     .. Intrinsic Functions .. */
 213 /*     .. */
 214 /*     .. Executable Statements .. */
 215
 216 /*     Test the input parameters. */
 217
 218     /* Parameter adjustments */
 219     --d__;
 220     --e;
 221     vt_dim1 = *ldvt;
 222     vt_offset = 1 + vt_dim1;
 223     vt -= vt_offset;
 224     u_dim1 = *ldu;
 225     u_offset = 1 + u_dim1;
 226     u -= u_offset;
 227     c_dim1 = *ldc;
 228     c_offset = 1 + c_dim1;
 229     c__ -= c_offset;
 230     --work;
 231
 232     /* Function Body */
 233     *info = 0;
 234     lower = lsame_(uplo, "L");
 235     if (! lsame_(uplo, "U") && ! lower) {
 236         *info = -1;
 237     } else if (*n < 0) {
 238         *info = -2;
 239     } else if (*ncvt < 0) {
 240         *info = -3;
 241     } else if (*nru < 0) {
 242         *info = -4;
 243     } else if (*ncc < 0) {
 244         *info = -5;
 245     } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) {
 246         *info = -9;
 247     } else if (*ldu < max(1,*nru)) {
 248         *info = -11;
 249     } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) {
 250         *info = -13;
 251     }
 252     if (*info != 0) {
 253         i__1 = -(*info);
 254         xerbla_("DBDSQR", &i__1);
 255         return 0;
 256     }
 257     if (*n == 0) {
 258         return 0;
 259     }
 260     if (*n == 1) {
 261         goto L160;
 262     }
 263
 264 /*     ROTATE is true if any singular vectors desired, false otherwise */
 265
 266     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
 267
 268 /*     If no singular vectors desired, use qd algorithm */
 269
 270     if (! rotate) {
 271         dlasq1_(n, &d__[1], &e[1], &work[1], info);
 272         return 0;
 273     }
 274
 275     nm1 = *n - 1;
 276     nm12 = nm1 + nm1;
 277     nm13 = nm12 + nm1;
 278     idir = 0;
 279
 280 /*     Get machine constants */
 281
 282     eps = dlamch_("Epsilon");
 283     unfl = dlamch_("Safe minimum");
 284
 285 /*     If matrix lower bidiagonal, rotate to be upper bidiagonal */
 286 /*     by applying Givens rotations on the left */
 287
 288     if (lower) {
 289         i__1 = *n - 1;
 290         for (i__ = 1; i__ <= i__1; ++i__) {
 291             dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 292             d__[i__] = r__;
 293             e[i__] = sn * d__[i__ + 1];
 294             d__[i__ + 1] = cs * d__[i__ + 1];
 295             work[i__] = cs;
 296             work[nm1 + i__] = sn;
 297 /* L10: */
 298         }
 299
 300 /*        Update singular vectors if desired */
 301
 302         if (*nru > 0) {
 303             dlasr_("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset],
 304                     ldu);
 305         }
 306         if (*ncc > 0) {
 307             dlasr_("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset],
 308                      ldc);
 309         }
 310     }
 311
 312 /*     Compute singular values to relative accuracy TOL */
 313 /*     (By setting TOL to be negative, algorithm will compute */
 314 /*     singular values to absolute accuracy ABS(TOL)*norm(input matrix)) */
 315
 316 /* Computing MAX */
 317 /* Computing MIN */
 318     d__3 = 100., d__4 = pow_dd(&eps, &c_b15);
 319     d__1 = 10., d__2 = min(d__3,d__4);
 320     tolmul = max(d__1,d__2);
 321     tol = tolmul * eps;
 322
 323 /*     Compute approximate maximum, minimum singular values */
 324
 325     smax = 0.;
 326     i__1 = *n;
 327     for (i__ = 1; i__ <= i__1; ++i__) {
 328 /* Computing MAX */
 329         d__2 = smax, d__3 = (d__1 = d__[i__], abs(d__1));
 330         smax = max(d__2,d__3);
 331 /* L20: */
 332     }
 333     i__1 = *n - 1;
 334     for (i__ = 1; i__ <= i__1; ++i__) {
 335 /* Computing MAX */
 336         d__2 = smax, d__3 = (d__1 = e[i__], abs(d__1));
 337         smax = max(d__2,d__3);
 338 /* L30: */
 339     }
 340     sminl = 0.;
 341     if (tol >= 0.) {
 342
 343 /*        Relative accuracy desired */
 344
 345         sminoa = abs(d__[1]);
 346         if (sminoa == 0.) {
 347             goto L50;
 348         }
 349         mu = sminoa;
 350         i__1 = *n;
 351         for (i__ = 2; i__ <= i__1; ++i__) {
 352             mu = (d__2 = d__[i__], abs(d__2)) * (mu / (mu + (d__1 = e[i__ - 1]
 353                     , abs(d__1))));
 354             sminoa = min(sminoa,mu);
 355             if (sminoa == 0.) {
 356                 goto L50;
 357             }
 358 /* L40: */
 359         }
 360 L50:
 361         sminoa /= sqrt((doublereal) (*n));
 362 /* Computing MAX */
 363         d__1 = tol * sminoa, d__2 = *n * 6 * *n * unfl;
 364         thresh = max(d__1,d__2);
 365     } else {
 366
 367 /*        Absolute accuracy desired */
 368
 369 /* Computing MAX */
 370         d__1 = abs(tol) * smax, d__2 = *n * 6 * *n * unfl;
 371         thresh = max(d__1,d__2);
 372     }
 373
 374 /*     Prepare for main iteration loop for the singular values */
 375 /*     (MAXIT is the maximum number of passes through the inner */
 376 /*     loop permitted before nonconvergence signalled.) */
 377
 378     maxit = *n * 6 * *n;
 379     iter = 0;
 380     oldll = -1;
 381     oldm = -1;
 382
 383 /*     M points to last element of unconverged part of matrix */
 384
 385     m = *n;
 386
 387 /*     Begin main iteration loop */
 388
 389 L60:
 390
 391 /*     Check for convergence or exceeding iteration count */
 392
 393     if (m <= 1) {
 394         goto L160;
 395     }
 396     if (iter > maxit) {
 397         goto L200;
 398     }
 399
 400 /*     Find diagonal block of matrix to work on */
 401
 402     if (tol < 0. && (d__1 = d__[m], abs(d__1)) <= thresh) {
 403         d__[m] = 0.;
 404     }
 405     smax = (d__1 = d__[m], abs(d__1));
 406     smin = smax;
 407     i__1 = m - 1;
 408     for (lll = 1; lll <= i__1; ++lll) {
 409         ll = m - lll;
 410         abss = (d__1 = d__[ll], abs(d__1));
 411         abse = (d__1 = e[ll], abs(d__1));
 412         if (tol < 0. && abss <= thresh) {
 413             d__[ll] = 0.;
 414         }
 415         if (abse <= thresh) {
 416             goto L80;
 417         }
 418         smin = min(smin,abss);
 419 /* Computing MAX */
 420         d__1 = max(smax,abss);
 421         smax = max(d__1,abse);
 422 /* L70: */
 423     }
 424     ll = 0;
 425     goto L90;
 426 L80:
 427     e[ll] = 0.;
 428
 429 /*     Matrix splits since E(LL) = 0 */
 430
 431     if (ll == m - 1) {
 432
 433 /*        Convergence of bottom singular value, return to top of loop */
 434
 435         --m;
 436         goto L60;
 437     }
 438 L90:
 439     ++ll;
 440
 441 /*     E(LL) through E(M-1) are nonzero, E(LL-1) is zero */
 442
 443     if (ll == m - 1) {
 444
 445 /*        2 by 2 block, handle separately */
 446
 447         dlasv2_(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr,
 448                  &sinl, &cosl);
 449         d__[m - 1] = sigmx;
 450         e[m - 1] = 0.;
 451         d__[m] = sigmn;
 452
 453 /*        Compute singular vectors, if desired */
 454
 455         if (*ncvt > 0) {
 456             drot_(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &
 457                     cosr, &sinr);
 458         }
 459         if (*nru > 0) {
 460             drot_(nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &
 461                     c__1, &cosl, &sinl);
 462         }
 463         if (*ncc > 0) {
 464             drot_(ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, &
 465                     cosl, &sinl);
 466         }
 467         m += -2;
 468         goto L60;
 469     }
 470
 471 /*     If working on new submatrix, choose shift direction */
 472 /*     (from larger end diagonal element towards smaller) */
 473
 474     if (ll > oldm || m < oldll) {
 475         if ((d__1 = d__[ll], abs(d__1)) >= (d__2 = d__[m], abs(d__2))) {
 476
 477 /*           Chase bulge from top (big end) to bottom (small end) */
 478
 479             idir = 1;
 480         } else {
 481
 482 /*           Chase bulge from bottom (big end) to top (small end) */
 483
 484             idir = 2;
 485         }
 486     }
 487
 488 /*     Apply convergence tests */
 489
 490     if (idir == 1) {
 491
 492 /*        Run convergence test in forward direction */
 493 /*        First apply standard test to bottom of matrix */
 494
 495         if ((d__2 = e[m - 1], abs(d__2)) <= abs(tol) * (d__1 = d__[m], abs(
 496                 d__1)) || tol < 0. && (d__3 = e[m - 1], abs(d__3)) <= thresh)
 497                 {
 498             e[m - 1] = 0.;
 499             goto L60;
 500         }
 501
 502         if (tol >= 0.) {
 503
 504 /*           If relative accuracy desired, */
 505 /*           apply convergence criterion forward */
 506
 507             mu = (d__1 = d__[ll], abs(d__1));
 508             sminl = mu;
 509             i__1 = m - 1;
 510             for (lll = ll; lll <= i__1; ++lll) {
 511                 if ((d__1 = e[lll], abs(d__1)) <= tol * mu) {
 512                     e[lll] = 0.;
 513                     goto L60;
 514                 }
 515                 mu = (d__2 = d__[lll + 1], abs(d__2)) * (mu / (mu + (d__1 = e[
 516                         lll], abs(d__1))));
 517                 sminl = min(sminl,mu);
 518 /* L100: */
 519             }
 520         }
 521
 522     } else {
 523
 524 /*        Run convergence test in backward direction */
 525 /*        First apply standard test to top of matrix */
 526
 527         if ((d__2 = e[ll], abs(d__2)) <= abs(tol) * (d__1 = d__[ll], abs(d__1)
 528                 ) || tol < 0. && (d__3 = e[ll], abs(d__3)) <= thresh) {
 529             e[ll] = 0.;
 530             goto L60;
 531         }
 532
 533         if (tol >= 0.) {
 534
 535 /*           If relative accuracy desired, */
 536 /*           apply convergence criterion backward */
 537
 538             mu = (d__1 = d__[m], abs(d__1));
 539             sminl = mu;
 540             i__1 = ll;
 541             for (lll = m - 1; lll >= i__1; --lll) {
 542                 if ((d__1 = e[lll], abs(d__1)) <= tol * mu) {
 543                     e[lll] = 0.;
 544                     goto L60;
 545                 }
 546                 mu = (d__2 = d__[lll], abs(d__2)) * (mu / (mu + (d__1 = e[lll]
 547                         , abs(d__1))));
 548                 sminl = min(sminl,mu);
 549 /* L110: */
 550             }
 551         }
 552     }
 553     oldll = ll;
 554     oldm = m;
 555
 556 /*     Compute shift.  First, test if shifting would ruin relative */
 557 /*     accuracy, and if so set the shift to zero. */
 558
 559 /* Computing MAX */
 560     d__1 = eps, d__2 = tol * .01;
 561     if (tol >= 0. && *n * tol * (sminl / smax) <= max(d__1,d__2)) {
 562
 563 /*        Use a zero shift to avoid loss of relative accuracy */
 564
 565         shift = 0.;
 566     } else {
 567
 568 /*        Compute the shift from 2-by-2 block at end of matrix */
 569
 570         if (idir == 1) {
 571             sll = (d__1 = d__[ll], abs(d__1));
 572             dlas2_(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
 573         } else {
 574             sll = (d__1 = d__[m], abs(d__1));
 575             dlas2_(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
 576         }
 577
 578 /*        Test if shift negligible, and if so set to zero */
 579
 580         if (sll > 0.) {
 581 /* Computing 2nd power */
 582             d__1 = shift / sll;
 583             if (d__1 * d__1 < eps) {
 584                 shift = 0.;
 585             }
 586         }
 587     }
 588
 589 /*     Increment iteration count */
 590
 591     iter = iter + m - ll;
 592
 593 /*     If SHIFT = 0, do simplified QR iteration */
 594
 595     if (shift == 0.) {
 596         if (idir == 1) {
 597
 598 /*           Chase bulge from top to bottom */
 599 /*           Save cosines and sines for later singular vector updates */
 600
 601             cs = 1.;
 602             oldcs = 1.;
 603             i__1 = m - 1;
 604             for (i__ = ll; i__ <= i__1; ++i__) {
 605                 d__1 = d__[i__] * cs;
 606                 dlartg_(&d__1, &e[i__], &cs, &sn, &r__);
 607                 if (i__ > ll) {
 608                     e[i__ - 1] = oldsn * r__;
 609                 }
 610                 d__1 = oldcs * r__;
 611                 d__2 = d__[i__ + 1] * sn;
 612                 dlartg_(&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
 613                 work[i__ - ll + 1] = cs;
 614                 work[i__ - ll + 1 + nm1] = sn;
 615                 work[i__ - ll + 1 + nm12] = oldcs;
 616                 work[i__ - ll + 1 + nm13] = oldsn;
 617 /* L120: */
 618             }
 619             h__ = d__[m] * cs;
 620             d__[m] = h__ * oldcs;
 621             e[m - 1] = h__ * oldsn;
 622
 623 /*           Update singular vectors */
 624
 625             if (*ncvt > 0) {
 626                 i__1 = m - ll + 1;
 627                 dlasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
 628                         ll + vt_dim1], ldvt);
 629             }
 630             if (*nru > 0) {
 631                 i__1 = m - ll + 1;
 632                 dlasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
 633                         + 1], &u[ll * u_dim1 + 1], ldu);
 634             }
 635             if (*ncc > 0) {
 636                 i__1 = m - ll + 1;
 637                 dlasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
 638                         + 1], &c__[ll + c_dim1], ldc);
 639             }
 640
 641 /*           Test convergence */
 642
 643             if ((d__1 = e[m - 1], abs(d__1)) <= thresh) {
 644                 e[m - 1] = 0.;
 645             }
 646
 647         } else {
 648
 649 /*           Chase bulge from bottom to top */
 650 /*           Save cosines and sines for later singular vector updates */
 651
 652             cs = 1.;
 653             oldcs = 1.;
 654             i__1 = ll + 1;
 655             for (i__ = m; i__ >= i__1; --i__) {
 656                 d__1 = d__[i__] * cs;
 657                 dlartg_(&d__1, &e[i__ - 1], &cs, &sn, &r__);
 658                 if (i__ < m) {
 659                     e[i__] = oldsn * r__;
 660                 }
 661                 d__1 = oldcs * r__;
 662                 d__2 = d__[i__ - 1] * sn;
 663                 dlartg_(&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
 664                 work[i__ - ll] = cs;
 665                 work[i__ - ll + nm1] = -sn;
 666                 work[i__ - ll + nm12] = oldcs;
 667                 work[i__ - ll + nm13] = -oldsn;
 668 /* L130: */
 669             }
 670             h__ = d__[ll] * cs;
 671             d__[ll] = h__ * oldcs;
 672             e[ll] = h__ * oldsn;
 673
 674 /*           Update singular vectors */
 675
 676             if (*ncvt > 0) {
 677                 i__1 = m - ll + 1;
 678                 dlasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
 679                         nm13 + 1], &vt[ll + vt_dim1], ldvt);
 680             }
 681             if (*nru > 0) {
 682                 i__1 = m - ll + 1;
 683                 dlasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
 684                          u_dim1 + 1], ldu);
 685             }
 686             if (*ncc > 0) {
 687                 i__1 = m - ll + 1;
 688                 dlasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
 689                         ll + c_dim1], ldc);
 690             }
 691
 692 /*           Test convergence */
 693
 694             if ((d__1 = e[ll], abs(d__1)) <= thresh) {
 695                 e[ll] = 0.;
 696             }
 697         }
 698     } else {
 699
 700 /*        Use nonzero shift */
 701
 702         if (idir == 1) {
 703
 704 /*           Chase bulge from top to bottom */
 705 /*           Save cosines and sines for later singular vector updates */
 706
 707             f = ((d__1 = d__[ll], abs(d__1)) - shift) * (d_sign(&c_b49, &d__[
 708                     ll]) + shift / d__[ll]);
 709             g = e[ll];
 710             i__1 = m - 1;
 711             for (i__ = ll; i__ <= i__1; ++i__) {
 712                 dlartg_(&f, &g, &cosr, &sinr, &r__);
 713                 if (i__ > ll) {
 714                     e[i__ - 1] = r__;
 715                 }
 716                 f = cosr * d__[i__] + sinr * e[i__];
 717                 e[i__] = cosr * e[i__] - sinr * d__[i__];
 718                 g = sinr * d__[i__ + 1];
 719                 d__[i__ + 1] = cosr * d__[i__ + 1];
 720                 dlartg_(&f, &g, &cosl, &sinl, &r__);
 721                 d__[i__] = r__;
 722                 f = cosl * e[i__] + sinl * d__[i__ + 1];
 723                 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
 724                 if (i__ < m - 1) {
 725                     g = sinl * e[i__ + 1];
 726                     e[i__ + 1] = cosl * e[i__ + 1];
 727                 }
 728                 work[i__ - ll + 1] = cosr;
 729                 work[i__ - ll + 1 + nm1] = sinr;
 730                 work[i__ - ll + 1 + nm12] = cosl;
 731                 work[i__ - ll + 1 + nm13] = sinl;
 732 /* L140: */
 733             }
 734             e[m - 1] = f;
 735
 736 /*           Update singular vectors */
 737
 738             if (*ncvt > 0) {
 739                 i__1 = m - ll + 1;
 740                 dlasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
 741                         ll + vt_dim1], ldvt);
 742             }
 743             if (*nru > 0) {
 744                 i__1 = m - ll + 1;
 745                 dlasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
 746                         + 1], &u[ll * u_dim1 + 1], ldu);
 747             }
 748             if (*ncc > 0) {
 749                 i__1 = m - ll + 1;
 750                 dlasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
 751                         + 1], &c__[ll + c_dim1], ldc);
 752             }
 753
 754 /*           Test convergence */
 755
 756             if ((d__1 = e[m - 1], abs(d__1)) <= thresh) {
 757                 e[m - 1] = 0.;
 758             }
 759
 760         } else {
 761
 762 /*           Chase bulge from bottom to top */
 763 /*           Save cosines and sines for later singular vector updates */
 764
 765             f = ((d__1 = d__[m], abs(d__1)) - shift) * (d_sign(&c_b49, &d__[m]
 766                     ) + shift / d__[m]);
 767             g = e[m - 1];
 768             i__1 = ll + 1;
 769             for (i__ = m; i__ >= i__1; --i__) {
 770                 dlartg_(&f, &g, &cosr, &sinr, &r__);
 771                 if (i__ < m) {
 772                     e[i__] = r__;
 773                 }
 774                 f = cosr * d__[i__] + sinr * e[i__ - 1];
 775                 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
 776                 g = sinr * d__[i__ - 1];
 777                 d__[i__ - 1] = cosr * d__[i__ - 1];
 778                 dlartg_(&f, &g, &cosl, &sinl, &r__);
 779                 d__[i__] = r__;
 780                 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
 781                 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
 782                 if (i__ > ll + 1) {
 783                     g = sinl * e[i__ - 2];
 784                     e[i__ - 2] = cosl * e[i__ - 2];
 785                 }
 786                 work[i__ - ll] = cosr;
 787                 work[i__ - ll + nm1] = -sinr;
 788                 work[i__ - ll + nm12] = cosl;
 789                 work[i__ - ll + nm13] = -sinl;
 790 /* L150: */
 791             }
 792             e[ll] = f;
 793
 794 /*           Test convergence */
 795
 796             if ((d__1 = e[ll], abs(d__1)) <= thresh) {
 797                 e[ll] = 0.;
 798             }
 799
 800 /*           Update singular vectors if desired */
 801
 802             if (*ncvt > 0) {
 803                 i__1 = m - ll + 1;
 804                 dlasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
 805                         nm13 + 1], &vt[ll + vt_dim1], ldvt);
 806             }
 807             if (*nru > 0) {
 808                 i__1 = m - ll + 1;
 809                 dlasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
 810                          u_dim1 + 1], ldu);
 811             }
 812             if (*ncc > 0) {
 813                 i__1 = m - ll + 1;
 814                 dlasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
 815                         ll + c_dim1], ldc);
 816             }
 817         }
 818     }
 819
 820 /*     QR iteration finished, go back and check convergence */
 821
 822     goto L60;
 823
 824 /*     All singular values converged, so make them positive */
 825
 826 L160:
 827     i__1 = *n;
 828     for (i__ = 1; i__ <= i__1; ++i__) {
 829         if (d__[i__] < 0.) {
 830             d__[i__] = -d__[i__];
 831
 832 /*           Change sign of singular vectors, if desired */
 833
 834             if (*ncvt > 0) {
 835                 dscal_(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
 836             }
 837         }
 838 /* L170: */
 839     }
 840
 841 /*     Sort the singular values into decreasing order (insertion sort on */
 842 /*     singular values, but only one transposition per singular vector) */
 843
 844     i__1 = *n - 1;
 845     for (i__ = 1; i__ <= i__1; ++i__) {
 846
 847 /*        Scan for smallest D(I) */
 848
 849         isub = 1;
 850         smin = d__[1];
 851         i__2 = *n + 1 - i__;
 852         for (j = 2; j <= i__2; ++j) {
 853             if (d__[j] <= smin) {
 854                 isub = j;
 855                 smin = d__[j];
 856             }
 857 /* L180: */
 858         }
 859         if (isub != *n + 1 - i__) {
 860
 861 /*           Swap singular values and vectors */
 862
 863             d__[isub] = d__[*n + 1 - i__];
 864             d__[*n + 1 - i__] = smin;
 865             if (*ncvt > 0) {
 866                 dswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ +
 867                         vt_dim1], ldvt);
 868             }
 869             if (*nru > 0) {
 870                 dswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) *
 871                         u_dim1 + 1], &c__1);
 872             }
 873             if (*ncc > 0) {
 874                 dswap_(ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ +
 875                         c_dim1], ldc);
 876             }
 877         }
 878 /* L190: */
 879     }
 880     goto L220;
 881
 882 /*     Maximum number of iterations exceeded, failure to converge */
 883
 884 L200:
 885     *info = 0;
 886     i__1 = *n - 1;
 887     for (i__ = 1; i__ <= i__1; ++i__) {
 888         if (e[i__] != 0.) {
 889             ++(*info);
 890         }
 891 /* L210: */
 892     }
 893 L220:
 894     return 0;
 895
 896 /*     End of DBDSQR */
 897
 898 } /* dbdsqr_ */