LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
zbdsqr.f
Go to the documentation of this file.
1*> \brief \b ZBDSQR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZBDSQR + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbdsqr.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbdsqr.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbdsqr.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
20* LDU, C, LDC, RWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO
24* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION D( * ), E( * ), RWORK( * )
28* COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> ZBDSQR computes the singular values and, optionally, the right and/or
38*> left singular vectors from the singular value decomposition (SVD) of
39*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
40*> zero-shift QR algorithm. The SVD of B has the form
41*>
42*> B = Q * S * P**H
43*>
44*> where S is the diagonal matrix of singular values, Q is an orthogonal
45*> matrix of left singular vectors, and P is an orthogonal matrix of
46*> right singular vectors. If left singular vectors are requested, this
47*> subroutine actually returns U*Q instead of Q, and, if right singular
48*> vectors are requested, this subroutine returns P**H*VT instead of
49*> P**H, for given complex input matrices U and VT. When U and VT are
50*> the unitary matrices that reduce a general matrix A to bidiagonal
51*> form: A = U*B*VT, as computed by ZGEBRD, then
52*>
53*> A = (U*Q) * S * (P**H*VT)
54*>
55*> is the SVD of A. Optionally, the subroutine may also compute Q**H*C
56*> for a given complex input matrix C.
57*>
58*> See "Computing Small Singular Values of Bidiagonal Matrices With
59*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
60*> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
61*> no. 5, pp. 873-912, Sept 1990) and
62*> "Accurate singular values and differential qd algorithms," by
63*> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
64*> Department, University of California at Berkeley, July 1992
65*> for a detailed description of the algorithm.
66*> \endverbatim
67*
68* Arguments:
69* ==========
70*
71*> \param[in] UPLO
72*> \verbatim
73*> UPLO is CHARACTER*1
74*> = 'U': B is upper bidiagonal;
75*> = 'L': B is lower bidiagonal.
76*> \endverbatim
77*>
78*> \param[in] N
79*> \verbatim
80*> N is INTEGER
81*> The order of the matrix B. N >= 0.
82*> \endverbatim
83*>
84*> \param[in] NCVT
85*> \verbatim
86*> NCVT is INTEGER
87*> The number of columns of the matrix VT. NCVT >= 0.
88*> \endverbatim
89*>
90*> \param[in] NRU
91*> \verbatim
92*> NRU is INTEGER
93*> The number of rows of the matrix U. NRU >= 0.
94*> \endverbatim
95*>
96*> \param[in] NCC
97*> \verbatim
98*> NCC is INTEGER
99*> The number of columns of the matrix C. NCC >= 0.
100*> \endverbatim
101*>
102*> \param[in,out] D
103*> \verbatim
104*> D is DOUBLE PRECISION array, dimension (N)
105*> On entry, the n diagonal elements of the bidiagonal matrix B.
106*> On exit, if INFO=0, the singular values of B in decreasing
107*> order.
108*> \endverbatim
109*>
110*> \param[in,out] E
111*> \verbatim
112*> E is DOUBLE PRECISION array, dimension (N-1)
113*> On entry, the N-1 offdiagonal elements of the bidiagonal
114*> matrix B.
115*> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
116*> will contain the diagonal and superdiagonal elements of a
117*> bidiagonal matrix orthogonally equivalent to the one given
118*> as input.
119*> \endverbatim
120*>
121*> \param[in,out] VT
122*> \verbatim
123*> VT is COMPLEX*16 array, dimension (LDVT, NCVT)
124*> On entry, an N-by-NCVT matrix VT.
125*> On exit, VT is overwritten by P**H * VT.
126*> Not referenced if NCVT = 0.
127*> \endverbatim
128*>
129*> \param[in] LDVT
130*> \verbatim
131*> LDVT is INTEGER
132*> The leading dimension of the array VT.
133*> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
134*> \endverbatim
135*>
136*> \param[in,out] U
137*> \verbatim
138*> U is COMPLEX*16 array, dimension (LDU, N)
139*> On entry, an NRU-by-N matrix U.
140*> On exit, U is overwritten by U * Q.
141*> Not referenced if NRU = 0.
142*> \endverbatim
143*>
144*> \param[in] LDU
145*> \verbatim
146*> LDU is INTEGER
147*> The leading dimension of the array U. LDU >= max(1,NRU).
148*> \endverbatim
149*>
150*> \param[in,out] C
151*> \verbatim
152*> C is COMPLEX*16 array, dimension (LDC, NCC)
153*> On entry, an N-by-NCC matrix C.
154*> On exit, C is overwritten by Q**H * C.
155*> Not referenced if NCC = 0.
156*> \endverbatim
157*>
158*> \param[in] LDC
159*> \verbatim
160*> LDC is INTEGER
161*> The leading dimension of the array C.
162*> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
163*> \endverbatim
164*>
165*> \param[out] RWORK
166*> \verbatim
167*> RWORK is DOUBLE PRECISION array, dimension (LRWORK)
168*> LRWORK = 4*N, if NCVT = NRU = NCC = 0, and
169*> LRWORK = 4*(N-1), otherwise
170*> \endverbatim
171*>
172*> \param[out] INFO
173*> \verbatim
174*> INFO is INTEGER
175*> = 0: successful exit
176*> < 0: If INFO = -i, the i-th argument had an illegal value
177*> > 0: the algorithm did not converge; D and E contain the
178*> elements of a bidiagonal matrix which is orthogonally
179*> similar to the input matrix B; if INFO = i, i
180*> elements of E have not converged to zero.
181*> \endverbatim
182*
183*> \par Internal Parameters:
184* =========================
185*>
186*> \verbatim
187*> TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
188*> TOLMUL controls the convergence criterion of the QR loop.
189*> If it is positive, TOLMUL*EPS is the desired relative
190*> precision in the computed singular values.
191*> If it is negative, abs(TOLMUL*EPS*sigma_max) is the
192*> desired absolute accuracy in the computed singular
193*> values (corresponds to relative accuracy
194*> abs(TOLMUL*EPS) in the largest singular value.
195*> abs(TOLMUL) should be between 1 and 1/EPS, and preferably
196*> between 10 (for fast convergence) and .1/EPS
197*> (for there to be some accuracy in the results).
198*> Default is to lose at either one eighth or 2 of the
199*> available decimal digits in each computed singular value
200*> (whichever is smaller).
201*>
202*> MAXITR INTEGER, default = 6
203*> MAXITR controls the maximum number of passes of the
204*> algorithm through its inner loop. The algorithms stops
205*> (and so fails to converge) if the number of passes
206*> through the inner loop exceeds MAXITR*N**2.
207*>
208*> \endverbatim
209*
210*> \par Note:
211* ===========
212*>
213*> \verbatim
214*> Bug report from Cezary Dendek.
215*> On November 3rd 2023, the INTEGER variable MAXIT = MAXITR*N**2 is
216*> removed since it can overflow pretty easily (for N larger or equal
217*> than 18,919). We instead use MAXITDIVN = MAXITR*N.
218*> \endverbatim
219*
220* Authors:
221* ========
222*
223*> \author Univ. of Tennessee
224*> \author Univ. of California Berkeley
225*> \author Univ. of Colorado Denver
226*> \author NAG Ltd.
227*
228*> \ingroup bdsqr
229*
230* =====================================================================
231 SUBROUTINE zbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
232 $ LDU, C, LDC, RWORK, INFO )
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,
283 $ zdrot,
284 $ zdscal, zlasr, zswap
285* ..
286* .. Intrinsic Functions ..
287 INTRINSIC abs, dble, max, min, sign, sqrt
288* ..
289* .. Executable Statements ..
290*
291* Test the input parameters.
292*
293 info = 0
294 lower = lsame( uplo, 'L' )
295 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
296 info = -1
297 ELSE IF( n.LT.0 ) THEN
298 info = -2
299 ELSE IF( ncvt.LT.0 ) THEN
300 info = -3
301 ELSE IF( nru.LT.0 ) THEN
302 info = -4
303 ELSE IF( ncc.LT.0 ) THEN
304 info = -5
305 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
306 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
307 info = -9
308 ELSE IF( ldu.LT.max( 1, nru ) ) THEN
309 info = -11
310 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
311 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
312 info = -13
313 END IF
314 IF( info.NE.0 ) THEN
315 CALL xerbla( 'ZBDSQR', -info )
316 RETURN
317 END IF
318 IF( n.EQ.0 )
319 $ RETURN
320 IF( n.EQ.1 )
321 $ GO TO 160
322*
323* ROTATE is true if any singular vectors desired, false otherwise
324*
325 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
326*
327* If no singular vectors desired, use qd algorithm
328*
329 IF( .NOT.rotate ) THEN
330 CALL dlasq1( n, d, e, rwork, info )
331*
332* If INFO equals 2, dqds didn't finish, try to finish
333*
334 IF( info .NE. 2 ) RETURN
335 info = 0
336 END IF
337*
338 nm1 = n - 1
339 nm12 = nm1 + nm1
340 nm13 = nm12 + nm1
341 idir = 0
342*
343* Get machine constants
344*
345 eps = dlamch( 'Epsilon' )
346 unfl = dlamch( 'Safe minimum' )
347*
348* If matrix lower bidiagonal, rotate to be upper bidiagonal
349* by applying Givens rotations on the left
350*
351 IF( lower ) THEN
352 DO 10 i = 1, n - 1
353 CALL dlartg( d( i ), e( i ), cs, sn, r )
354 d( i ) = r
355 e( i ) = sn*d( i+1 )
356 d( i+1 ) = cs*d( i+1 )
357 rwork( i ) = cs
358 rwork( nm1+i ) = sn
359 10 CONTINUE
360*
361* Update singular vectors if desired
362*
363 IF( nru.GT.0 )
364 $ CALL zlasr( 'R', 'V', 'F', nru, n, rwork( 1 ),
365 $ rwork( n ),
366 $ u, ldu )
367 IF( ncc.GT.0 )
368 $ CALL zlasr( 'L', 'V', 'F', n, ncc, rwork( 1 ),
369 $ rwork( n ),
370 $ c, ldc )
371 END IF
372*
373* Compute singular values to relative accuracy TOL
374* (By setting TOL to be negative, algorithm will compute
375* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
376*
377 tolmul = max( ten, min( hndrd, eps**meigth ) )
378 tol = tolmul*eps
379*
380* Compute approximate maximum, minimum singular values
381*
382 smax = zero
383 DO 20 i = 1, n
384 smax = max( smax, abs( d( i ) ) )
385 20 CONTINUE
386 DO 30 i = 1, n - 1
387 smax = max( smax, abs( e( i ) ) )
388 30 CONTINUE
389 smin = zero
390 IF( tol.GE.zero ) THEN
391*
392* Relative accuracy desired
393*
394 sminoa = abs( d( 1 ) )
395 IF( sminoa.EQ.zero )
396 $ GO TO 50
397 mu = sminoa
398 DO 40 i = 2, n
399 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
400 sminoa = min( sminoa, mu )
401 IF( sminoa.EQ.zero )
402 $ GO TO 50
403 40 CONTINUE
404 50 CONTINUE
405 sminoa = sminoa / sqrt( dble( n ) )
406 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
407 ELSE
408*
409* Absolute accuracy desired
410*
411 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
412 END IF
413*
414* Prepare for main iteration loop for the singular values
415* (MAXIT is the maximum number of passes through the inner
416* loop permitted before nonconvergence signalled.)
417*
418 maxitdivn = maxitr*n
419 iterdivn = 0
420 iter = -1
421 oldll = -1
422 oldm = -1
423*
424* M points to last element of unconverged part of matrix
425*
426 m = n
427*
428* Begin main iteration loop
429*
430 60 CONTINUE
431*
432* Check for convergence or exceeding iteration count
433*
434 IF( m.LE.1 )
435 $ GO TO 160
436 IF( iter.GE.n ) THEN
437 iter = iter - n
438 iterdivn = iterdivn + 1
439 IF( iterdivn.GE.maxitdivn )
440 $ GO TO 200
441 END IF
442*
443* Find diagonal block of matrix to work on
444*
445 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
446 $ d( m ) = zero
447 smax = abs( d( m ) )
448 DO 70 lll = 1, m - 1
449 ll = m - lll
450 abss = abs( d( ll ) )
451 abse = abs( e( ll ) )
452 IF( tol.LT.zero .AND. abss.LE.thresh )
453 $ d( ll ) = zero
454 IF( abse.LE.thresh )
455 $ GO TO 80
456 smax = max( smax, abss, abse )
457 70 CONTINUE
458 ll = 0
459 GO TO 90
460 80 CONTINUE
461 e( ll ) = zero
462*
463* Matrix splits since E(LL) = 0
464*
465 IF( ll.EQ.m-1 ) THEN
466*
467* Convergence of bottom singular value, return to top of loop
468*
469 m = m - 1
470 GO TO 60
471 END IF
472 90 CONTINUE
473 ll = ll + 1
474*
475* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
476*
477 IF( ll.EQ.m-1 ) THEN
478*
479* 2 by 2 block, handle separately
480*
481 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
482 $ cosr, sinl, cosl )
483 d( m-1 ) = sigmx
484 e( m-1 ) = zero
485 d( m ) = sigmn
486*
487* Compute singular vectors, if desired
488*
489 IF( ncvt.GT.0 )
490 $ CALL zdrot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
491 $ cosr, sinr )
492 IF( nru.GT.0 )
493 $ CALL zdrot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
494 $ sinl )
495 IF( ncc.GT.0 )
496 $ CALL zdrot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
497 $ sinl )
498 m = m - 2
499 GO TO 60
500 END IF
501*
502* If working on new submatrix, choose shift direction
503* (from larger end diagonal element towards smaller)
504*
505 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
506 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
507*
508* Chase bulge from top (big end) to bottom (small end)
509*
510 idir = 1
511 ELSE
512*
513* Chase bulge from bottom (big end) to top (small end)
514*
515 idir = 2
516 END IF
517 END IF
518*
519* Apply convergence tests
520*
521 IF( idir.EQ.1 ) THEN
522*
523* Run convergence test in forward direction
524* First apply standard test to bottom of matrix
525*
526 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
527 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
528 e( m-1 ) = zero
529 GO TO 60
530 END IF
531*
532 IF( tol.GE.zero ) THEN
533*
534* If relative accuracy desired,
535* apply convergence criterion forward
536*
537 mu = abs( d( ll ) )
538 smin = mu
539 DO 100 lll = ll, m - 1
540 IF( abs( e( lll ) ).LE.tol*mu ) THEN
541 e( lll ) = zero
542 GO TO 60
543 END IF
544 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
545 smin = min( smin, mu )
546 100 CONTINUE
547 END IF
548*
549 ELSE
550*
551* Run convergence test in backward direction
552* First apply standard test to top of matrix
553*
554 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
555 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
556 e( ll ) = zero
557 GO TO 60
558 END IF
559*
560 IF( tol.GE.zero ) THEN
561*
562* If relative accuracy desired,
563* apply convergence criterion backward
564*
565 mu = abs( d( m ) )
566 smin = mu
567 DO 110 lll = m - 1, ll, -1
568 IF( abs( e( lll ) ).LE.tol*mu ) THEN
569 e( lll ) = zero
570 GO TO 60
571 END IF
572 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
573 smin = min( smin, mu )
574 110 CONTINUE
575 END IF
576 END IF
577 oldll = ll
578 oldm = m
579*
580* Compute shift. First, test if shifting would ruin relative
581* accuracy, and if so set the shift to zero.
582*
583 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
584 $ max( eps, hndrth*tol ) ) THEN
585*
586* Use a zero shift to avoid loss of relative accuracy
587*
588 shift = zero
589 ELSE
590*
591* Compute the shift from 2-by-2 block at end of matrix
592*
593 IF( idir.EQ.1 ) THEN
594 sll = abs( d( ll ) )
595 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
596 ELSE
597 sll = abs( d( m ) )
598 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
599 END IF
600*
601* Test if shift negligible, and if so set to zero
602*
603 IF( sll.GT.zero ) THEN
604 IF( ( shift / sll )**2.LT.eps )
605 $ shift = zero
606 END IF
607 END IF
608*
609* Increment iteration count
610*
611 iter = iter + m - ll
612*
613* If SHIFT = 0, do simplified QR iteration
614*
615 IF( shift.EQ.zero ) THEN
616 IF( idir.EQ.1 ) THEN
617*
618* Chase bulge from top to bottom
619* Save cosines and sines for later singular vector updates
620*
621 cs = one
622 oldcs = one
623 DO 120 i = ll, m - 1
624 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
625 IF( i.GT.ll )
626 $ e( i-1 ) = oldsn*r
627 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn,
628 $ d( i ) )
629 rwork( i-ll+1 ) = cs
630 rwork( i-ll+1+nm1 ) = sn
631 rwork( i-ll+1+nm12 ) = oldcs
632 rwork( i-ll+1+nm13 ) = oldsn
633 120 CONTINUE
634 h = d( m )*cs
635 d( m ) = h*oldcs
636 e( m-1 ) = h*oldsn
637*
638* Update singular vectors
639*
640 IF( ncvt.GT.0 )
641 $ CALL zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
642 $ rwork( n ), vt( ll, 1 ), ldvt )
643 IF( nru.GT.0 )
644 $ CALL zlasr( 'R', 'V', 'F', nru, m-ll+1,
645 $ rwork( nm12+1 ),
646 $ rwork( nm13+1 ), u( 1, ll ), ldu )
647 IF( ncc.GT.0 )
648 $ CALL zlasr( 'L', 'V', 'F', m-ll+1, ncc,
649 $ rwork( nm12+1 ),
650 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
651*
652* Test convergence
653*
654 IF( abs( e( m-1 ) ).LE.thresh )
655 $ e( m-1 ) = zero
656*
657 ELSE
658*
659* Chase bulge from bottom to top
660* Save cosines and sines for later singular vector updates
661*
662 cs = one
663 oldcs = one
664 DO 130 i = m, ll + 1, -1
665 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
666 IF( i.LT.m )
667 $ e( i ) = oldsn*r
668 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn,
669 $ d( i ) )
670 rwork( i-ll ) = cs
671 rwork( i-ll+nm1 ) = -sn
672 rwork( i-ll+nm12 ) = oldcs
673 rwork( i-ll+nm13 ) = -oldsn
674 130 CONTINUE
675 h = d( ll )*cs
676 d( ll ) = h*oldcs
677 e( ll ) = h*oldsn
678*
679* Update singular vectors
680*
681 IF( ncvt.GT.0 )
682 $ CALL zlasr( 'L', 'V', 'B', m-ll+1, ncvt,
683 $ rwork( nm12+1 ),
684 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
685 IF( nru.GT.0 )
686 $ CALL zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
687 $ rwork( n ), u( 1, ll ), ldu )
688 IF( ncc.GT.0 )
689 $ CALL zlasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
690 $ rwork( n ), c( ll, 1 ), ldc )
691*
692* Test convergence
693*
694 IF( abs( e( ll ) ).LE.thresh )
695 $ e( ll ) = zero
696 END IF
697 ELSE
698*
699* Use nonzero shift
700*
701 IF( idir.EQ.1 ) THEN
702*
703* Chase bulge from top to bottom
704* Save cosines and sines for later singular vector updates
705*
706 f = ( abs( d( ll ) )-shift )*
707 $ ( sign( one, d( ll ) )+shift / d( ll ) )
708 g = e( ll )
709 DO 140 i = ll, m - 1
710 CALL dlartg( f, g, cosr, sinr, r )
711 IF( i.GT.ll )
712 $ e( i-1 ) = r
713 f = cosr*d( i ) + sinr*e( i )
714 e( i ) = cosr*e( i ) - sinr*d( i )
715 g = sinr*d( i+1 )
716 d( i+1 ) = cosr*d( i+1 )
717 CALL dlartg( f, g, cosl, sinl, r )
718 d( i ) = r
719 f = cosl*e( i ) + sinl*d( i+1 )
720 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
721 IF( i.LT.m-1 ) THEN
722 g = sinl*e( i+1 )
723 e( i+1 ) = cosl*e( i+1 )
724 END IF
725 rwork( i-ll+1 ) = cosr
726 rwork( i-ll+1+nm1 ) = sinr
727 rwork( i-ll+1+nm12 ) = cosl
728 rwork( i-ll+1+nm13 ) = sinl
729 140 CONTINUE
730 e( m-1 ) = f
731*
732* Update singular vectors
733*
734 IF( ncvt.GT.0 )
735 $ CALL zlasr( 'L', 'V', 'F', m-ll+1, ncvt, rwork( 1 ),
736 $ rwork( n ), vt( ll, 1 ), ldvt )
737 IF( nru.GT.0 )
738 $ CALL zlasr( 'R', 'V', 'F', nru, m-ll+1,
739 $ rwork( nm12+1 ),
740 $ rwork( nm13+1 ), u( 1, ll ), ldu )
741 IF( ncc.GT.0 )
742 $ CALL zlasr( 'L', 'V', 'F', m-ll+1, ncc,
743 $ rwork( nm12+1 ),
744 $ rwork( nm13+1 ), c( ll, 1 ), ldc )
745*
746* Test convergence
747*
748 IF( abs( e( m-1 ) ).LE.thresh )
749 $ e( m-1 ) = zero
750*
751 ELSE
752*
753* Chase bulge from bottom to top
754* Save cosines and sines for later singular vector updates
755*
756 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
757 $ d( m ) )
758 g = e( m-1 )
759 DO 150 i = m, ll + 1, -1
760 CALL dlartg( f, g, cosr, sinr, r )
761 IF( i.LT.m )
762 $ e( i ) = r
763 f = cosr*d( i ) + sinr*e( i-1 )
764 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
765 g = sinr*d( i-1 )
766 d( i-1 ) = cosr*d( i-1 )
767 CALL dlartg( f, g, cosl, sinl, r )
768 d( i ) = r
769 f = cosl*e( i-1 ) + sinl*d( i-1 )
770 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
771 IF( i.GT.ll+1 ) THEN
772 g = sinl*e( i-2 )
773 e( i-2 ) = cosl*e( i-2 )
774 END IF
775 rwork( i-ll ) = cosr
776 rwork( i-ll+nm1 ) = -sinr
777 rwork( i-ll+nm12 ) = cosl
778 rwork( i-ll+nm13 ) = -sinl
779 150 CONTINUE
780 e( ll ) = f
781*
782* Test convergence
783*
784 IF( abs( e( ll ) ).LE.thresh )
785 $ e( ll ) = zero
786*
787* Update singular vectors if desired
788*
789 IF( ncvt.GT.0 )
790 $ CALL zlasr( 'L', 'V', 'B', m-ll+1, ncvt,
791 $ rwork( nm12+1 ),
792 $ rwork( nm13+1 ), vt( ll, 1 ), ldvt )
793 IF( nru.GT.0 )
794 $ CALL zlasr( 'R', 'V', 'B', nru, m-ll+1, rwork( 1 ),
795 $ rwork( n ), u( 1, ll ), ldu )
796 IF( ncc.GT.0 )
797 $ CALL zlasr( 'L', 'V', 'B', m-ll+1, ncc, rwork( 1 ),
798 $ rwork( n ), c( ll, 1 ), ldc )
799 END IF
800 END IF
801*
802* QR iteration finished, go back and check convergence
803*
804 GO TO 60
805*
806* All singular values converged, so make them positive
807*
808 160 CONTINUE
809 DO 170 i = 1, n
810 IF( d( i ).EQ.zero ) THEN
811*
812* Avoid -ZERO
813*
814 d( i ) = zero
815 END IF
816 IF( d( i ).LT.zero ) THEN
817 d( i ) = -d( i )
818*
819* Change sign of singular vectors, if desired
820*
821 IF( ncvt.GT.0 )
822 $ CALL zdscal( ncvt, negone, vt( i, 1 ), ldvt )
823 END IF
824 170 CONTINUE
825*
826* Sort the singular values into decreasing order (insertion sort on
827* singular values, but only one transposition per singular vector)
828*
829 DO 190 i = 1, n - 1
830*
831* Scan for smallest D(I)
832*
833 isub = 1
834 smin = d( 1 )
835 DO 180 j = 2, n + 1 - i
836 IF( d( j ).LE.smin ) THEN
837 isub = j
838 smin = d( j )
839 END IF
840 180 CONTINUE
841 IF( isub.NE.n+1-i ) THEN
842*
843* Swap singular values and vectors
844*
845 d( isub ) = d( n+1-i )
846 d( n+1-i ) = smin
847 IF( ncvt.GT.0 )
848 $ CALL zswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
849 $ ldvt )
850 IF( nru.GT.0 )
851 $ CALL zswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
852 IF( ncc.GT.0 )
853 $ CALL zswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ),
854 $ ldc )
855 END IF
856 190 CONTINUE
857 GO TO 220
858*
859* Maximum number of iterations exceeded, failure to converge
860*
861 200 CONTINUE
862 info = 0
863 DO 210 i = 1, n - 1
864 IF( e( i ).NE.zero )
865 $ info = info + 1
866 210 CONTINUE
867 220 CONTINUE
868 RETURN
869*
870* End of ZBDSQR
871*
872 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
ZBDSQR
Definition zbdsqr.f:233
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:103
subroutine dlasq1(n, d, e, work, info)
DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition dlasq1.f:106
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:198
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:134
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