LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cbdsqr ( character  UPLO,
integer  N,
integer  NCVT,
integer  NRU,
integer  NCC,
real, dimension( * )  D,
real, dimension( * )  E,
complex, dimension( ldvt, * )  VT,
integer  LDVT,
complex, dimension( ldu, * )  U,
integer  LDU,
complex, dimension( ldc, * )  C,
integer  LDC,
real, dimension( * )  RWORK,
integer  INFO 
)

CBDSQR

Download CBDSQR + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CBDSQR computes the singular values and, optionally, the right and/or
 left singular vectors from the singular value decomposition (SVD) of
 a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
 zero-shift QR algorithm.  The SVD of B has the form
 
    B = Q * S * P**H
 
 where S is the diagonal matrix of singular values, Q is an orthogonal
 matrix of left singular vectors, and P is an orthogonal matrix of
 right singular vectors.  If left singular vectors are requested, this
 subroutine actually returns U*Q instead of Q, and, if right singular
 vectors are requested, this subroutine returns P**H*VT instead of
 P**H, for given complex input matrices U and VT.  When U and VT are
 the unitary matrices that reduce a general matrix A to bidiagonal
 form: A = U*B*VT, as computed by CGEBRD, then
 
    A = (U*Q) * S * (P**H*VT)
 
 is the SVD of A.  Optionally, the subroutine may also compute Q**H*C
 for a given complex input matrix C.

 See "Computing  Small Singular Values of Bidiagonal Matrices With
 Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
 LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
 no. 5, pp. 873-912, Sept 1990) and
 "Accurate singular values and differential qd algorithms," by
 B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
 Department, University of California at Berkeley, July 1992
 for a detailed description of the algorithm.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  B is upper bidiagonal;
          = 'L':  B is lower bidiagonal.
[in]N
          N is INTEGER
          The order of the matrix B.  N >= 0.
[in]NCVT
          NCVT is INTEGER
          The number of columns of the matrix VT. NCVT >= 0.
[in]NRU
          NRU is INTEGER
          The number of rows of the matrix U. NRU >= 0.
[in]NCC
          NCC is INTEGER
          The number of columns of the matrix C. NCC >= 0.
[in,out]D
          D is REAL array, dimension (N)
          On entry, the n diagonal elements of the bidiagonal matrix B.
          On exit, if INFO=0, the singular values of B in decreasing
          order.
[in,out]E
          E is REAL array, dimension (N-1)
          On entry, the N-1 offdiagonal elements of the bidiagonal
          matrix B.
          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
          will contain the diagonal and superdiagonal elements of a
          bidiagonal matrix orthogonally equivalent to the one given
          as input.
[in,out]VT
          VT is COMPLEX array, dimension (LDVT, NCVT)
          On entry, an N-by-NCVT matrix VT.
          On exit, VT is overwritten by P**H * VT.
          Not referenced if NCVT = 0.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.
          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
[in,out]U
          U is COMPLEX array, dimension (LDU, N)
          On entry, an NRU-by-N matrix U.
          On exit, U is overwritten by U * Q.
          Not referenced if NRU = 0.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= max(1,NRU).
[in,out]C
          C is COMPLEX array, dimension (LDC, NCC)
          On entry, an N-by-NCC matrix C.
          On exit, C is overwritten by Q**H * C.
          Not referenced if NCC = 0.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C.
          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
