LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ cbdsqr()

 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

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.```

Definition at line 220 of file cbdsqr.f.

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