LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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

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.```
Note:
```  Bug report from Cezary Dendek.
On November 3rd 2023, 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 231 of file zbdsqr.f.

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