LAPACK 3.12.0
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*> \htmlonly
9*> Download SBDSQR + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sbdsqr.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sbdsqr.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sbdsqr.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SBDSQR( 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* REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
30* $ VT( LDVT, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SBDSQR 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 SGEBRD, 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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:
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, srot,
290 $ sscal, sswap, xerbla
291* ..
292* .. Intrinsic Functions ..
293 INTRINSIC abs, max, min, real, sign, sqrt
294* ..
295* .. Executable Statements ..
296*
297* Test the input parameters.
298*
299 info = 0
300 lower = lsame( uplo, 'L' )
301 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
302 info = -1
303 ELSE IF( n.LT.0 ) THEN
304 info = -2
305 ELSE IF( ncvt.LT.0 ) THEN
306 info = -3
307 ELSE IF( nru.LT.0 ) THEN
308 info = -4
309 ELSE IF( ncc.LT.0 ) THEN
310 info = -5
311 ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
312 $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
313 info = -9
314 ELSE IF( ldu.LT.max( 1, nru ) ) THEN
315 info = -11
316 ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
317 $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
318 info = -13
319 END IF
320 IF( info.NE.0 ) THEN
321 CALL xerbla( 'SBDSQR', -info )
322 RETURN
323 END IF
324 IF( n.EQ.0 )
325 $ RETURN
326 IF( n.EQ.1 )
327 $ GO TO 160
328*
329* ROTATE is true if any singular vectors desired, false otherwise
330*
331 rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
332*
333* If no singular vectors desired, use qd algorithm
334*
335 IF( .NOT.rotate ) THEN
336 CALL slasq1( n, d, e, work, info )
337*
338* If INFO equals 2, dqds didn't finish, try to finish
339*
340 IF( info .NE. 2 ) RETURN
341 info = 0
342 END IF
343*
344 nm1 = n - 1
345 nm12 = nm1 + nm1
346 nm13 = nm12 + nm1
347 idir = 0
348*
349* Get machine constants
350*
351 eps = slamch( 'Epsilon' )
352 unfl = slamch( 'Safe minimum' )
353*
354* If matrix lower bidiagonal, rotate to be upper bidiagonal
355* by applying Givens rotations on the left
356*
357 IF( lower ) THEN
358 DO 10 i = 1, n - 1
359 CALL slartg( d( i ), e( i ), cs, sn, r )
360 d( i ) = r
361 e( i ) = sn*d( i+1 )
362 d( i+1 ) = cs*d( i+1 )
363 work( i ) = cs
364 work( nm1+i ) = sn
365 10 CONTINUE
366*
367* Update singular vectors if desired
368*
369 IF( nru.GT.0 )
370 $ CALL slasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
371 $ ldu )
372 IF( ncc.GT.0 )
373 $ CALL slasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), c,
374 $ ldc )
375 END IF
376*
377* Compute singular values to relative accuracy TOL
378* (By setting TOL to be negative, algorithm will compute
379* singular values to absolute accuracy ABS(TOL)*norm(input matrix))
380*
381 tolmul = max( ten, min( hndrd, eps**meigth ) )
382 tol = tolmul*eps
383*
384* Compute approximate maximum, minimum singular values
385*
386 smax = zero
387 DO 20 i = 1, n
388 smax = max( smax, abs( d( i ) ) )
389 20 CONTINUE
390 DO 30 i = 1, n - 1
391 smax = max( smax, abs( e( i ) ) )
392 30 CONTINUE
393 smin = zero
394 IF( tol.GE.zero ) THEN
395*
396* Relative accuracy desired
397*
398 sminoa = abs( d( 1 ) )
399 IF( sminoa.EQ.zero )
400 $ GO TO 50
401 mu = sminoa
402 DO 40 i = 2, n
403 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
404 sminoa = min( sminoa, mu )
405 IF( sminoa.EQ.zero )
406 $ GO TO 50
407 40 CONTINUE
408 50 CONTINUE
409 sminoa = sminoa / sqrt( real( n ) )
410 thresh = max( tol*sminoa, maxitr*(n*(n*unfl)) )
411 ELSE
412*
413* Absolute accuracy desired
414*
415 thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
416 END IF
417*
418* Prepare for main iteration loop for the singular values
419* (MAXIT is the maximum number of passes through the inner
420* loop permitted before nonconvergence signalled.)
421*
422 maxitdivn = maxitr*n
423 iterdivn = 0
424 iter = -1
425 oldll = -1
426 oldm = -1
427*
428* M points to last element of unconverged part of matrix
429*
430 m = n
431*
432* Begin main iteration loop
433*
434 60 CONTINUE
435*
436* Check for convergence or exceeding iteration count
437*
438 IF( m.LE.1 )
439 $ GO TO 160
440*
441 IF( iter.GE.n ) THEN
442 iter = iter - n
443 iterdivn = iterdivn + 1
444 IF( iterdivn.GE.maxitdivn )
445 $ GO TO 200
446 END IF
447*
448* Find diagonal block of matrix to work on
449*
450 IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
451 $ d( m ) = zero
452 smax = abs( d( m ) )
453 DO 70 lll = 1, m - 1
454 ll = m - lll
455 abss = abs( d( ll ) )
456 abse = abs( e( ll ) )
457 IF( tol.LT.zero .AND. abss.LE.thresh )
458 $ d( ll ) = zero
459 IF( abse.LE.thresh )
460 $ GO TO 80
461 smax = max( smax, abss, abse )
462 70 CONTINUE
463 ll = 0
464 GO TO 90
465 80 CONTINUE
466 e( ll ) = zero
467*
468* Matrix splits since E(LL) = 0
469*
470 IF( ll.EQ.m-1 ) THEN
471*
472* Convergence of bottom singular value, return to top of loop
473*
474 m = m - 1
475 GO TO 60
476 END IF
477 90 CONTINUE
478 ll = ll + 1
479*
480* E(LL) through E(M-1) are nonzero, E(LL-1) is zero
481*
482 IF( ll.EQ.m-1 ) THEN
483*
484* 2 by 2 block, handle separately
485*
486 CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
487 $ cosr, sinl, cosl )
488 d( m-1 ) = sigmx
489 e( m-1 ) = zero
490 d( m ) = sigmn
491*
492* Compute singular vectors, if desired
493*
494 IF( ncvt.GT.0 )
495 $ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
496 $ sinr )
497 IF( nru.GT.0 )
498 $ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
499 IF( ncc.GT.0 )
500 $ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
501 $ sinl )
502 m = m - 2
503 GO TO 60
504 END IF
505*
506* If working on new submatrix, choose shift direction
507* (from larger end diagonal element towards smaller)
508*
509 IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
510 IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
511*
512* Chase bulge from top (big end) to bottom (small end)
513*
514 idir = 1
515 ELSE
516*
517* Chase bulge from bottom (big end) to top (small end)
518*
519 idir = 2
520 END IF
521 END IF
522*
523* Apply convergence tests
524*
525 IF( idir.EQ.1 ) THEN
526*
527* Run convergence test in forward direction
528* First apply standard test to bottom of matrix
529*
530 IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
531 $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
532 e( m-1 ) = zero
533 GO TO 60
534 END IF
535*
536 IF( tol.GE.zero ) THEN
537*
538* If relative accuracy desired,
539* apply convergence criterion forward
540*
541 mu = abs( d( ll ) )
542 smin = mu
543 DO 100 lll = ll, m - 1
544 IF( abs( e( lll ) ).LE.tol*mu ) THEN
545 e( lll ) = zero
546 GO TO 60
547 END IF
548 mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
549 smin = min( smin, mu )
550 100 CONTINUE
551 END IF
552*
553 ELSE
554*
555* Run convergence test in backward direction
556* First apply standard test to top of matrix
557*
558 IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
559 $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
560 e( ll ) = zero
561 GO TO 60
562 END IF
563*
564 IF( tol.GE.zero ) THEN
565*
566* If relative accuracy desired,
567* apply convergence criterion backward
568*
569 mu = abs( d( m ) )
570 smin = mu
571 DO 110 lll = m - 1, ll, -1
572 IF( abs( e( lll ) ).LE.tol*mu ) THEN
573 e( lll ) = zero
574 GO TO 60
575 END IF
576 mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
577 smin = min( smin, mu )
578 110 CONTINUE
579 END IF
580 END IF
581 oldll = ll
582 oldm = m
583*
584* Compute shift. First, test if shifting would ruin relative
585* accuracy, and if so set the shift to zero.
586*
587 IF( tol.GE.zero .AND. n*tol*( smin / smax ).LE.
588 $ max( eps, hndrth*tol ) ) THEN
589*
590* Use a zero shift to avoid loss of relative accuracy
591*
592 shift = zero
593 ELSE
594*
595* Compute the shift from 2-by-2 block at end of matrix
596*
597 IF( idir.EQ.1 ) THEN
598 sll = abs( d( ll ) )
599 CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
600 ELSE
601 sll = abs( d( m ) )
602 CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
603 END IF
604*
605* Test if shift negligible, and if so set to zero
606*
607 IF( sll.GT.zero ) THEN
608 IF( ( shift / sll )**2.LT.eps )
609 $ shift = zero
610 END IF
611 END IF
612*
613* Increment iteration count
614*
615 iter = iter + m - ll
616*
617* If SHIFT = 0, do simplified QR iteration
618*
619 IF( shift.EQ.zero ) THEN
620 IF( idir.EQ.1 ) THEN
621*
622* Chase bulge from top to bottom
623* Save cosines and sines for later singular vector updates
624*
625 cs = one
626 oldcs = one
627 DO 120 i = ll, m - 1
628 CALL slartg( d( i )*cs, e( i ), cs, sn, r )
629 IF( i.GT.ll )
630 $ e( i-1 ) = oldsn*r
631 CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
632 work( i-ll+1 ) = cs
633 work( i-ll+1+nm1 ) = sn
634 work( i-ll+1+nm12 ) = oldcs
635 work( i-ll+1+nm13 ) = oldsn
636 120 CONTINUE
637 h = d( m )*cs
638 d( m ) = h*oldcs
639 e( m-1 ) = h*oldsn
640*
641* Update singular vectors
642*
643 IF( ncvt.GT.0 )
644 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
645 $ work( n ), vt( ll, 1 ), ldvt )
646 IF( nru.GT.0 )
647 $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
648 $ work( nm13+1 ), u( 1, ll ), ldu )
649 IF( ncc.GT.0 )
650 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
651 $ work( nm13+1 ), c( ll, 1 ), ldc )
652*
653* Test convergence
654*
655 IF( abs( e( m-1 ) ).LE.thresh )
656 $ e( m-1 ) = zero
657*
658 ELSE
659*
660* Chase bulge from bottom to top
661* Save cosines and sines for later singular vector updates
662*
663 cs = one
664 oldcs = one
665 DO 130 i = m, ll + 1, -1
666 CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
667 IF( i.LT.m )
668 $ e( i ) = oldsn*r
669 CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
670 work( i-ll ) = cs
671 work( i-ll+nm1 ) = -sn
672 work( i-ll+nm12 ) = oldcs
673 work( 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 slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
683 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
684 IF( nru.GT.0 )
685 $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
686 $ work( n ), u( 1, ll ), ldu )
687 IF( ncc.GT.0 )
688 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
689 $ work( n ), c( ll, 1 ), ldc )
690*
691* Test convergence
692*
693 IF( abs( e( ll ) ).LE.thresh )
694 $ e( ll ) = zero
695 END IF
696 ELSE
697*
698* Use nonzero shift
699*
700 IF( idir.EQ.1 ) THEN
701*
702* Chase bulge from top to bottom
703* Save cosines and sines for later singular vector updates
704*
705 f = ( abs( d( ll ) )-shift )*
706 $ ( sign( one, d( ll ) )+shift / d( ll ) )
707 g = e( ll )
708 DO 140 i = ll, m - 1
709 CALL slartg( f, g, cosr, sinr, r )
710 IF( i.GT.ll )
711 $ e( i-1 ) = r
712 f = cosr*d( i ) + sinr*e( i )
713 e( i ) = cosr*e( i ) - sinr*d( i )
714 g = sinr*d( i+1 )
715 d( i+1 ) = cosr*d( i+1 )
716 CALL slartg( f, g, cosl, sinl, r )
717 d( i ) = r
718 f = cosl*e( i ) + sinl*d( i+1 )
719 d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
720 IF( i.LT.m-1 ) THEN
721 g = sinl*e( i+1 )
722 e( i+1 ) = cosl*e( i+1 )
723 END IF
724 work( i-ll+1 ) = cosr
725 work( i-ll+1+nm1 ) = sinr
726 work( i-ll+1+nm12 ) = cosl
727 work( i-ll+1+nm13 ) = sinl
728 140 CONTINUE
729 e( m-1 ) = f
730*
731* Update singular vectors
732*
733 IF( ncvt.GT.0 )
734 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
735 $ work( n ), vt( ll, 1 ), ldvt )
736 IF( nru.GT.0 )
737 $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
738 $ work( nm13+1 ), u( 1, ll ), ldu )
739 IF( ncc.GT.0 )
740 $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
741 $ work( nm13+1 ), c( ll, 1 ), ldc )
742*
743* Test convergence
744*
745 IF( abs( e( m-1 ) ).LE.thresh )
746 $ e( m-1 ) = zero
747*
748 ELSE
749*
750* Chase bulge from bottom to top
751* Save cosines and sines for later singular vector updates
752*
753 f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
754 $ d( m ) )
755 g = e( m-1 )
756 DO 150 i = m, ll + 1, -1
757 CALL slartg( f, g, cosr, sinr, r )
758 IF( i.LT.m )
759 $ e( i ) = r
760 f = cosr*d( i ) + sinr*e( i-1 )
761 e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
762 g = sinr*d( i-1 )
763 d( i-1 ) = cosr*d( i-1 )
764 CALL slartg( f, g, cosl, sinl, r )
765 d( i ) = r
766 f = cosl*e( i-1 ) + sinl*d( i-1 )
767 d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
768 IF( i.GT.ll+1 ) THEN
769 g = sinl*e( i-2 )
770 e( i-2 ) = cosl*e( i-2 )
771 END IF
772 work( i-ll ) = cosr
773 work( i-ll+nm1 ) = -sinr
774 work( i-ll+nm12 ) = cosl
775 work( i-ll+nm13 ) = -sinl
776 150 CONTINUE
777 e( ll ) = f
778*
779* Test convergence
780*
781 IF( abs( e( ll ) ).LE.thresh )
782 $ e( ll ) = zero
783*
784* Update singular vectors if desired
785*
786 IF( ncvt.GT.0 )
787 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
788 $ work( nm13+1 ), vt( ll, 1 ), ldvt )
789 IF( nru.GT.0 )
790 $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
791 $ work( n ), u( 1, ll ), ldu )
792 IF( ncc.GT.0 )
793 $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
794 $ work( n ), c( ll, 1 ), ldc )
795 END IF
796 END IF
797*
798* QR iteration finished, go back and check convergence
799*
800 GO TO 60
801*
802* All singular values converged, so make them positive
803*
804 160 CONTINUE
805 DO 170 i = 1, n
806 IF( d( i ).LT.zero ) THEN
807 d( i ) = -d( i )
808*
809* Change sign of singular vectors, if desired
810*
811 IF( ncvt.GT.0 )
812 $ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
813 END IF
814 170 CONTINUE
815*
816* Sort the singular values into decreasing order (insertion sort on
817* singular values, but only one transposition per singular vector)
818*
819 DO 190 i = 1, n - 1
820*
821* Scan for smallest D(I)
822*
823 isub = 1
824 smin = d( 1 )
825 DO 180 j = 2, n + 1 - i
826 IF( d( j ).LE.smin ) THEN
827 isub = j
828 smin = d( j )
829 END IF
830 180 CONTINUE
831 IF( isub.NE.n+1-i ) THEN
832*
833* Swap singular values and vectors
834*
835 d( isub ) = d( n+1-i )
836 d( n+1-i ) = smin
837 IF( ncvt.GT.0 )
838 $ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
839 $ ldvt )
840 IF( nru.GT.0 )
841 $ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
842 IF( ncc.GT.0 )
843 $ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
844 END IF
845 190 CONTINUE
846 GO TO 220
847*
848* Maximum number of iterations exceeded, failure to converge
849*
850 200 CONTINUE
851 info = 0
852 DO 210 i = 1, n - 1
853 IF( e( i ).NE.zero )
854 $ info = info + 1
855 210 CONTINUE
856 220 CONTINUE
857 RETURN
858*
859* End of SBDSQR
860*
861 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:105
subroutine slasq1(n, d, e, work, info)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition slasq1.f:108
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:199
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:136
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