[out]RWORK
          RWORK is REAL array, dimension (4*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  If INFO = -i, the i-th argument had an illegal value
          > 0:  the algorithm did not converge; D and E contain the
                elements of a bidiagonal matrix which is orthogonally
                similar to the input matrix B;  if INFO = i, i
                elements of E have not converged to zero.
Internal Parameters:
  TOLMUL  REAL, default = max(10,min(100,EPS**(-1/8)))
          TOLMUL controls the convergence criterion of the QR loop.
          If it is positive, TOLMUL*EPS is the desired relative
             precision in the computed singular values.
          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
             desired absolute accuracy in the computed singular
             values (corresponds to relative accuracy
             abs(TOLMUL*EPS) in the largest singular value.
          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
             between 10 (for fast convergence) and .1/EPS
             (for there to be some accuracy in the results).
          Default is to lose at either one eighth or 2 of the
             available decimal digits in each computed singular value
             (whichever is smaller).

  MAXITR  INTEGER, default = 6
          MAXITR controls the maximum number of passes of the
          algorithm through its inner loop. The algorithms stops
          (and so fails to converge) if the number of passes
          through the inner loop exceeds MAXITR*N**2.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 224 of file cbdsqr.f.

224 *
225 * -- LAPACK computational routine (version 3.6.0) --
226 * -- LAPACK is a software package provided by Univ. of Tennessee, --
227 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
228 * November 2015
229 *
230 * .. Scalar Arguments ..
231  CHARACTER uplo
232  INTEGER info, ldc, ldu, ldvt, n, ncc, ncvt, nru
233 * ..
234 * .. Array Arguments ..
235  REAL d( * ), e( * ), rwork( * )
236  COMPLEX c( ldc, * ), u( ldu, * ), vt( ldvt, * )
237 * ..
238 *
239 * =====================================================================
240 *
241 * .. Parameters ..
242  REAL zero
243  parameter ( zero = 0.0e0 )
244  REAL one
245  parameter ( one = 1.0e0 )
246  REAL negone
247  parameter ( negone = -1.0e0 )
248  REAL hndrth
249  parameter ( hndrth = 0.01e0 )
250  REAL ten
251  parameter ( ten = 10.0e0 )
252  REAL hndrd
253  parameter ( hndrd = 100.0e0 )
254  REAL meigth
255  parameter ( meigth = -0.125e0 )
256  INTEGER maxitr
257  parameter ( maxitr = 6 )
258 * ..
259 * .. Local Scalars ..
260  LOGICAL lower, rotate
261  INTEGER i, idir, isub, iter, j, ll, lll, m, maxit, nm1,
262  $ nm12, nm13, oldll, oldm
263  REAL abse, abss, cosl, cosr, cs, eps, f, g, h, mu,
264  $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
265  $ sinr, sll, smax, smin, sminl, sminoa,
266  $ sn, thresh, tol, tolmul, unfl
267 * ..
268 * .. External Functions ..
269  LOGICAL lsame
270  REAL slamch
271  EXTERNAL lsame, slamch
272 * ..
273 * .. External Subroutines ..
274  EXTERNAL clasr, csrot, csscal, cswap, slartg, slas2,
275  $ slasq1, slasv2, xerbla
276 * ..
277 * .. Intrinsic Functions ..
278  INTRINSIC abs, max, min, REAL, sign, sqrt
279 * ..
280 * .. Executable Statements ..
281 *
282 * Test the input parameters.
283 *
284  info = 0
285  lower = lsame( uplo, 'L' )
286  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
287  info = -1
288  ELSE IF( n.LT.0 ) THEN
289  info = -2
290  ELSE IF( ncvt.LT.0 ) THEN
291  info = -3
292  ELSE IF( nru.LT.0 ) THEN
293  info = -4
294  ELSE IF( ncc.LT.0 ) THEN
295  info = -5
296  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
297  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
298  info = -9
299  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
300  info = -11
301  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
302  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
303  info = -13
304  END IF
305  IF( info.NE.0 ) THEN
306  CALL xerbla( 'CBDSQR', -info )
307  RETURN
308  END IF
309  IF( n.EQ.0 )
310  $ RETURN
311  IF( n.EQ.1 )
312  $ GO TO 160
313 *
314 * ROTATE is true if any singular vectors desired, false otherwise
315 *
316  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
317 *
318 * If no singular vectors desired, use qd algorithm
319 *
320  IF( .NOT.rotate ) THEN
321  CALL slasq1( n, d, e, rwork, info )
322 *
323 * If INFO equals 2, dqds didn't finish, try to finish
324 *
325  IF( info .NE. 2 ) RETURN
326  info = 0
327  END IF
328 *
329  nm1 = n - 1
330  nm12 = nm1 + nm1
331  nm13 = nm12 + nm1
332  idir = 0
333 *
334 * Get machine constants
335 *
336  eps = slamch( 'Epsilon' )
337  unfl = slamch( 'Safe minimum' )
338 *
339 * If matrix lower bidiagonal, rotate to be upper bidiagonal
340 * by applying Givens rotations on the left
341 *
342  IF( lower ) THEN
343  DO 10 i = 1, n - 1
344  CALL slartg( d( i ), e( i ), cs, sn, r )
345  d( i ) = r
346  e( i ) = sn*d( i+1 )
347  d( i+1 ) = cs*d( i+1 )
348  rwork( i ) = cs
349  rwork( nm1+i ) = sn
350  10 CONTINUE
351 *
352 * Update singular vectors if desired
353 *
354  IF( nru.GT.0 )
355  $ CALL clasr( 'R', 'V', 'F', nru, n, rwork( 1 ), rwork( n ),
356  $ u, ldu )
357  IF( ncc.GT.0 )
358  $ CALL clasr( 'L', 'V', 'F', n, ncc, rwork( 1 ), rwork( n ),
359  $ c, ldc )
360  END IF
361 *
362 * Compute singular values to relative accuracy TOL
363 * (By setting TOL to be negative, algorithm will compute
364 * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
365 *
366  tolmul = max( ten, min( hndrd, eps**meigth ) )
367  tol = tolmul*eps
368 *
369 * Compute approximate maximum, minimum singular values
370 *
371  smax = zero
372  DO 20 i = 1, n
373  smax = max( smax, abs( d( i ) ) )
374  20 CONTINUE
375  DO 30 i = 1, n - 1
376  smax = max( smax, abs( e( i ) ) )
377  30 CONTINUE
378  sminl = zero
379  IF( tol.GE.zero ) THEN
380 *
381 * Relative accuracy desired
382 *
383  sminoa = abs( d( 1 ) )
384  IF( sminoa.EQ.zero )
385  $ GO TO 50
386  mu = sminoa
387  DO 40 i = 2, n
388  mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
389  sminoa = min( sminoa, mu )
390  IF( sminoa.EQ.zero )
391  $ GO TO 50
392  40 CONTINUE
393  50 CONTINUE
394  sminoa = sminoa / sqrt( REAL( N ) )
395  thresh = max( tol*sminoa, maxitr*n*n*unfl )
396  ELSE
397 *
398 * Absolute accuracy desired
399 *
400  thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
401  END IF
402 *
403 * Prepare for main iteration loop for the singular values
404 * (MAXIT is the maximum number of passes through the inner
405 * loop permitted before nonconvergence signalled.)
406 *
407  maxit = maxitr*n*n
408  iter = 0
409  oldll = -1
410  oldm = -1
411 *
412 * M points to last element of unconverged part of matrix
413 *
414  m = n
415 *
416 * Begin main iteration loop
417 *
418  60 CONTINUE
419 *
420 * Check for convergence or exceeding iteration count
421 *
422  IF( m.LE.1 )
423  $ GO TO 160
424  IF( iter.GT.maxit )
425  $ GO TO 200
426 *
427 * Find diagonal block of matrix to work on
428 *
429  IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
430  $ d( m ) = zero
431  smax = abs( d( m ) )
432  smin = smax
433  DO 70 lll = 1, m - 1
434  ll = m - lll
435  abss = abs( d( ll ) )
436  abse = abs( e( ll ) )
437  IF( tol.LT.zero .AND. abss.LE.thresh )
438  $ d( ll ) = zero
439  IF( abse.LE.thresh )
440  $ GO TO 80
441  smin = min( smin, abss )
442  smax = max( smax, abss, abse )
443  70 CONTINUE
444  ll = 0
445  GO TO 90
446  80 CONTINUE
447  e( ll ) = zero
448 *
449 * Matrix splits since E(LL) = 0
450 *
451  IF( ll.EQ.m-1 ) THEN
452 *
453 * Convergence of bottom singular value, return to top of loop
454 *
455  m = m - 1
456  GO TO 60
457  END IF
458  90 CONTINUE
459  ll = ll + 1
460 *
461 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
462 *
463  IF( ll.EQ.m-1 ) THEN
464 *
465 * 2 by 2 block, handle separately
466 *
467  CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
468  $ cosr, sinl, cosl )
469  d( m-1 ) = sigmx
470  e( m-1 ) = zero
471  d( m ) = sigmn
472 *
473 * Compute singular vectors, if desired
474 *
475  IF( ncvt.GT.0 )
476  $ CALL csrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
477  $ cosr, sinr )
478  IF( nru.GT.0 )
479  $ CALL csrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
480  IF( ncc.GT.0 )
481  $ CALL csrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
482  $ sinl )
483  m = m - 2
484  GO TO 60
485  END IF
486 *
487 * If working on new submatrix, choose shift direction
488 * (from larger end diagonal element towards smaller)
489 *
490  IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
491  IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
492 *
493 * Chase bulge from top (big end) to bottom (small end)
494 *
495  idir = 1
496  ELSE
497 *
498 * Chase bulge from bottom (big end) to top (small end)
499 *
500  idir = 2
501  END IF
502  END IF
503 *
504 * Apply convergence tests
505 *
506  IF( idir.EQ.1 ) THEN
507 *
508 * Run convergence test in forward direction
509 * First apply standard test to bottom of matrix
510 *
511  IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
512  $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
513  e( m-1 ) = zero
514  GO TO 60
515  END IF
516 *
517  IF( tol.GE.zero ) THEN
518 *
519 * If relative accuracy desired,
520 * apply convergence criterion forward
521 *
522  mu = abs( d( ll ) )
523  sminl = mu
524  DO 100 lll = ll, m - 1
525  IF( abs( e( lll ) ).LE.tol*mu ) THEN
526  e( lll ) = zero
527  GO TO 60
528  END IF
529  mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
530  sminl = min( sminl, mu )
531  100 CONTINUE
532  END IF
533 *
534  ELSE
535 *
536 * Run convergence test in backward direction
537 * First apply standard test to top of matrix
538 *
539  IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
540  $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
541  e( ll ) = zero
542  GO TO 60
543  END IF
544 *
545  IF( tol.GE.zero ) THEN
546 *
547 * If relative accuracy desired,
548 * apply convergence criterion backward
549 *
550  mu = abs( d( m ) )
551  sminl = mu
552  DO 110 lll = m - 1, ll, -1
553  IF( abs( e( lll ) ).LE.tol*mu ) THEN
554  e( lll ) = zero
555  GO TO 60
556  END IF
557  mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
558  sminl = min( sminl, mu )
559  110 CONTINUE
560  END IF
561  END IF
562  oldll = ll
563  oldm = m
564 *
565 * Compute shift. First, test if shifting would ruin relative
566 * accuracy, and if so set the shift to zero.
567 *
568  IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
569  $ max( eps, hndrth*tol ) ) THEN
570 *
571 * Use a zero shift to avoid loss of relative accuracy
572 *
573  shift = zero
574  ELSE
575 *
576 * Compute the shift from 2-by-2 block at end of matrix
577 *
578  IF( idir.EQ.1 ) THEN
579  sll = abs( d( ll ) )
580  CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
581  ELSE
582  sll = abs( d( m ) )
583  CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
584  END IF
585 *
586 * Test if shift negligible, and if so set to zero
587 *
588  IF( sll.GT.zero ) THEN
589  IF( ( shift / sll )**2.LT.eps )
590  $ shift = zero
591  END IF
592  END IF
593 *
594 * Increment iteration count
595 *
596  iter = iter + m - ll
597 *
598 * If SHIFT = 0, do simplified QR iteration
599 *
600  IF( shift.EQ.zero ) THEN
601  IF( idir.EQ.1 ) THEN
602 *
603 * Chase bulge from top to bottom
604 * Save cosines and sines for later singular vector updates
605 *
606  cs = one
607  oldcs = one
608  DO 120 i = ll, m - 1
609  CALL slartg( d( i )*cs, e( i ), cs, sn, r )
610  IF( i.GT.ll )
611  $ e( i-1 ) = oldsn*r
612  CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
613  rwork( i-ll+1 ) = cs
614  rwork( i-ll+1+nm1 ) = sn
615  rwork( i-ll+1+nm12 ) = oldcs
616  rwork( i-ll+1+nm13 ) = oldsn
617  120 CONTINUE
618  h = d( m )*cs
619  d( m ) = h*oldcs
620  e( m-1 ) = h*oldsn
621 *
622 * Update singular vectors
623 *
624  IF( ncvt.GT.0 )
625  $ CALL clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
626  $ rwork( n ), vt( ll, 1 ), ldvt )
627  IF( nru.GT.0 )
628  $ CALL clasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),
629  $ rwork( nm13+1 ), u( 1, ll ), ldu )
630  IF( ncc.GT.0 )
631  $ CALL clasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),
632  $ rwork( nm13+1 ), c( ll, 1 ), ldc )
633 *
634 * Test convergence
635 *
636  IF( abs( e( m-1 ) ).LE.thresh )
637  $ e( m-1 ) = zero
638 *
639  ELSE
640 *
641 * Chase bulge from bottom to top
642 * Save cosines and sines for later singular vector updates
643 *
644  cs = one
645  oldcs = one
646  DO 130 i = m, ll + 1, -1
647  CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
648  IF( i.LT.m )
649  $ e( i ) = oldsn*r
650  CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
651  rwork( i-ll ) = cs
652  rwork( i-ll+nm1 ) = -sn
653  rwork( i-ll+nm12 ) = oldcs
654  rwork( i-ll+nm13 ) = -oldsn
655  130 CONTINUE
656  h = d( ll )*cs
657  d( ll ) = h*oldcs
658  e( ll ) = h*oldsn
659 *
660 * Update singular vectors
661 *
662  IF( ncvt.GT.0 )
663  $ CALL clasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),
664  $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
665  IF( nru.GT.0 )
666  $ CALL clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
667  $ rwork( n ), u( 1, ll ), ldu )
668  IF( ncc.GT.0 )
669  $ CALL clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
670  $ rwork( n ), c( ll, 1 ), ldc )
671 *
672 * Test convergence
673 *
674  IF( abs( e( ll ) ).LE.thresh )
675  $ e( ll ) = zero
676  END IF
677  ELSE
678 *
679 * Use nonzero shift
680 *
681  IF( idir.EQ.1 ) THEN
682 *
683 * Chase bulge from top to bottom
684 * Save cosines and sines for later singular vector updates
685 *
686  f = ( abs( d( ll ) )-shift )*
687  $ ( sign( one, d( ll ) )+shift / d( ll ) )
688  g = e( ll )
689  DO 140 i = ll, m - 1
690  CALL slartg( f, g, cosr, sinr, r )
691  IF( i.GT.ll )
692  $ e( i-1 ) = r
693  f = cosr*d( i ) + sinr*e( i )
694  e( i ) = cosr*e( i ) - sinr*d( i )
695  g = sinr*d( i+1 )
696  d( i+1 ) = cosr*d( i+1 )
697  CALL slartg( f, g, cosl, sinl, r )
698  d( i ) = r
699  f = cosl*e( i ) + sinl*d( i+1 )
700  d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
701  IF( i.LT.m-1 ) THEN
702  g = sinl*e( i+1 )
703  e( i+1 ) = cosl*e( i+1 )
704  END IF
705  rwork( i-ll+1 ) = cosr
706  rwork( i-ll+1+nm1 ) = sinr
707  rwork( i-ll+1+nm12 ) = cosl
708  rwork( i-ll+1+nm13 ) = sinl
709  140 CONTINUE
710  e( m-1 ) = f
711 *
712 * Update singular vectors
713 *
714  IF( ncvt.GT.0 )
715  $ CALL clasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
716  $ rwork( n ), vt( ll, 1 ), ldvt )
717  IF( nru.GT.0 )
718  $ CALL clasr( 'R', 'V', 'F', nru, m-ll+1, rwork( nm12+1 ),
719  $ rwork( nm13+1 ), u( 1, ll ), ldu )
720  IF( ncc.GT.0 )
721  $ CALL clasr( 'L', 'V', 'F', m-ll+1, ncc, rwork( nm12+1 ),
722  $ rwork( nm13+1 ), c( ll, 1 ), ldc )
723 *
724 * Test convergence
725 *
726  IF( abs( e( m-1 ) ).LE.thresh )
727  $ e( m-1 ) = zero
728 *
729  ELSE
730 *
731 * Chase bulge from bottom to top
732 * Save cosines and sines for later singular vector updates
733 *
734  f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
735  $ d( m ) )
736  g = e( m-1 )
737  DO 150 i = m, ll + 1, -1
738  CALL slartg( f, g, cosr, sinr, r )
739  IF( i.LT.m )
740  $ e( i ) = r
741  f = cosr*d( i ) + sinr*e( i-1 )
742  e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
743  g = sinr*d( i-1 )
744  d( i-1 ) = cosr*d( i-1 )
745  CALL slartg( f, g, cosl, sinl, r )
746  d( i ) = r
747  f = cosl*e( i-1 ) + sinl*d( i-1 )
748  d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
749  IF( i.GT.ll+1 ) THEN
750  g = sinl*e( i-2 )
751  e( i-2 ) = cosl*e( i-2 )
752  END IF
753  rwork( i-ll ) = cosr
754  rwork( i-ll+nm1 ) = -sinr
755  rwork( i-ll+nm12 ) = cosl
756  rwork( i-ll+nm13 ) = -sinl
757  150 CONTINUE
758  e( ll ) = f
759 *
760 * Test convergence
761 *
762  IF( abs( e( ll ) ).LE.thresh )
763  $ e( ll ) = zero
764 *
765 * Update singular vectors if desired
766 *
767  IF( ncvt.GT.0 )
768  $ CALL clasr( 'L', 'V', 'B', m-ll+1, ncvt, rwork( nm12+1 ),
769  $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
770  IF( nru.GT.0 )
771  $ CALL clasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
772  $ rwork( n ), u( 1, ll ), ldu )
773  IF( ncc.GT.0 )
774  $ CALL clasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
775  $ rwork( n ), c( ll, 1 ), ldc )
776  END IF
777  END IF
778 *
779 * QR iteration finished, go back and check convergence
780 *
781  GO TO 60
782 *
783 * All singular values converged, so make them positive
784 *
785  160 CONTINUE
786  DO 170 i = 1, n
787  IF( d( i ).LT.zero ) THEN
788  d( i ) = -d( i )
789 *
790 * Change sign of singular vectors, if desired
791 *
792  IF( ncvt.GT.0 )
793  $ CALL csscal( ncvt, negone, vt( i, 1 ), ldvt )
794  END IF
795  170 CONTINUE
796 *
797 * Sort the singular values into decreasing order (insertion sort on
798 * singular values, but only one transposition per singular vector)
799 *
800  DO 190 i = 1, n - 1
801 *
802 * Scan for smallest D(I)
803 *
804  isub = 1
805  smin = d( 1 )
806  DO 180 j = 2, n + 1 - i
807  IF( d( j ).LE.smin ) THEN
808  isub = j
809  smin = d( j )
810  END IF
811  180 CONTINUE
812  IF( isub.NE.n+1-i ) THEN
813 *
814 * Swap singular values and vectors
815 *
816  d( isub ) = d( n+1-i )
817  d( n+1-i ) = smin
818  IF( ncvt.GT.0 )
819  $ CALL cswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
820  $ ldvt )
821  IF( nru.GT.0 )
822  $ CALL cswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
823  IF( ncc.GT.0 )
824  $ CALL cswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
825  END IF
826  190 CONTINUE
827  GO TO 220
828 *
829 * Maximum number of iterations exceeded, failure to converge
830 *
831  200 CONTINUE
832  info = 0
833  DO 210 i = 1, n - 1
834  IF( e( i ).NE.zero )
835  $ info = info + 1
836  210 CONTINUE
837  220 CONTINUE
838  RETURN
839 *
840 * End of CBDSQR
841 *
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:109
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99
subroutine slasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: slasv2.f:140
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: clasr.f:202
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:100
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine slasq1(N, D, E, WORK, INFO)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition: slasq1.f:110
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: