git.maemo.org Git - opencv/blob - 3rdparty/lapack/sbdsqr.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 real c_b49 = 1.f;
   8 static real c_b72 = -1.f;
   9
  10 /* Subroutine */ int sbdsqr_(char *uplo, integer *n, integer *ncvt, integer *
  11         nru, integer *ncc, real *d__, real *e, real *vt, integer *ldvt, real *
  12         u, integer *ldu, real *c__, integer *ldc, real *work, integer *info)
  13 {
  14     /* System generated locals */
  15     integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
  16             i__2;
  17     real r__1, r__2, r__3, r__4;
  18     doublereal d__1;
  19
  20     /* Builtin functions */
  21     double pow_dd(doublereal *, doublereal *), sqrt(doublereal), r_sign(real *
  22             , real *);
  23
  24     /* Local variables */
  25     real f, g, h__;
  26     integer i__, j, m;
  27     real r__, cs;
  28     integer ll;
  29     real sn, mu;
  30     integer nm1, nm12, nm13, lll;
  31     real eps, sll, tol, abse;
  32     integer idir;
  33     real abss;
  34     integer oldm;
  35     real cosl;
  36     integer isub, iter;
  37     real unfl, sinl, cosr, smin, smax, sinr;
  38     extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
  39             integer *, real *, real *), slas2_(real *, real *, real *, real *,
  40              real *);
  41     extern logical lsame_(char *, char *);
  42     real oldcs;
  43     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
  44     integer oldll;
  45     real shift, sigmn, oldsn;
  46     integer maxit;
  47     real sminl;
  48     extern /* Subroutine */ int slasr_(char *, char *, char *, integer *,
  49             integer *, real *, real *, real *, integer *);
  50     real sigmx;
  51     logical lower;
  52     extern /* Subroutine */ int sswap_(integer *, real *, integer *, real *,
  53             integer *), slasq1_(integer *, real *, real *, real *, integer *),
  54              slasv2_(real *, real *, real *, real *, real *, real *, real *,
  55             real *, real *);
  56     extern doublereal slamch_(char *);
  57     extern /* Subroutine */ int xerbla_(char *, integer *);
  58     real sminoa;
  59     extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real *
  60 );
  61     real thresh;
  62     logical rotate;
  63     real tolmul;
  64
  65
  66 /*  -- LAPACK routine (version 3.1.1) -- */
  67 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  68 /*     January 2007 */
  69
  70 /*     .. Scalar Arguments .. */
  71 /*     .. */
  72 /*     .. Array Arguments .. */
  73 /*     .. */
  74
  75 /*  Purpose */
  76 /*  ======= */
  77
  78 /*  SBDSQR computes the singular values and, optionally, the right and/or */
  79 /*  left singular vectors from the singular value decomposition (SVD) of */
  80 /*  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit */
  81 /*  zero-shift QR algorithm.  The SVD of B has the form */
  82
  83 /*     B = Q * S * P**T */
  84
  85 /*  where S is the diagonal matrix of singular values, Q is an orthogonal */
  86 /*  matrix of left singular vectors, and P is an orthogonal matrix of */
  87 /*  right singular vectors.  If left singular vectors are requested, this */
  88 /*  subroutine actually returns U*Q instead of Q, and, if right singular */
  89 /*  vectors are requested, this subroutine returns P**T*VT instead of */
  90 /*  P**T, for given real input matrices U and VT.  When U and VT are the */
  91 /*  orthogonal matrices that reduce a general matrix A to bidiagonal */
  92 /*  form:  A = U*B*VT, as computed by SGEBRD, then */
  93
  94 /*     A = (U*Q) * S * (P**T*VT) */
  95
  96 /*  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C */
  97 /*  for a given real input matrix C. */
  98
  99 /*  See "Computing  Small Singular Values of Bidiagonal Matrices With */
 100 /*  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
 101 /*  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, */
 102 /*  no. 5, pp. 873-912, Sept 1990) and */
 103 /*  "Accurate singular values and differential qd algorithms," by */
 104 /*  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics */
 105 /*  Department, University of California at Berkeley, July 1992 */
 106 /*  for a detailed description of the algorithm. */
 107
 108 /*  Arguments */
 109 /*  ========= */
 110
 111 /*  UPLO    (input) CHARACTER*1 */
 112 /*          = 'U':  B is upper bidiagonal; */
 113 /*          = 'L':  B is lower bidiagonal. */
 114
 115 /*  N       (input) INTEGER */
 116 /*          The order of the matrix B.  N >= 0. */
 117
 118 /*  NCVT    (input) INTEGER */
 119 /*          The number of columns of the matrix VT. NCVT >= 0. */
 120
 121 /*  NRU     (input) INTEGER */
 122 /*          The number of rows of the matrix U. NRU >= 0. */
 123
 124 /*  NCC     (input) INTEGER */
 125 /*          The number of columns of the matrix C. NCC >= 0. */
 126
 127 /*  D       (input/output) REAL array, dimension (N) */
 128 /*          On entry, the n diagonal elements of the bidiagonal matrix B. */
 129 /*          On exit, if INFO=0, the singular values of B in decreasing */
 130 /*          order. */
 131
 132 /*  E       (input/output) REAL array, dimension (N-1) */
 133 /*          On entry, the N-1 offdiagonal elements of the bidiagonal */
 134 /*          matrix B. */
 135 /*          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E */
 136 /*          will contain the diagonal and superdiagonal elements of a */
 137 /*          bidiagonal matrix orthogonally equivalent to the one given */
 138 /*          as input. */
 139
 140 /*  VT      (input/output) REAL array, dimension (LDVT, NCVT) */
 141 /*          On entry, an N-by-NCVT matrix VT. */
 142 /*          On exit, VT is overwritten by P**T * VT. */
 143 /*          Not referenced if NCVT = 0. */
 144
 145 /*  LDVT    (input) INTEGER */
 146 /*          The leading dimension of the array VT. */
 147 /*          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. */
 148
 149 /*  U       (input/output) REAL array, dimension (LDU, N) */
 150 /*          On entry, an NRU-by-N matrix U. */
 151 /*          On exit, U is overwritten by U * Q. */
 152 /*          Not referenced if NRU = 0. */
 153
 154 /*  LDU     (input) INTEGER */
 155 /*          The leading dimension of the array U.  LDU >= max(1,NRU). */
 156
 157 /*  C       (input/output) REAL array, dimension (LDC, NCC) */
 158 /*          On entry, an N-by-NCC matrix C. */
 159 /*          On exit, C is overwritten by Q**T * C. */
 160 /*          Not referenced if NCC = 0. */
 161
 162 /*  LDC     (input) INTEGER */
 163 /*          The leading dimension of the array C. */
 164 /*          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. */
 165
 166 /*  WORK    (workspace) REAL array, dimension (2*N) */
 167 /*          if NCVT = NRU = NCC = 0, (max(1, 4*N)) otherwise */
 168
 169 /*  INFO    (output) INTEGER */
 170 /*          = 0:  successful exit */
 171 /*          < 0:  If INFO = -i, the i-th argument had an illegal value */
 172 /*          > 0:  the algorithm did not converge; D and E contain the */
 173 /*                elements of a bidiagonal matrix which is orthogonally */
 174 /*                similar to the input matrix B;  if INFO = i, i */
 175 /*                elements of E have not converged to zero. */
 176
 177 /*  Internal Parameters */
 178 /*  =================== */
 179
 180 /*  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8))) */
 181 /*          TOLMUL controls the convergence criterion of the QR loop. */
 182 /*          If it is positive, TOLMUL*EPS is the desired relative */
 183 /*             precision in the computed singular values. */
 184 /*          If it is negative, abs(TOLMUL*EPS*sigma_max) is the */
 185 /*             desired absolute accuracy in the computed singular */
 186 /*             values (corresponds to relative accuracy */
 187 /*             abs(TOLMUL*EPS) in the largest singular value. */
 188 /*          abs(TOLMUL) should be between 1 and 1/EPS, and preferably */
 189 /*             between 10 (for fast convergence) and .1/EPS */
 190 /*             (for there to be some accuracy in the results). */
 191 /*          Default is to lose at either one eighth or 2 of the */
 192 /*             available decimal digits in each computed singular value */
 193 /*             (whichever is smaller). */
 194
 195 /*  MAXITR  INTEGER, default = 6 */
 196 /*          MAXITR controls the maximum number of passes of the */
 197 /*          algorithm through its inner loop. The algorithms stops */
 198 /*          (and so fails to converge) if the number of passes */
 199 /*          through the inner loop exceeds MAXITR*N**2. */
 200
 201 /*  ===================================================================== */
 202
 203 /*     .. Parameters .. */
 204 /*     .. */
 205 /*     .. Local Scalars .. */
 206 /*     .. */
 207 /*     .. External Functions .. */
 208 /*     .. */
 209 /*     .. External Subroutines .. */
 210 /*     .. */
 211 /*     .. Intrinsic Functions .. */
 212 /*     .. */
 213 /*     .. Executable Statements .. */
 214
 215 /*     Test the input parameters. */
 216
 217     /* Parameter adjustments */
 218     --d__;
 219     --e;
 220     vt_dim1 = *ldvt;
 221     vt_offset = 1 + vt_dim1;
 222     vt -= vt_offset;
 223     u_dim1 = *ldu;
 224     u_offset = 1 + u_dim1;
 225     u -= u_offset;
 226     c_dim1 = *ldc;
 227     c_offset = 1 + c_dim1;
 228     c__ -= c_offset;
 229     --work;
 230
 231     /* Function Body */
 232     *info = 0;
 233     lower = lsame_(uplo, "L");
 234     if (! lsame_(uplo, "U") && ! lower) {
 235         *info = -1;
 236     } else if (*n < 0) {
 237         *info = -2;
 238     } else if (*ncvt < 0) {
 239         *info = -3;
 240     } else if (*nru < 0) {
 241         *info = -4;
 242     } else if (*ncc < 0) {
 243         *info = -5;
 244     } else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n)) {
 245         *info = -9;
 246     } else if (*ldu < max(1,*nru)) {
 247         *info = -11;
 248     } else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n)) {
 249         *info = -13;
 250     }
 251     if (*info != 0) {
 252         i__1 = -(*info);
 253         xerbla_("SBDSQR", &i__1);
 254         return 0;
 255     }
 256     if (*n == 0) {
 257         return 0;
 258     }
 259     if (*n == 1) {
 260         goto L160;
 261     }
 262
 263 /*     ROTATE is true if any singular vectors desired, false otherwise */
 264
 265     rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
 266
 267 /*     If no singular vectors desired, use qd algorithm */
 268
 269     if (! rotate) {
 270         slasq1_(n, &d__[1], &e[1], &work[1], info);
 271         return 0;
 272     }
 273
 274     nm1 = *n - 1;
 275     nm12 = nm1 + nm1;
 276     nm13 = nm12 + nm1;
 277     idir = 0;
 278
 279 /*     Get machine constants */
 280
 281     eps = slamch_("Epsilon");
 282     unfl = slamch_("Safe minimum");
 283
 284 /*     If matrix lower bidiagonal, rotate to be upper bidiagonal */
 285 /*     by applying Givens rotations on the left */
 286
 287     if (lower) {
 288         i__1 = *n - 1;
 289         for (i__ = 1; i__ <= i__1; ++i__) {
 290             slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
 291             d__[i__] = r__;
 292             e[i__] = sn * d__[i__ + 1];
 293             d__[i__ + 1] = cs * d__[i__ + 1];
 294             work[i__] = cs;
 295             work[nm1 + i__] = sn;
 296 /* L10: */
 297         }
 298
 299 /*        Update singular vectors if desired */
 300
 301         if (*nru > 0) {
 302             slasr_("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset],
 303                     ldu);
 304         }
 305         if (*ncc > 0) {
 306             slasr_("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset],
 307                      ldc);
 308         }
 309     }
 310
 311 /*     Compute singular values to relative accuracy TOL */
 312 /*     (By setting TOL to be negative, algorithm will compute */
 313 /*     singular values to absolute accuracy ABS(TOL)*norm(input matrix)) */
 314
 315 /* Computing MAX */
 316 /* Computing MIN */
 317     d__1 = (doublereal) eps;
 318     r__3 = 100.f, r__4 = pow_dd(&d__1, &c_b15);
 319     r__1 = 10.f, r__2 = dmin(r__3,r__4);
 320     tolmul = dmax(r__1,r__2);
 321     tol = tolmul * eps;
 322
 323 /*     Compute approximate maximum, minimum singular values */
 324
 325     smax = 0.f;
 326     i__1 = *n;
 327     for (i__ = 1; i__ <= i__1; ++i__) {
 328 /* Computing MAX */
 329         r__2 = smax, r__3 = (r__1 = d__[i__], dabs(r__1));
 330         smax = dmax(r__2,r__3);
 331 /* L20: */
 332     }
 333     i__1 = *n - 1;
 334     for (i__ = 1; i__ <= i__1; ++i__) {
 335 /* Computing MAX */
 336         r__2 = smax, r__3 = (r__1 = e[i__], dabs(r__1));
 337         smax = dmax(r__2,r__3);
 338 /* L30: */
 339     }
 340     sminl = 0.f;
 341     if (tol >= 0.f) {
 342
 343 /*        Relative accuracy desired */
 344
 345         sminoa = dabs(d__[1]);
 346         if (sminoa == 0.f) {
 347             goto L50;
 348         }
 349         mu = sminoa;
 350         i__1 = *n;
 351         for (i__ = 2; i__ <= i__1; ++i__) {
 352             mu = (r__2 = d__[i__], dabs(r__2)) * (mu / (mu + (r__1 = e[i__ -
 353                     1], dabs(r__1))));
 354             sminoa = dmin(sminoa,mu);
 355             if (sminoa == 0.f) {
 356                 goto L50;
 357             }
 358 /* L40: */
 359         }
 360 L50:
 361         sminoa /= sqrt((real) (*n));
 362 /* Computing MAX */
 363         r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl;
 364         thresh = dmax(r__1,r__2);
 365     } else {
 366
 367 /*        Absolute accuracy desired */
 368
 369 /* Computing MAX */
 370         r__1 = dabs(tol) * smax, r__2 = *n * 6 * *n * unfl;
 371         thresh = dmax(r__1,r__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.f && (r__1 = d__[m], dabs(r__1)) <= thresh) {
 403         d__[m] = 0.f;
 404     }
 405     smax = (r__1 = d__[m], dabs(r__1));
 406     smin = smax;
 407     i__1 = m - 1;
 408     for (lll = 1; lll <= i__1; ++lll) {
 409         ll = m - lll;
 410         abss = (r__1 = d__[ll], dabs(r__1));
 411         abse = (r__1 = e[ll], dabs(r__1));
 412         if (tol < 0.f && abss <= thresh) {
 413             d__[ll] = 0.f;
 414         }
 415         if (abse <= thresh) {
 416             goto L80;
 417         }
 418         smin = dmin(smin,abss);
 419 /* Computing MAX */
 420         r__1 = max(smax,abss);
 421         smax = dmax(r__1,abse);
 422 /* L70: */
 423     }
 424     ll = 0;
 425     goto L90;
 426 L80:
 427     e[ll] = 0.f;
 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         slasv2_(&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.f;
 451         d__[m] = sigmn;
 452
 453 /*        Compute singular vectors, if desired */
 454
 455         if (*ncvt > 0) {
 456             srot_(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &
 457                     cosr, &sinr);
 458         }
 459         if (*nru > 0) {
 460             srot_(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             srot_(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 ((r__1 = d__[ll], dabs(r__1)) >= (r__2 = d__[m], dabs(r__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 ((r__2 = e[m - 1], dabs(r__2)) <= dabs(tol) * (r__1 = d__[m], dabs(
 496                 r__1)) || tol < 0.f && (r__3 = e[m - 1], dabs(r__3)) <=
 497                 thresh) {
 498             e[m - 1] = 0.f;
 499             goto L60;
 500         }
 501
 502         if (tol >= 0.f) {
 503
 504 /*           If relative accuracy desired, */
 505 /*           apply convergence criterion forward */
 506
 507             mu = (r__1 = d__[ll], dabs(r__1));
 508             sminl = mu;
 509             i__1 = m - 1;
 510             for (lll = ll; lll <= i__1; ++lll) {
 511                 if ((r__1 = e[lll], dabs(r__1)) <= tol * mu) {
 512                     e[lll] = 0.f;
 513                     goto L60;
 514                 }
 515                 mu = (r__2 = d__[lll + 1], dabs(r__2)) * (mu / (mu + (r__1 =
 516                         e[lll], dabs(r__1))));
 517                 sminl = dmin(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 ((r__2 = e[ll], dabs(r__2)) <= dabs(tol) * (r__1 = d__[ll], dabs(
 528                 r__1)) || tol < 0.f && (r__3 = e[ll], dabs(r__3)) <= thresh) {
 529             e[ll] = 0.f;
 530             goto L60;
 531         }
 532
 533         if (tol >= 0.f) {
 534
 535 /*           If relative accuracy desired, */
 536 /*           apply convergence criterion backward */
 537
 538             mu = (r__1 = d__[m], dabs(r__1));
 539             sminl = mu;
 540             i__1 = ll;
 541             for (lll = m - 1; lll >= i__1; --lll) {
 542                 if ((r__1 = e[lll], dabs(r__1)) <= tol * mu) {
 543                     e[lll] = 0.f;
 544                     goto L60;
 545                 }
 546                 mu = (r__2 = d__[lll], dabs(r__2)) * (mu / (mu + (r__1 = e[
 547                         lll], dabs(r__1))));
 548                 sminl = dmin(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     r__1 = eps, r__2 = tol * .01f;
 561     if (tol >= 0.f && *n * tol * (sminl / smax) <= dmax(r__1,r__2)) {
 562
 563 /*        Use a zero shift to avoid loss of relative accuracy */
 564
 565         shift = 0.f;
 566     } else {
 567
 568 /*        Compute the shift from 2-by-2 block at end of matrix */
 569
 570         if (idir == 1) {
 571             sll = (r__1 = d__[ll], dabs(r__1));
 572             slas2_(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
 573         } else {
 574             sll = (r__1 = d__[m], dabs(r__1));
 575             slas2_(&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.f) {
 581 /* Computing 2nd power */
 582             r__1 = shift / sll;
 583             if (r__1 * r__1 < eps) {
 584                 shift = 0.f;
 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.f) {
 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.f;
 602             oldcs = 1.f;
 603             i__1 = m - 1;
 604             for (i__ = ll; i__ <= i__1; ++i__) {
 605                 r__1 = d__[i__] * cs;
 606                 slartg_(&r__1, &e[i__], &cs, &sn, &r__);
 607                 if (i__ > ll) {
 608                     e[i__ - 1] = oldsn * r__;
 609                 }
 610                 r__1 = oldcs * r__;
 611                 r__2 = d__[i__ + 1] * sn;
 612                 slartg_(&r__1, &r__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                 slasr_("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                 slasr_("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                 slasr_("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 ((r__1 = e[m - 1], dabs(r__1)) <= thresh) {
 644                 e[m - 1] = 0.f;
 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.f;
 653             oldcs = 1.f;
 654             i__1 = ll + 1;
 655             for (i__ = m; i__ >= i__1; --i__) {
 656                 r__1 = d__[i__] * cs;
 657                 slartg_(&r__1, &e[i__ - 1], &cs, &sn, &r__);
 658                 if (i__ < m) {
 659                     e[i__] = oldsn * r__;
 660                 }
 661                 r__1 = oldcs * r__;
 662                 r__2 = d__[i__ - 1] * sn;
 663                 slartg_(&r__1, &r__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                 slasr_("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                 slasr_("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                 slasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
 689                         ll + c_dim1], ldc);
 690             }
 691
 692 /*           Test convergence */
 693
 694             if ((r__1 = e[ll], dabs(r__1)) <= thresh) {
 695                 e[ll] = 0.f;
 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 = ((r__1 = d__[ll], dabs(r__1)) - shift) * (r_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                 slartg_(&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                 slartg_(&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                 slasr_("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                 slasr_("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                 slasr_("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 ((r__1 = e[m - 1], dabs(r__1)) <= thresh) {
 757                 e[m - 1] = 0.f;
 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 = ((r__1 = d__[m], dabs(r__1)) - shift) * (r_sign(&c_b49, &d__[
 766                     m]) + shift / d__[m]);
 767             g = e[m - 1];
 768             i__1 = ll + 1;
 769             for (i__ = m; i__ >= i__1; --i__) {
 770                 slartg_(&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                 slartg_(&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 ((r__1 = e[ll], dabs(r__1)) <= thresh) {
 797                 e[ll] = 0.f;
 798             }
 799
 800 /*           Update singular vectors if desired */
 801
 802             if (*ncvt > 0) {
 803                 i__1 = m - ll + 1;
 804                 slasr_("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                 slasr_("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                 slasr_("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.f) {
 830             d__[i__] = -d__[i__];
 831
 832 /*           Change sign of singular vectors, if desired */
 833
 834             if (*ncvt > 0) {
 835                 sscal_(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                 sswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ +
 867                         vt_dim1], ldvt);
 868             }
 869             if (*nru > 0) {
 870                 sswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) *
 871                         u_dim1 + 1], &c__1);
 872             }
 873             if (*ncc > 0) {
 874                 sswap_(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.f) {
 889             ++(*info);
 890         }
 891 /* L210: */
 892     }
 893 L220:
 894     return 0;
 895
 896 /*     End of SBDSQR */
 897
 898 } /* sbdsqr_ */