LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sbdsqr.f
Go to the documentation of this file.
1*> \brief \b SBDSQR
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SBDSQR + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sbdsqr.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sbdsqr.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sbdsqr.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
20* LDU, C, LDC, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER UPLO
24* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
25* ..
26* .. Array Arguments ..
27* REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
28* $ VT( LDVT, * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SBDSQR 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**T
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**T*VT instead of
49*> P**T, for given real input matrices U and VT. When U and VT are the
50*> orthogonal matrices that reduce a general matrix A to bidiagonal
51*> form: A = U*B*VT, as computed by SGEBRD, then
52*>
53*> A = (U*Q) * S * (P**T*VT)
54*>
55*> is the SVD of A. Optionally, the subroutine may also compute Q**T*C
56*> for a given real 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 REAL 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 REAL 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 REAL array, dimension (LDVT, NCVT)
124*> On entry, an N-by-NCVT matrix VT.
125*> On exit, VT is overwritten by P**T * 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 REAL 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 REAL array, dimension (LDC, NCC)
153*> On entry, an N-by-NCC matrix C.
154*> On exit, C is overwritten by Q**T * 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] WORK
166*> \verbatim
167*> WORK is REAL array, dimension (LWORK)
168*> LWORK = 4*N, if NCVT = NRU = NCC = 0, and
169*> LWORK = 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:
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 REAL, 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*> \endverbatim
216*
217*> \par Note:
218* ===========
219*>
220*> \verbatim
221*> Bug report from Cezary Dendek.
222*> On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
223*> removed since it can overflow pretty easily (for N larger or equal
224*> than 18,919). We instead use MAXITDIVN = MAXITR*N.
225*> \endverbatim
226*
227* Authors:
228* ========
229*
230*> \author Univ. of Tennessee
231*> \author Univ. of California Berkeley
232*> \author Univ. of Colorado Denver
233*> \author NAG Ltd.
234*
235*> \ingroup bdsqr
236*
237* =====================================================================
238 SUBROUTINE sbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
239 $ LDU, C, LDC, WORK, INFO )
240*
241* -- LAPACK computational routine --
242* -- LAPACK is a software package provided by Univ. of Tennessee, --
243* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244*
245* .. Scalar Arguments ..
246 CHARACTER UPLO
247 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
248* ..
249* .. Array Arguments ..
250 REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
251 $ vt( ldvt, * ), work( * )
252* ..
253*
254* =====================================================================
255*
256* .. Parameters ..
257 REAL ZERO
258 parameter( zero = 0.0e0 )
259 REAL ONE
260 parameter( one = 1.0e0 )
261 REAL NEGONE
262 parameter( negone = -1.0e0 )
263 REAL HNDRTH
264 parameter( hndrth = 0.01e0 )
265 REAL TEN
266 parameter( ten = 10.0e0 )
267 REAL HNDRD
268 parameter( hndrd = 100.0e0 )
269 REAL MEIGTH
270 parameter( meigth = -0.125e0 )
271 INTEGER MAXITR
272 parameter( maxitr = 6 )
273* ..
274* .. Local Scalars ..
275 LOGICAL LOWER, ROTATE
276 INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
277 $ maxitdivn, nm1, nm12, nm13, oldll, oldm
278 REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
279 $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
280 $ sinr, sll, smax, smin, sminoa,
281 $ sn, thresh, tol, tolmul, unfl
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 REAL SLAMCH
286 EXTERNAL lsame, slamch
287* ..
288* .. External Subroutines ..
289 EXTERNAL slartg, slas2, slasq1, slasr, slasv2,
290 $ srot,
291 $ sscal, sswap, xerbla
292* ..
293* .. Intrinsic Functions ..
294 INTRINSIC abs, max, min, real, 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( 'SBDSQR', -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 slasq1( 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 = slamch( 'Epsilon' )
353 unfl = slamch( '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 slartg( 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 slasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ),
372 $ u,
373 $ ldu )
374 IF( ncc.GT.0 )
375 $ CALL slasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ),
376 $ c,
377 $ ldc )
378 END IF
379*
380* Compute singular values to relative accuracy TOL
381* (By setting TOL to be negative, algorithm will compute
382* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
383*
384 tolmul = max( ten, min( hndrd, eps**meigth ) )
385 tol = tolmul*eps
386*
387* Compute approximate maximum, minimum singular values
388*
389 smax = zero
390 DO 20 i = 1, n
391 smax = max( smax, abs( d( i ) ) )
392 20 CONTINUE
393 DO 30 i = 1, n - 1
394 smax = max( smax, abs( e( i ) ) )
395 30 CONTINUE
396 smin = zero
397 IF( tol.GE.zero ) THEN
398*
399* Relative accuracy desired
400*
401 sminoa = abs( d( 1 ) )
402 IF( sminoa.EQ.zero )
403 $ GO TO 50
404 mu = sminoa
405 DO 40 i = 2, n
406 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
407 sminoa = min( sminoa, mu )
408 IF( sminoa.EQ.zero )
409 $ GO TO 50
410 40 CONTINUE
411 50 CONTINUE
412 sminoa = sminoa / sqrt( real( n ) )
413 thresh = max( tol*sminoa, maxitr*(real( n )*(real( n )*unfl)) )
414 ELSE
415*
416* Absolute accuracy desired
417*
418 thresh = max( abs( tol )*smax, maxitr*(real( n )*
419 $ (real( n )*unfl)) )
420 END IF
421*
422* Prepare for main iteration loop for the singular values
423* (MAXIT is the maximum number of passes through the inner
424* loop permitted before nonconvergence signalled.)
425*
426 maxitdivn = maxitr*n
427 iterdivn = 0
428 iter = -1
429 oldll = -1
430 oldm = -1
431*
432* M points to last element of unconverged part of matrix
433*
434 m = n
435*
436* Begin main iteration loop
437*
438 60 CONTINUE
439*
440* Check for convergence or exceeding iteration count
441*
442 IF( m.LE.1 )
443 $ GO TO 160
444*
445 IF( iter.GE.n ) THEN
446 iter = iter - n
447 iterdivn = iterdivn + 1
448 IF( iterdivn.GE.maxitdivn )
449 $ GO TO 200
450 END IF
451*
452* Find diagonal block of matrix to work on
453*
454 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
455 $ d( m ) = zero
456 smax = abs( d( m ) )
457 DO 70 lll = 1, m - 1
458 ll = m - lll
459 abss = abs( d( ll ) )
460 abse = abs( e( ll ) )
461 IF( tol.LT.zero .AND. abss.LE.thresh )
462 $ d( ll ) = zero
463 IF( abse.LE.thresh )
464 $ GO TO 80
465 smax = max( smax, abss, abse )
466 70 CONTINUE
467 ll = 0
468 GO TO 90
469 80 CONTINUE
470 e( ll ) = zero
471*
472* Matrix splits since E(LL) = 0
473*
474 IF( ll.EQ.m-1 ) THEN
475*
476* Convergence of bottom singular value, return to top of loop
477*
478 m = m - 1
479 GO TO 60
480 END IF
481 90 CONTINUE
482 ll = ll + 1
483*
484* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
485*
486 IF( ll.EQ.m-1 ) THEN
487*
488* 2 by 2 block, handle separately
489*
490 CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
491 $ cosr, sinl, cosl )
492 d( m-1 ) = sigmx
493 e( m-1 ) = zero
494 d( m ) = sigmn
495*
496* Compute singular vectors, if desired
497*
498 IF( ncvt.GT.0 )
499 $ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt,
500 $ cosr,
501 $ sinr )
502 IF( nru.GT.0 )
503 $ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl,
504 $ sinl )
505 IF( ncc.GT.0 )
506 $ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
507 $ sinl )
508 m = m - 2
509 GO TO 60
510 END IF
511*
512* If working on new submatrix, choose shift direction
513* (from larger end diagonal element towards smaller)
514*
515 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
516 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
517*
518* Chase bulge from top (big end) to bottom (small end)
519*
520 idir = 1
521 ELSE
522*
523* Chase bulge from bottom (big end) to top (small end)
524*
525 idir = 2
526 END IF
527 END IF
528*
529* Apply convergence tests
530*
531 IF( idir.EQ.1 ) THEN
532*
533* Run convergence test in forward direction
534* First apply standard test to bottom of matrix
535*
536 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
537 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
538 e( m-1 ) = 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 forward
546*
547 mu = abs( d( ll ) )
548 smin = mu
549 DO 100 lll = ll, m - 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+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
555 smin = min( smin, mu )
556 100 CONTINUE
557 END IF
558*
559 ELSE
560*
561* Run convergence test in backward direction
562* First apply standard test to top of matrix
563*
564 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
565 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
566 e( ll ) = zero
567 GO TO 60
568 END IF
569*
570 IF( tol.GE.zero ) THEN
571*
572* If relative accuracy desired,
573* apply convergence criterion backward
574*
575 mu = abs( d( m ) )
576 smin = mu
577 DO 110 lll = m - 1, ll, -1
578 IF( abs( e( lll ) ).LE.tol*mu ) THEN
579 e( lll ) = zero
580 GO TO 60
581 END IF
582 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
583 smin = min( smin, mu )
584 110 CONTINUE
585 END IF
586 END IF
587 oldll = ll
588 oldm = m
589*
590* Compute shift. First, test if shifting would ruin relative
591* accuracy, and if so set the shift to zero.
592*
593 IF( tol.GE.zero .AND. real( n )*tol*( smin / smax ).LE.
594 $ max( eps, hndrth*tol ) ) THEN
595*
596* Use a zero shift to avoid loss of relative accuracy
597*
598 shift = zero
599 ELSE
600*
601* Compute the shift from 2-by-2 block at end of matrix
602*
603 IF( idir.EQ.1 ) THEN
604 sll = abs( d( ll ) )
605 CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
606 ELSE
607 sll = abs( d( m ) )
608 CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
609 END IF
610*
611* Test if shift negligible, and if so set to zero
612*
613 IF( sll.GT.zero ) THEN
614 IF( ( shift / sll )**2.LT.eps )
615 $ shift = zero
616 END IF
617 END IF
618*
619* Increment iteration count
620*
621 iter = iter + m - ll
622*
623* If SHIFT = 0, do simplified QR iteration
624*
625 IF( shift.EQ.zero ) THEN
626 IF( idir.EQ.1 ) THEN
627*
628* Chase bulge from top to bottom
629* Save cosines and sines for later singular vector updates
630*
631 cs = one
632 oldcs = one
633 DO 120 i = ll, m - 1
634 CALL slartg( d( i )*cs, e( i ), cs, sn, r )
635 IF( i.GT.ll )
636 $ e( i-1 ) = oldsn*r
637 CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn,
638 $ d( i ) )
639 work( i-ll+1 ) = cs
640 work( i-ll+1+nm1 ) = sn
641 work( i-ll+1+nm12 ) = oldcs
642 work( i-ll+1+nm13 ) = oldsn
643 120 CONTINUE
644 h = d( m )*cs
645 d( m ) = h*oldcs
646 e( m-1 ) = h*oldsn
647*
648* Update singular vectors
649*
650 IF( ncvt.GT.0 )
651 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
652 $ work( n ), vt( ll, 1 ), ldvt )
653 IF( nru.GT.0 )
654 $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1,
655 $ work( nm12+1 ),
656 $ work( nm13+1 ), u( 1, ll ), ldu )
657 IF( ncc.GT.0 )
658 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc,
659 $ work( nm12+1 ),
660 $ work( nm13+1 ), c( ll, 1 ), ldc )
661*
662* Test convergence
663*
664 IF( abs( e( m-1 ) ).LE.thresh )
665 $ e( m-1 ) = zero
666*
667 ELSE
668*
669* Chase bulge from bottom to top
670* Save cosines and sines for later singular vector updates
671*
672 cs = one
673 oldcs = one
674 DO 130 i = m, ll + 1, -1
675 CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
676 IF( i.LT.m )
677 $ e( i ) = oldsn*r
678 CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn,
679 $ d( i ) )
680 work( i-ll ) = cs
681 work( i-ll+nm1 ) = -sn
682 work( i-ll+nm12 ) = oldcs
683 work( i-ll+nm13 ) = -oldsn
684 130 CONTINUE
685 h = d( ll )*cs
686 d( ll ) = h*oldcs
687 e( ll ) = h*oldsn
688*
689* Update singular vectors
690*
691 IF( ncvt.GT.0 )
692 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt,
693 $ work( nm12+1 ),
694 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
695 IF( nru.GT.0 )
696 $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
697 $ work( n ), u( 1, ll ), ldu )
698 IF( ncc.GT.0 )
699 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
700 $ work( n ), c( ll, 1 ), ldc )
701*
702* Test convergence
703*
704 IF( abs( e( ll ) ).LE.thresh )
705 $ e( ll ) = zero
706 END IF
707 ELSE
708*
709* Use nonzero shift
710*
711 IF( idir.EQ.1 ) THEN
712*
713* Chase bulge from top to bottom
714* Save cosines and sines for later singular vector updates
715*
716 f = ( abs( d( ll ) )-shift )*
717 $ ( sign( one, d( ll ) )+shift / d( ll ) )
718 g = e( ll )
719 DO 140 i = ll, m - 1
720 CALL slartg( f, g, cosr, sinr, r )
721 IF( i.GT.ll )
722 $ e( i-1 ) = r
723 f = cosr*d( i ) + sinr*e( i )
724 e( i ) = cosr*e( i ) - sinr*d( i )
725 g = sinr*d( i+1 )
726 d( i+1 ) = cosr*d( i+1 )
727 CALL slartg( f, g, cosl, sinl, r )
728 d( i ) = r
729 f = cosl*e( i ) + sinl*d( i+1 )
730 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
731 IF( i.LT.m-1 ) THEN
732 g = sinl*e( i+1 )
733 e( i+1 ) = cosl*e( i+1 )
734 END IF
735 work( i-ll+1 ) = cosr
736 work( i-ll+1+nm1 ) = sinr
737 work( i-ll+1+nm12 ) = cosl
738 work( i-ll+1+nm13 ) = sinl
739 140 CONTINUE
740 e( m-1 ) = f
741*
742* Update singular vectors
743*
744 IF( ncvt.GT.0 )
745 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
746 $ work( n ), vt( ll, 1 ), ldvt )
747 IF( nru.GT.0 )
748 $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1,
749 $ work( nm12+1 ),
750 $ work( nm13+1 ), u( 1, ll ), ldu )
751 IF( ncc.GT.0 )
752 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc,
753 $ work( nm12+1 ),
754 $ work( nm13+1 ), c( ll, 1 ), ldc )
755*
756* Test convergence
757*
758 IF( abs( e( m-1 ) ).LE.thresh )
759 $ e( m-1 ) = zero
760*
761 ELSE
762*
763* Chase bulge from bottom to top
764* Save cosines and sines for later singular vector updates
765*
766 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
767 $ d( m ) )
768 g = e( m-1 )
769 DO 150 i = m, ll + 1, -1
770 CALL slartg( f, g, cosr, sinr, r )
771 IF( i.LT.m )
772 $ e( i ) = r
773 f = cosr*d( i ) + sinr*e( i-1 )
774 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
775 g = sinr*d( i-1 )
776 d( i-1 ) = cosr*d( i-1 )
777 CALL slartg( f, g, cosl, sinl, r )
778 d( i ) = r
779 f = cosl*e( i-1 ) + sinl*d( i-1 )
780 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
781 IF( i.GT.ll+1 ) THEN
782 g = sinl*e( i-2 )
783 e( i-2 ) = cosl*e( i-2 )
784 END IF
785 work( i-ll ) = cosr
786 work( i-ll+nm1 ) = -sinr
787 work( i-ll+nm12 ) = cosl
788 work( i-ll+nm13 ) = -sinl
789 150 CONTINUE
790 e( ll ) = f
791*
792* Test convergence
793*
794 IF( abs( e( ll ) ).LE.thresh )
795 $ e( ll ) = zero
796*
797* Update singular vectors if desired
798*
799 IF( ncvt.GT.0 )
800 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt,
801 $ work( nm12+1 ),
802 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
803 IF( nru.GT.0 )
804 $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
805 $ work( n ), u( 1, ll ), ldu )
806 IF( ncc.GT.0 )
807 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
808 $ work( n ), c( ll, 1 ), ldc )
809 END IF
810 END IF
811*
812* QR iteration finished, go back and check convergence
813*
814 GO TO 60
815*
816* All singular values converged, so make them positive
817*
818 160 CONTINUE
819 DO 170 i = 1, n
820 IF( d( i ).EQ.zero ) THEN
821*
822* Avoid -ZERO
823*
824 d( i ) = zero
825 END IF
826 IF( d( i ).LT.zero ) THEN
827
828 d( i ) = -d( i )
829*
830* Change sign of singular vectors, if desired
831*
832 IF( ncvt.GT.0 )
833 $ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
834 END IF
835 170 CONTINUE
836*
837* Sort the singular values into decreasing order (insertion sort on
838* singular values, but only one transposition per singular vector)
839*
840 DO 190 i = 1, n - 1
841*
842* Scan for smallest D(I)
843*
844 isub = 1
845 smin = d( 1 )
846 DO 180 j = 2, n + 1 - i
847 IF( d( j ).LE.smin ) THEN
848 isub = j
849 smin = d( j )
850 END IF
851 180 CONTINUE
852 IF( isub.NE.n+1-i ) THEN
853*
854* Swap singular values and vectors
855*
856 d( isub ) = d( n+1-i )
857 d( n+1-i ) = smin
858 IF( ncvt.GT.0 )
859 $ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
860 $ ldvt )
861 IF( nru.GT.0 )
862 $ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
863 IF( ncc.GT.0 )
864 $ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ),
865 $ ldc )
866 END IF
867 190 CONTINUE
868 GO TO 220
869*
870* Maximum number of iterations exceeded, failure to converge
871*
872 200 CONTINUE
873 info = 0
874 DO 210 i = 1, n - 1
875 IF( e( i ).NE.zero )
876 $ info = info + 1
877 210 CONTINUE
878 220 CONTINUE
879 RETURN
880*
881* End of SBDSQR
882*
883 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
SBDSQR
Definition sbdsqr.f:240
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:111
subroutine slas2(f, g, h, ssmin, ssmax)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition slas2.f:103
subroutine slasq1(n, d, e, work, info)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition slasq1.f:106
subroutine slasr(side, pivot, direct, m, n, c, s, a, lda)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition slasr.f:197
subroutine slasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition slasv2.f:134
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82