LAPACK 3.11.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

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

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:111
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: