LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \htmlonly
9*> Download ZBDSQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zbdsqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zbdsqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zbdsqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
22* LDU, C, LDC, RWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER UPLO
26* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION D( * ), E( * ), RWORK( * )
30* COMPLEX*16 C( LDC, * ), U( LDU, * ), VT( LDVT, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> ZBDSQR computes the singular values and, optionally, the right and/or
40*> left singular vectors from the singular value decomposition (SVD) of
41*> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
42*> zero-shift QR algorithm. The SVD of B has the form
43*>
44*> B = Q * S * P**H
45*>
46*> where S is the diagonal matrix of singular values, Q is an orthogonal
47*> matrix of left singular vectors, and P is an orthogonal matrix of
48*> right singular vectors. If left singular vectors are requested, this
49*> subroutine actually returns U*Q instead of Q, and, if right singular
50*> vectors are requested, this subroutine returns P**H*VT instead of
51*> P**H, for given complex input matrices U and VT. When U and VT are
52*> the unitary matrices that reduce a general matrix A to bidiagonal
53*> form: A = U*B*VT, as computed by ZGEBRD, then
54*>
55*> A = (U*Q) * S * (P**H*VT)
56*>
57*> is the SVD of A. Optionally, the subroutine may also compute Q**H*C
58*> for a given complex input matrix C.
59*>
60*> See "Computing Small Singular Values of Bidiagonal Matrices With
61*> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
62*> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
63*> no. 5, pp. 873-912, Sept 1990) and
64*> "Accurate singular values and differential qd algorithms," by
65*> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
66*> Department, University of California at Berkeley, July 1992
67*> for a detailed description of the algorithm.
68*> \endverbatim
69*
70* Arguments:
71* ==========
72*
73*> \param[in] UPLO
74*> \verbatim
75*> UPLO is CHARACTER*1
76*> = 'U': B is upper bidiagonal;
77*> = 'L': B is lower bidiagonal.
78*> \endverbatim
79*>
80*> \param[in] N
81*> \verbatim
82*> N is INTEGER
83*> The order of the matrix B. N >= 0.
84*> \endverbatim
85*>
86*> \param[in] NCVT
87*> \verbatim
88*> NCVT is INTEGER
89*> The number of columns of the matrix VT. NCVT >= 0.
90*> \endverbatim
91*>
92*> \param[in] NRU
93*> \verbatim
94*> NRU is INTEGER
95*> The number of rows of the matrix U. NRU >= 0.
96*> \endverbatim
97*>
98*> \param[in] NCC
99*> \verbatim
100*> NCC is INTEGER
101*> The number of columns of the matrix C. NCC >= 0.
102*> \endverbatim
103*>
104*> \param[in,out] D
105*> \verbatim
106*> D is DOUBLE PRECISION array, dimension (N)
107*> On entry, the n diagonal elements of the bidiagonal matrix B.
108*> On exit, if INFO=0, the singular values of B in decreasing
109*> order.
110*> \endverbatim
111*>
112*> \param[in,out] E
113*> \verbatim
114*> E is DOUBLE PRECISION array, dimension (N-1)
115*> On entry, the N-1 offdiagonal elements of the bidiagonal
116*> matrix B.
117*> On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
118*> will contain the diagonal and superdiagonal elements of a
119*> bidiagonal matrix orthogonally equivalent to the one given
120*> as input.
121*> \endverbatim
122*>
123*> \param[in,out] VT
124*> \verbatim
125*> VT is COMPLEX*16 array, dimension (LDVT, NCVT)
126*> On entry, an N-by-NCVT matrix VT.
127*> On exit, VT is overwritten by P**H * VT.
128*> Not referenced if NCVT = 0.
129*> \endverbatim
130*>
131*> \param[in] LDVT
132*> \verbatim
133*> LDVT is INTEGER
134*> The leading dimension of the array VT.
135*> LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
136*> \endverbatim
137*>
138*> \param[in,out] U
139*> \verbatim
140*> U is COMPLEX*16 array, dimension (LDU, N)
141*> On entry, an NRU-by-N matrix U.
142*> On exit, U is overwritten by U * Q.
143*> Not referenced if NRU = 0.
144*> \endverbatim
145*>
146*> \param[in] LDU
147*> \verbatim
148*> LDU is INTEGER
149*> The leading dimension of the array U. LDU >= max(1,NRU).
150*> \endverbatim
151*>
152*> \param[in,out] C
153*> \verbatim
154*> C is COMPLEX*16 array, dimension (LDC, NCC)
155*> On entry, an N-by-NCC matrix C.
156*> On exit, C is overwritten by Q**H * C.
157*> Not referenced if NCC = 0.
158*> \endverbatim
159*>
160*> \param[in] LDC
161*> \verbatim
162*> LDC is INTEGER
163*> The leading dimension of the array C.
164*> LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
165*> \endverbatim
166*>
167*> \param[out] RWORK
168*> \verbatim
169*> RWORK is DOUBLE PRECISION array, dimension (4*N)
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, 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*
464* Convergence of bottom singular value, return to top of loop
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*
853 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zbdsqr
subroutine zbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
ZBDSQR
Definition zbdsqr.f:233
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlas2
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:105
dlasq1
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
zlasr
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
dlasv2
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
zdrot
subroutine zdrot(n, zx, incx, zy, incy, c, s)
ZDROT
Definition zdrot.f:98
zdscal
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81