LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
cbdsqr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cbdsqr (UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
 CBDSQR More...
 

Function/Subroutine Documentation

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 (2*N)
          if NCVT = NRU = NCC = 0, (max(1, 4*N-4)) otherwise
[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 2011

Definition at line 225 of file cbdsqr.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: