LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zbdsqr()

subroutine zbdsqr ( character  UPLO,
integer  N,
integer  NCVT,
integer  NRU,
integer  NCC,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
complex*16, dimension( ldvt, * )  VT,
integer  LDVT,
complex*16, dimension( ldu, * )  U,
integer  LDU,
complex*16, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZBDSQR

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

Purpose:
 ZBDSQR 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 ZGEBRD, 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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*16 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 DOUBLE PRECISION 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  DOUBLE PRECISION, 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.

Definition at line 220 of file zbdsqr.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  DOUBLE PRECISION D( * ), E( * ), RWORK( * )
233  COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
234 * ..
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239  DOUBLE PRECISION ZERO
240  parameter( zero = 0.0d0 )
241  DOUBLE PRECISION ONE
242  parameter( one = 1.0d0 )
243  DOUBLE PRECISION NEGONE
244  parameter( negone = -1.0d0 )
245  DOUBLE PRECISION HNDRTH
246  parameter( hndrth = 0.01d0 )
247  DOUBLE PRECISION TEN
248  parameter( ten = 10.0d0 )
249  DOUBLE PRECISION HNDRD
250  parameter( hndrd = 100.0d0 )
251  DOUBLE PRECISION MEIGTH
252  parameter( meigth = -0.125d0 )
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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
268  EXTERNAL lsame, dlamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL dlartg, dlas2, dlasq1, dlasv2, xerbla, zdrot,
272  $ zdscal, zlasr, zswap
273 * ..
274 * .. Intrinsic Functions ..
275  INTRINSIC abs, dble, max, min, 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( 'ZBDSQR', -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 dlasq1( 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 = dlamch( 'Epsilon' )
334  unfl = dlamch( '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 dlartg( 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 zlasr( 'R', 'V', 'F', nru, n, rwork( 1 ), rwork( n ),
353  $ u, ldu )
354  IF( ncc.GT.0 )
355  $ CALL zlasr( '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( dble( 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 dlasv2( 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 zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
474  $ cosr, sinr )
475  IF( nru.GT.0 )
476  $ CALL zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
477  IF( ncc.GT.0 )
478  $ CALL zdrot( 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 dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
578  ELSE
579  sll = abs( d( m ) )
580  CALL dlas2( 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 dlartg( d( i )*cs, e( i ), cs, sn, r )
607  IF( i.GT.ll )
608  $ e( i-1 ) = oldsn*r
609  CALL dlartg( 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 zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
623  $ rwork( n ), vt( ll, 1 ), ldvt )
624  IF( nru.GT.0 )
625  $ CALL zlasr( '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 zlasr( '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 dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
645  IF( i.LT.m )
646  $ e( i ) = oldsn*r
647  CALL dlartg( 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 zlasr( '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 zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
664  $ rwork( n ), u( 1, ll ), ldu )
665  IF( ncc.GT.0 )
666  $ CALL zlasr( '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 dlartg( 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 dlartg( 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 zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
713  $ rwork( n ), vt( ll, 1 ), ldvt )
714  IF( nru.GT.0 )
715  $ CALL zlasr( '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 zlasr( '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 dlartg( 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 dlartg( 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 zlasr( '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 zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
769  $ rwork( n ), u( 1, ll ), ldu )
770  IF( ncc.GT.0 )
771  $ CALL zlasr( '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 zdscal( 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 zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
817  $ ldvt )
818  IF( nru.GT.0 )
819  $ CALL zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
820  IF( ncc.GT.0 )
821  $ CALL zswap( 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 ZBDSQR
838 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: dlas2.f:107
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition: dlasv2.f:138
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dlasq1(N, D, E, WORK, INFO)
DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition: dlasq1.f:108
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zdrot(N, ZX, INCX, ZY, INCY, C, S)
ZDROT
Definition: zdrot.f:98
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: zlasr.f:200
Here is the call graph for this function:
Here is the caller graph for this function: