 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ sbdsqr()

 subroutine sbdsqr ( character UPLO, integer N, integer NCVT, integer NRU, integer NCC, real, dimension( * ) D, real, dimension( * ) E, real, dimension( ldvt, * ) VT, integer LDVT, real, dimension( ldu, * ) U, integer LDU, real, dimension( ldc, * ) C, integer LDC, real, dimension( * ) WORK, integer INFO )

SBDSQR

Purpose:
SBDSQR 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**T

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**T*VT instead of
P**T, for given real input matrices U and VT.  When U and VT are the
orthogonal matrices that reduce a general matrix A to bidiagonal
form:  A = U*B*VT, as computed by SGEBRD, then

A = (U*Q) * S * (P**T*VT)

is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
for a given real 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 REAL array, dimension (LDVT, NCVT) On entry, an N-by-NCVT matrix VT. On exit, VT is overwritten by P**T * 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 REAL 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 REAL array, dimension (LDC, NCC) On entry, an N-by-NCC matrix C. On exit, C is overwritten by Q**T * 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] WORK WORK 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: if NCVT = NRU = NCC = 0, = 1, a split was marked by a positive value in E = 2, current block of Z not diagonalized after 30*N iterations (in inner while loop) = 3, termination criterion of outer while loop not met (program created more than N unreduced blocks) else NCVT = NRU = NCC = 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.
Note:
Bug report from Cezary Dendek.
On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
removed since it can overflow pretty easily (for N larger or equal
than 18,919). We instead use MAXITDIVN = MAXITR*N.

Definition at line 238 of file sbdsqr.f.

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