LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dbdsqr.f
Go to the documentation of this file.
1*> \brief \b DBDSQR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DBDSQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
22* LDU, C, LDC, WORK, 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 C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30* $ VT( LDVT, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DBDSQR 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**T
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**T*VT instead of
51*> P**T, for given real input matrices U and VT. When U and VT are the
52*> orthogonal matrices that reduce a general matrix A to bidiagonal
53*> form: A = U*B*VT, as computed by DGEBRD, then
54*>
55*> A = (U*Q) * S * (P**T*VT)
56*>
57*> is the SVD of A. Optionally, the subroutine may also compute Q**T*C
58*> for a given real 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 DOUBLE PRECISION array, dimension (LDVT, NCVT)
126*> On entry, an N-by-NCVT matrix VT.
127*> On exit, VT is overwritten by P**T * 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDC, NCC)
155*> On entry, an N-by-NCC matrix C.
156*> On exit, C is overwritten by Q**T * 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] WORK
168*> \verbatim
169*> WORK is DOUBLE PRECISION array, dimension (4*(N-1))
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:
178*> if NCVT = NRU = NCC = 0,
179*> = 1, a split was marked by a positive value in E
180*> = 2, current block of Z not diagonalized after 30*N
181*> iterations (in inner while loop)
182*> = 3, termination criterion of outer while loop not met
183*> (program created more than N unreduced blocks)
184*> else NCVT = NRU = NCC = 0,
185*> the algorithm did not converge; D and E contain the
186*> elements of a bidiagonal matrix which is orthogonally
187*> similar to the input matrix B; if INFO = i, i
188*> elements of E have not converged to zero.
189*> \endverbatim
190*
191*> \par Internal Parameters:
192* =========================
193*>
194*> \verbatim
195*> TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
196*> TOLMUL controls the convergence criterion of the QR loop.
197*> If it is positive, TOLMUL*EPS is the desired relative
198*> precision in the computed singular values.
199*> If it is negative, abs(TOLMUL*EPS*sigma_max) is the
200*> desired absolute accuracy in the computed singular
201*> values (corresponds to relative accuracy
202*> abs(TOLMUL*EPS) in the largest singular value.
203*> abs(TOLMUL) should be between 1 and 1/EPS, and preferably
204*> between 10 (for fast convergence) and .1/EPS
205*> (for there to be some accuracy in the results).
206*> Default is to lose at either one eighth or 2 of the
207*> available decimal digits in each computed singular value
208*> (whichever is smaller).
209*>
210*> MAXITR INTEGER, default = 6
211*> MAXITR controls the maximum number of passes of the
212*> algorithm through its inner loop. The algorithms stops
213*> (and so fails to converge) if the number of passes
214*> through the inner loop exceeds MAXITR*N**2.
215*>
216*> \endverbatim
217*
218*> \par Note:
219* ===========
220*>
221*> \verbatim
222*> Bug report from Cezary Dendek.
223*> On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
224*> removed since it can overflow pretty easily (for N larger or equal
225*> than 18,919). We instead use MAXITDIVN = MAXITR*N.
226*> \endverbatim
227*
228* Authors:
229* ========
230*
231*> \author Univ. of Tennessee
232*> \author Univ. of California Berkeley
233*> \author Univ. of Colorado Denver
234*> \author NAG Ltd.
235*
236*> \ingroup bdsqr
237*
238* =====================================================================
239 SUBROUTINE dbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
240 $ LDU, C, LDC, WORK, INFO )
241*
242* -- LAPACK computational routine --
243* -- LAPACK is a software package provided by Univ. of Tennessee, --
244* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
245*
246* .. Scalar Arguments ..
247 CHARACTER UPLO
248 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
249* ..
250* .. Array Arguments ..
251 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
252 $ vt( ldvt, * ), work( * )
253* ..
254*
255* =====================================================================
256*
257* .. Parameters ..
258 DOUBLE PRECISION ZERO
259 parameter( zero = 0.0d0 )
260 DOUBLE PRECISION ONE
261 parameter( one = 1.0d0 )
262 DOUBLE PRECISION NEGONE
263 parameter( negone = -1.0d0 )
264 DOUBLE PRECISION HNDRTH
265 parameter( hndrth = 0.01d0 )
266 DOUBLE PRECISION TEN
267 parameter( ten = 10.0d0 )
268 DOUBLE PRECISION HNDRD
269 parameter( hndrd = 100.0d0 )
270 DOUBLE PRECISION MEIGTH
271 parameter( meigth = -0.125d0 )
272 INTEGER MAXITR
273 parameter( maxitr = 6 )
274* ..
275* .. Local Scalars ..
276 LOGICAL LOWER, ROTATE
277 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
278 $ maxitdivn, nm1, nm12, nm13, oldll, oldm
279 DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
280 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
281 $ sinr, sll, smax, smin, sminoa,
282 $ sn, thresh, tol, tolmul, unfl
283* ..
284* .. External Functions ..
285 LOGICAL LSAME
286 DOUBLE PRECISION DLAMCH
287 EXTERNAL lsame, dlamch
288* ..
289* .. External Subroutines ..
290 EXTERNAL dlartg, dlas2, dlasq1, dlasr, dlasv2, drot,
291 $ dscal, dswap, xerbla
292* ..
293* .. Intrinsic Functions ..
294 INTRINSIC abs, dble, max, min, sign, sqrt
295* ..
296* .. Executable Statements ..
297*
298* Test the input parameters.
299*
300 info = 0
301 lower = lsame( uplo, 'L' )
302 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
303 info = -1
304 ELSE IF( n.LT.0 ) THEN
305 info = -2
306 ELSE IF( ncvt.LT.0 ) THEN
307 info = -3
308 ELSE IF( nru.LT.0 ) THEN
309 info = -4
310 ELSE IF( ncc.LT.0 ) THEN
311 info = -5
312 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
313 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
314 info = -9
315 ELSE IF( ldu.LT.max( 1, nru ) ) THEN
316 info = -11
317 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
318 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
319 info = -13
320 END IF
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'DBDSQR', -info )
323 RETURN
324 END IF
325 IF( n.EQ.0 )
326 $ RETURN
327 IF( n.EQ.1 )
328 $ GO TO 160
329*
330* ROTATE is true if any singular vectors desired, false otherwise
331*
332 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
333*
334* If no singular vectors desired, use qd algorithm
335*
336 IF( .NOT.rotate ) THEN
337 CALL dlasq1( n, d, e, work, info )
338*
339* If INFO equals 2, dqds didn't finish, try to finish
340*
341 IF( info .NE. 2 ) RETURN
342 info = 0
343 END IF
344*
345 nm1 = n - 1
346 nm12 = nm1 + nm1
347 nm13 = nm12 + nm1
348 idir = 0
349*
350* Get machine constants
351*
352 eps = dlamch( 'Epsilon' )
353 unfl = dlamch( 'Safe minimum' )
354*
355* If matrix lower bidiagonal, rotate to be upper bidiagonal
356* by applying Givens rotations on the left
357*
358 IF( lower ) THEN
359 DO 10 i = 1, n - 1
360 CALL dlartg( d( i ), e( i ), cs, sn, r )
361 d( i ) = r
362 e( i ) = sn*d( i+1 )
363 d( i+1 ) = cs*d( i+1 )
364 work( i ) = cs
365 work( nm1+i ) = sn
366 10 CONTINUE
367*
368* Update singular vectors if desired
369*
370 IF( nru.GT.0 )
371 $ CALL dlasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
372 $ ldu )
373 IF( ncc.GT.0 )
374 $ CALL dlasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), c,
375 $ ldc )
376 END IF
377*
378* Compute singular values to relative accuracy TOL
379* (By setting TOL to be negative, algorithm will compute
380* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
381*
382 tolmul = max( ten, min( hndrd, eps**meigth ) )
383 tol = tolmul*eps
384*
385* Compute approximate maximum, minimum singular values
386*
387 smax = zero
388 DO 20 i = 1, n
389 smax = max( smax, abs( d( i ) ) )
390 20 CONTINUE
391 DO 30 i = 1, n - 1
392 smax = max( smax, abs( e( i ) ) )
393 30 CONTINUE
394 smin = zero
395 IF( tol.GE.zero ) THEN
396*
397* Relative accuracy desired
398*
399 sminoa = abs( d( 1 ) )
400 IF( sminoa.EQ.zero )
401 $ GO TO 50
402 mu = sminoa
403 DO 40 i = 2, n
404 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
405 sminoa = min( sminoa, mu )
406 IF( sminoa.EQ.zero )
407 $ GO TO 50
408 40 CONTINUE
409 50 CONTINUE
410 sminoa = sminoa / sqrt( dble( n ) )
411 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
412 ELSE
413*
414* Absolute accuracy desired
415*
416 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
417 END IF
418*
419* Prepare for main iteration loop for the singular values
420* (MAXIT is the maximum number of passes through the inner
421* loop permitted before nonconvergence signalled.)
422*
423 maxitdivn = maxitr*n
424 iterdivn = 0
425 iter = -1
426 oldll = -1
427 oldm = -1
428*
429* M points to last element of unconverged part of matrix
430*
431 m = n
432*
433* Begin main iteration loop
434*
435 60 CONTINUE
436*
437* Check for convergence or exceeding iteration count
438*
439 IF( m.LE.1 )
440 $ GO TO 160
441*
442 IF( iter.GE.n ) THEN
443 iter = iter - n
444 iterdivn = iterdivn + 1
445 IF( iterdivn.GE.maxitdivn )
446 $ GO TO 200
447 END IF
448*
449* Find diagonal block of matrix to work on
450*
451 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
452 $ d( m ) = zero
453 smax = abs( d( m ) )
454 DO 70 lll = 1, m - 1
455 ll = m - lll
456 abss = abs( d( ll ) )
457 abse = abs( e( ll ) )
458 IF( tol.LT.zero .AND. abss.LE.thresh )
459 $ d( ll ) = zero
460 IF( abse.LE.thresh )
461 $ GO TO 80
462 smax = max( smax, abss, abse )
463 70 CONTINUE
464 ll = 0
465 GO TO 90
466 80 CONTINUE
467 e( ll ) = zero
468*
469* Matrix splits since E(LL) = 0
470*
471 IF( ll.EQ.m-1 ) THEN
472*
473* Convergence of bottom singular value, return to top of loop
474*
475 m = m - 1
476 GO TO 60
477 END IF
478 90 CONTINUE
479 ll = ll + 1
480*
481* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
482*
483 IF( ll.EQ.m-1 ) THEN
484*
485* 2 by 2 block, handle separately
486*
487 CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
488 $ cosr, sinl, cosl )
489 d( m-1 ) = sigmx
490 e( m-1 ) = zero
491 d( m ) = sigmn
492*
493* Compute singular vectors, if desired
494*
495 IF( ncvt.GT.0 )
496 $ CALL drot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
497 $ sinr )
498 IF( nru.GT.0 )
499 $ CALL drot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
500 IF( ncc.GT.0 )
501 $ CALL drot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
502 $ sinl )
503 m = m - 2
504 GO TO 60
505 END IF
506*
507* If working on new submatrix, choose shift direction
508* (from larger end diagonal element towards smaller)
509*
510 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
511 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
512*
513* Chase bulge from top (big end) to bottom (small end)
514*
515 idir = 1
516 ELSE
517*
518* Chase bulge from bottom (big end) to top (small end)
519*
520 idir = 2
521 END IF
522 END IF
523*
524* Apply convergence tests
525*
526 IF( idir.EQ.1 ) THEN
527*
528* Run convergence test in forward direction
529* First apply standard test to bottom of matrix
530*
531 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
532 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
533 e( m-1 ) = zero
534 GO TO 60
535 END IF
536*
537 IF( tol.GE.zero ) THEN
538*
539* If relative accuracy desired,
540* apply convergence criterion forward
541*
542 mu = abs( d( ll ) )
543 smin = mu
544 DO 100 lll = ll, m - 1
545 IF( abs( e( lll ) ).LE.tol*mu ) THEN
546 e( lll ) = zero
547 GO TO 60
548 END IF
549 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
550 smin = min( smin, mu )
551 100 CONTINUE
552 END IF
553*
554 ELSE
555*
556* Run convergence test in backward direction
557* First apply standard test to top of matrix
558*
559 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
560 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
561 e( ll ) = zero
562 GO TO 60
563 END IF
564*
565 IF( tol.GE.zero ) THEN
566*
567* If relative accuracy desired,
568* apply convergence criterion backward
569*
570 mu = abs( d( m ) )
571 smin = mu
572 DO 110 lll = m - 1, ll, -1
573 IF( abs( e( lll ) ).LE.tol*mu ) THEN
574 e( lll ) = zero
575 GO TO 60
576 END IF
577 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
578 smin = min( smin, mu )
579 110 CONTINUE
580 END IF
581 END IF
582 oldll = ll
583 oldm = m
584*
585* Compute shift. First, test if shifting would ruin relative
586* accuracy, and if so set the shift to zero.
587*
588 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
589 $ max( eps, hndrth*tol ) ) THEN
590*
591* Use a zero shift to avoid loss of relative accuracy
592*
593 shift = zero
594 ELSE
595*
596* Compute the shift from 2-by-2 block at end of matrix
597*
598 IF( idir.EQ.1 ) THEN
599 sll = abs( d( ll ) )
600 CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
601 ELSE
602 sll = abs( d( m ) )
603 CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
604 END IF
605*
606* Test if shift negligible, and if so set to zero
607*
608 IF( sll.GT.zero ) THEN
609 IF( ( shift / sll )**2.LT.eps )
610 $ shift = zero
611 END IF
612 END IF
613*
614* Increment iteration count
615*
616 iter = iter + m - ll
617*
618* If SHIFT = 0, do simplified QR iteration
619*
620 IF( shift.EQ.zero ) THEN
621 IF( idir.EQ.1 ) THEN
622*
623* Chase bulge from top to bottom
624* Save cosines and sines for later singular vector updates
625*
626 cs = one
627 oldcs = one
628 DO 120 i = ll, m - 1
629 CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
630 IF( i.GT.ll )
631 $ e( i-1 ) = oldsn*r
632 CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
633 work( i-ll+1 ) = cs
634 work( i-ll+1+nm1 ) = sn
635 work( i-ll+1+nm12 ) = oldcs
636 work( i-ll+1+nm13 ) = oldsn
637 120 CONTINUE
638 h = d( m )*cs
639 d( m ) = h*oldcs
640 e( m-1 ) = h*oldsn
641*
642* Update singular vectors
643*
644 IF( ncvt.GT.0 )
645 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
646 $ work( n ), vt( ll, 1 ), ldvt )
647 IF( nru.GT.0 )
648 $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
649 $ work( nm13+1 ), u( 1, ll ), ldu )
650 IF( ncc.GT.0 )
651 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
652 $ work( nm13+1 ), c( ll, 1 ), ldc )
653*
654* Test convergence
655*
656 IF( abs( e( m-1 ) ).LE.thresh )
657 $ e( m-1 ) = zero
658*
659 ELSE
660*
661* Chase bulge from bottom to top
662* Save cosines and sines for later singular vector updates
663*
664 cs = one
665 oldcs = one
666 DO 130 i = m, ll + 1, -1
667 CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
668 IF( i.LT.m )
669 $ e( i ) = oldsn*r
670 CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
671 work( i-ll ) = cs
672 work( i-ll+nm1 ) = -sn
673 work( i-ll+nm12 ) = oldcs
674 work( i-ll+nm13 ) = -oldsn
675 130 CONTINUE
676 h = d( ll )*cs
677 d( ll ) = h*oldcs
678 e( ll ) = h*oldsn
679*
680* Update singular vectors
681*
682 IF( ncvt.GT.0 )
683 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
684 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
685 IF( nru.GT.0 )
686 $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
687 $ work( n ), u( 1, ll ), ldu )
688 IF( ncc.GT.0 )
689 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
690 $ work( 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 work( i-ll+1 ) = cosr
726 work( i-ll+1+nm1 ) = sinr
727 work( i-ll+1+nm12 ) = cosl
728 work( 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 dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
736 $ work( n ), vt( ll, 1 ), ldvt )
737 IF( nru.GT.0 )
738 $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
739 $ work( nm13+1 ), u( 1, ll ), ldu )
740 IF( ncc.GT.0 )
741 $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
742 $ work( nm13+1 ), c( ll, 1 ), ldc )
743*
744* Test convergence
745*
746 IF( abs( e( m-1 ) ).LE.thresh )
747 $ e( m-1 ) = zero
748*
749 ELSE
750*
751* Chase bulge from bottom to top
752* Save cosines and sines for later singular vector updates
753*
754 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
755 $ d( m ) )
756 g = e( m-1 )
757 DO 150 i = m, ll + 1, -1
758 CALL dlartg( f, g, cosr, sinr, r )
759 IF( i.LT.m )
760 $ e( i ) = r
761 f = cosr*d( i ) + sinr*e( i-1 )
762 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
763 g = sinr*d( i-1 )
764 d( i-1 ) = cosr*d( i-1 )
765 CALL dlartg( f, g, cosl, sinl, r )
766 d( i ) = r
767 f = cosl*e( i-1 ) + sinl*d( i-1 )
768 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
769 IF( i.GT.ll+1 ) THEN
770 g = sinl*e( i-2 )
771 e( i-2 ) = cosl*e( i-2 )
772 END IF
773 work( i-ll ) = cosr
774 work( i-ll+nm1 ) = -sinr
775 work( i-ll+nm12 ) = cosl
776 work( i-ll+nm13 ) = -sinl
777 150 CONTINUE
778 e( ll ) = f
779*
780* Test convergence
781*
782 IF( abs( e( ll ) ).LE.thresh )
783 $ e( ll ) = zero
784*
785* Update singular vectors if desired
786*
787 IF( ncvt.GT.0 )
788 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
789 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
790 IF( nru.GT.0 )
791 $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
792 $ work( n ), u( 1, ll ), ldu )
793 IF( ncc.GT.0 )
794 $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
795 $ work( n ), c( ll, 1 ), ldc )
796 END IF
797 END IF
798*
799* QR iteration finished, go back and check convergence
800*
801 GO TO 60
802*
803* All singular values converged, so make them positive
804*
805 160 CONTINUE
806 DO 170 i = 1, n
807 IF( d( i ).LT.zero ) THEN
808 d( i ) = -d( i )
809*
810* Change sign of singular vectors, if desired
811*
812 IF( ncvt.GT.0 )
813 $ CALL dscal( ncvt, negone, vt( i, 1 ), ldvt )
814 END IF
815 170 CONTINUE
816*
817* Sort the singular values into decreasing order (insertion sort on
818* singular values, but only one transposition per singular vector)
819*
820 DO 190 i = 1, n - 1
821*
822* Scan for smallest D(I)
823*
824 isub = 1
825 smin = d( 1 )
826 DO 180 j = 2, n + 1 - i
827 IF( d( j ).LE.smin ) THEN
828 isub = j
829 smin = d( j )
830 END IF
831 180 CONTINUE
832 IF( isub.NE.n+1-i ) THEN
833*
834* Swap singular values and vectors
835*
836 d( isub ) = d( n+1-i )
837 d( n+1-i ) = smin
838 IF( ncvt.GT.0 )
839 $ CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
840 $ ldvt )
841 IF( nru.GT.0 )
842 $ CALL dswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
843 IF( ncc.GT.0 )
844 $ CALL dswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
845 END IF
846 190 CONTINUE
847 GO TO 220
848*
849* Maximum number of iterations exceeded, failure to converge
850*
851 200 CONTINUE
852 info = 0
853 DO 210 i = 1, n - 1
854 IF( e( i ).NE.zero )
855 $ info = info + 1
856 210 CONTINUE
857 220 CONTINUE
858 RETURN
859*
860* End of DBDSQR
861*
862 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dbdsqr
subroutine dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DBDSQR
Definition dbdsqr.f:241
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
dlasr
subroutine dlasr(side, pivot, direct, m, n, c, s, a, lda)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition dlasr.f:199
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
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82