LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
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 *> \date June 2017
236 *
237 *> \ingroup auxOTHERcomputational
238 *
239 * =====================================================================
240  SUBROUTINE sbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
241  $ LDU, C, LDC, WORK, INFO )
242 *
243 * -- LAPACK computational routine (version 3.7.1) --
244 * -- LAPACK is a software package provided by Univ. of Tennessee, --
245 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246 * June 2017
247 *
248 * .. Scalar Arguments ..
249  CHARACTER UPLO
250  INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
251 * ..
252 * .. Array Arguments ..
253  REAL C( LDC, * ), D( * ), E( * ), U( LDU, * ),
254  $ vt( ldvt, * ), work( * )
255 * ..
256 *
257 * =====================================================================
258 *
259 * .. Parameters ..
260  REAL ZERO
261  parameter( zero = 0.0e0 )
262  REAL ONE
263  parameter( one = 1.0e0 )
264  REAL NEGONE
265  parameter( negone = -1.0e0 )
266  REAL HNDRTH
267  parameter( hndrth = 0.01e0 )
268  REAL TEN
269  parameter( ten = 10.0e0 )
270  REAL HNDRD
271  parameter( hndrd = 100.0e0 )
272  REAL MEIGTH
273  parameter( meigth = -0.125e0 )
274  INTEGER MAXITR
275  parameter( maxitr = 6 )
276 * ..
277 * .. Local Scalars ..
278  LOGICAL LOWER, ROTATE
279  INTEGER I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
280  $ maxitdivn, nm1, nm12, nm13, oldll, oldm
281  REAL ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
282  $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
283  $ sinr, sll, smax, smin, sminl, sminoa,
284  $ sn, thresh, tol, tolmul, unfl
285 * ..
286 * .. External Functions ..
287  LOGICAL LSAME
288  REAL SLAMCH
289  EXTERNAL lsame, slamch
290 * ..
291 * .. External Subroutines ..
292  EXTERNAL slartg, slas2, slasq1, slasr, slasv2, srot,
293  $ sscal, sswap, xerbla
294 * ..
295 * .. Intrinsic Functions ..
296  INTRINSIC abs, max, min, real, sign, sqrt
297 * ..
298 * .. Executable Statements ..
299 *
300 * Test the input parameters.
301 *
302  info = 0
303  lower = lsame( uplo, 'L' )
304  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
305  info = -1
306  ELSE IF( n.LT.0 ) THEN
307  info = -2
308  ELSE IF( ncvt.LT.0 ) THEN
309  info = -3
310  ELSE IF( nru.LT.0 ) THEN
311  info = -4
312  ELSE IF( ncc.LT.0 ) THEN
313  info = -5
314  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
315  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
316  info = -9
317  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
318  info = -11
319  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
320  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
321  info = -13
322  END IF
323  IF( info.NE.0 ) THEN
324  CALL xerbla( 'SBDSQR', -info )
325  RETURN
326  END IF
327  IF( n.EQ.0 )
328  $ RETURN
329  IF( n.EQ.1 )
330  $ GO TO 160
331 *
332 * ROTATE is true if any singular vectors desired, false otherwise
333 *
334  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
335 *
336 * If no singular vectors desired, use qd algorithm
337 *
338  IF( .NOT.rotate ) THEN
339  CALL slasq1( n, d, e, work, info )
340 *
341 * If INFO equals 2, dqds didn't finish, try to finish
342 *
343  IF( info .NE. 2 ) RETURN
344  info = 0
345  END IF
346 *
347  nm1 = n - 1
348  nm12 = nm1 + nm1
349  nm13 = nm12 + nm1
350  idir = 0
351 *
352 * Get machine constants
353 *
354  eps = slamch( 'Epsilon' )
355  unfl = slamch( 'Safe minimum' )
356 *
357 * If matrix lower bidiagonal, rotate to be upper bidiagonal
358 * by applying Givens rotations on the left
359 *
360  IF( lower ) THEN
361  DO 10 i = 1, n - 1
362  CALL slartg( d( i ), e( i ), cs, sn, r )
363  d( i ) = r
364  e( i ) = sn*d( i+1 )
365  d( i+1 ) = cs*d( i+1 )
366  work( i ) = cs
367  work( nm1+i ) = sn
368  10 CONTINUE
369 *
370 * Update singular vectors if desired
371 *
372  IF( nru.GT.0 )
373  $ CALL slasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
374  $ ldu )
375  IF( ncc.GT.0 )
376  $ CALL slasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), 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  sminl = 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*(n*(n*unfl)) )
414  ELSE
415 *
416 * Absolute accuracy desired
417 *
418  thresh = max( abs( tol )*smax, maxitr*(n*(n*unfl)) )
419  END IF
420 *
421 * Prepare for main iteration loop for the singular values
422 * (MAXIT is the maximum number of passes through the inner
423 * loop permitted before nonconvergence signalled.)
424 *
425  maxitdivn = maxitr*n
426  iterdivn = 0
427  iter = -1
428  oldll = -1
429  oldm = -1
430 *
431 * M points to last element of unconverged part of matrix
432 *
433  m = n
434 *
435 * Begin main iteration loop
436 *
437  60 CONTINUE
438 *
439 * Check for convergence or exceeding iteration count
440 *
441  IF( m.LE.1 )
442  $ GO TO 160
443 *
444  IF( iter.GE.n ) THEN
445  iter = iter - n
446  iterdivn = iterdivn + 1
447  IF( iterdivn.GE.maxitdivn )
448  $ GO TO 200
449  END IF
450 *
451 * Find diagonal block of matrix to work on
452 *
453  IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
454  $ d( m ) = zero
455  smax = abs( d( m ) )
456  smin = smax
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  smin = min( smin, abss )
466  smax = max( smax, abss, abse )
467  70 CONTINUE
468  ll = 0
469  GO TO 90
470  80 CONTINUE
471  e( ll ) = zero
472 *
473 * Matrix splits since E(LL) = 0
474 *
475  IF( ll.EQ.m-1 ) THEN
476 *
477 * Convergence of bottom singular value, return to top of loop
478 *
479  m = m - 1
480  GO TO 60
481  END IF
482  90 CONTINUE
483  ll = ll + 1
484 *
485 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
486 *
487  IF( ll.EQ.m-1 ) THEN
488 *
489 * 2 by 2 block, handle separately
490 *
491  CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
492  $ cosr, sinl, cosl )
493  d( m-1 ) = sigmx
494  e( m-1 ) = zero
495  d( m ) = sigmn
496 *
497 * Compute singular vectors, if desired
498 *
499  IF( ncvt.GT.0 )
500  $ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
501  $ sinr )
502  IF( nru.GT.0 )
503  $ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
504  IF( ncc.GT.0 )
505  $ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
506  $ sinl )
507  m = m - 2
508  GO TO 60
509  END IF
510 *
511 * If working on new submatrix, choose shift direction
512 * (from larger end diagonal element towards smaller)
513 *
514  IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
515  IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
516 *
517 * Chase bulge from top (big end) to bottom (small end)
518 *
519  idir = 1
520  ELSE
521 *
522 * Chase bulge from bottom (big end) to top (small end)
523 *
524  idir = 2
525  END IF
526  END IF
527 *
528 * Apply convergence tests
529 *
530  IF( idir.EQ.1 ) THEN
531 *
532 * Run convergence test in forward direction
533 * First apply standard test to bottom of matrix
534 *
535  IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
536  $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
537  e( m-1 ) = zero
538  GO TO 60
539  END IF
540 *
541  IF( tol.GE.zero ) THEN
542 *
543 * If relative accuracy desired,
544 * apply convergence criterion forward
545 *
546  mu = abs( d( ll ) )
547  sminl = mu
548  DO 100 lll = ll, m - 1
549  IF( abs( e( lll ) ).LE.tol*mu ) THEN
550  e( lll ) = zero
551  GO TO 60
552  END IF
553  mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
554  sminl = min( sminl, mu )
555  100 CONTINUE
556  END IF
557 *
558  ELSE
559 *
560 * Run convergence test in backward direction
561 * First apply standard test to top of matrix
562 *
563  IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
564  $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
565  e( ll ) = zero
566  GO TO 60
567  END IF
568 *
569  IF( tol.GE.zero ) THEN
570 *
571 * If relative accuracy desired,
572 * apply convergence criterion backward
573 *
574  mu = abs( d( m ) )
575  sminl = mu
576  DO 110 lll = m - 1, ll, -1
577  IF( abs( e( lll ) ).LE.tol*mu ) THEN
578  e( lll ) = zero
579  GO TO 60
580  END IF
581  mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
582  sminl = min( sminl, mu )
583  110 CONTINUE
584  END IF
585  END IF
586  oldll = ll
587  oldm = m
588 *
589 * Compute shift. First, test if shifting would ruin relative
590 * accuracy, and if so set the shift to zero.
591 *
592  IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
593  $ max( eps, hndrth*tol ) ) THEN
594 *
595 * Use a zero shift to avoid loss of relative accuracy
596 *
597  shift = zero
598  ELSE
599 *
600 * Compute the shift from 2-by-2 block at end of matrix
601 *
602  IF( idir.EQ.1 ) THEN
603  sll = abs( d( ll ) )
604  CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
605  ELSE
606  sll = abs( d( m ) )
607  CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
608  END IF
609 *
610 * Test if shift negligible, and if so set to zero
611 *
612  IF( sll.GT.zero ) THEN
613  IF( ( shift / sll )**2.LT.eps )
614  $ shift = zero
615  END IF
616  END IF
617 *
618 * Increment iteration count
619 *
620  iter = iter + m - ll
621 *
622 * If SHIFT = 0, do simplified QR iteration
623 *
624  IF( shift.EQ.zero ) THEN
625  IF( idir.EQ.1 ) THEN
626 *
627 * Chase bulge from top to bottom
628 * Save cosines and sines for later singular vector updates
629 *
630  cs = one
631  oldcs = one
632  DO 120 i = ll, m - 1
633  CALL slartg( d( i )*cs, e( i ), cs, sn, r )
634  IF( i.GT.ll )
635  $ e( i-1 ) = oldsn*r
636  CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
637  work( i-ll+1 ) = cs
638  work( i-ll+1+nm1 ) = sn
639  work( i-ll+1+nm12 ) = oldcs
640  work( i-ll+1+nm13 ) = oldsn
641  120 CONTINUE
642  h = d( m )*cs
643  d( m ) = h*oldcs
644  e( m-1 ) = h*oldsn
645 *
646 * Update singular vectors
647 *
648  IF( ncvt.GT.0 )
649  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
650  $ work( n ), vt( ll, 1 ), ldvt )
651  IF( nru.GT.0 )
652  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
653  $ work( nm13+1 ), u( 1, ll ), ldu )
654  IF( ncc.GT.0 )
655  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
656  $ work( nm13+1 ), c( ll, 1 ), ldc )
657 *
658 * Test convergence
659 *
660  IF( abs( e( m-1 ) ).LE.thresh )
661  $ e( m-1 ) = zero
662 *
663  ELSE
664 *
665 * Chase bulge from bottom to top
666 * Save cosines and sines for later singular vector updates
667 *
668  cs = one
669  oldcs = one
670  DO 130 i = m, ll + 1, -1
671  CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
672  IF( i.LT.m )
673  $ e( i ) = oldsn*r
674  CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
675  work( i-ll ) = cs
676  work( i-ll+nm1 ) = -sn
677  work( i-ll+nm12 ) = oldcs
678  work( i-ll+nm13 ) = -oldsn
679  130 CONTINUE
680  h = d( ll )*cs
681  d( ll ) = h*oldcs
682  e( ll ) = h*oldsn
683 *
684 * Update singular vectors
685 *
686  IF( ncvt.GT.0 )
687  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
688  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
689  IF( nru.GT.0 )
690  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
691  $ work( n ), u( 1, ll ), ldu )
692  IF( ncc.GT.0 )
693  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
694  $ work( n ), c( ll, 1 ), ldc )
695 *
696 * Test convergence
697 *
698  IF( abs( e( ll ) ).LE.thresh )
699  $ e( ll ) = zero
700  END IF
701  ELSE
702 *
703 * Use nonzero shift
704 *
705  IF( idir.EQ.1 ) THEN
706 *
707 * Chase bulge from top to bottom
708 * Save cosines and sines for later singular vector updates
709 *
710  f = ( abs( d( ll ) )-shift )*
711  $ ( sign( one, d( ll ) )+shift / d( ll ) )
712  g = e( ll )
713  DO 140 i = ll, m - 1
714  CALL slartg( f, g, cosr, sinr, r )
715  IF( i.GT.ll )
716  $ e( i-1 ) = r
717  f = cosr*d( i ) + sinr*e( i )
718  e( i ) = cosr*e( i ) - sinr*d( i )
719  g = sinr*d( i+1 )
720  d( i+1 ) = cosr*d( i+1 )
721  CALL slartg( f, g, cosl, sinl, r )
722  d( i ) = r
723  f = cosl*e( i ) + sinl*d( i+1 )
724  d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
725  IF( i.LT.m-1 ) THEN
726  g = sinl*e( i+1 )
727  e( i+1 ) = cosl*e( i+1 )
728  END IF
729  work( i-ll+1 ) = cosr
730  work( i-ll+1+nm1 ) = sinr
731  work( i-ll+1+nm12 ) = cosl
732  work( i-ll+1+nm13 ) = sinl
733  140 CONTINUE
734  e( m-1 ) = f
735 *
736 * Update singular vectors
737 *
738  IF( ncvt.GT.0 )
739  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
740  $ work( n ), vt( ll, 1 ), ldvt )
741  IF( nru.GT.0 )
742  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
743  $ work( nm13+1 ), u( 1, ll ), ldu )
744  IF( ncc.GT.0 )
745  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
746  $ work( nm13+1 ), c( ll, 1 ), ldc )
747 *
748 * Test convergence
749 *
750  IF( abs( e( m-1 ) ).LE.thresh )
751  $ e( m-1 ) = zero
752 *
753  ELSE
754 *
755 * Chase bulge from bottom to top
756 * Save cosines and sines for later singular vector updates
757 *
758  f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
759  $ d( m ) )
760  g = e( m-1 )
761  DO 150 i = m, ll + 1, -1
762  CALL slartg( f, g, cosr, sinr, r )
763  IF( i.LT.m )
764  $ e( i ) = r
765  f = cosr*d( i ) + sinr*e( i-1 )
766  e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
767  g = sinr*d( i-1 )
768  d( i-1 ) = cosr*d( i-1 )
769  CALL slartg( f, g, cosl, sinl, r )
770  d( i ) = r
771  f = cosl*e( i-1 ) + sinl*d( i-1 )
772  d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
773  IF( i.GT.ll+1 ) THEN
774  g = sinl*e( i-2 )
775  e( i-2 ) = cosl*e( i-2 )
776  END IF
777  work( i-ll ) = cosr
778  work( i-ll+nm1 ) = -sinr
779  work( i-ll+nm12 ) = cosl
780  work( i-ll+nm13 ) = -sinl
781  150 CONTINUE
782  e( ll ) = f
783 *
784 * Test convergence
785 *
786  IF( abs( e( ll ) ).LE.thresh )
787  $ e( ll ) = zero
788 *
789 * Update singular vectors if desired
790 *
791  IF( ncvt.GT.0 )
792  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
793  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
794  IF( nru.GT.0 )
795  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
796  $ work( n ), u( 1, ll ), ldu )
797  IF( ncc.GT.0 )
798  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
799  $ work( n ), c( ll, 1 ), ldc )
800  END IF
801  END IF
802 *
803 * QR iteration finished, go back and check convergence
804 *
805  GO TO 60
806 *
807 * All singular values converged, so make them positive
808 *
809  160 CONTINUE
810  DO 170 i = 1, n
811  IF( d( i ).LT.zero ) THEN
812  d( i ) = -d( i )
813 *
814 * Change sign of singular vectors, if desired
815 *
816  IF( ncvt.GT.0 )
817  $ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
818  END IF
819  170 CONTINUE
820 *
821 * Sort the singular values into decreasing order (insertion sort on
822 * singular values, but only one transposition per singular vector)
823 *
824  DO 190 i = 1, n - 1
825 *
826 * Scan for smallest D(I)
827 *
828  isub = 1
829  smin = d( 1 )
830  DO 180 j = 2, n + 1 - i
831  IF( d( j ).LE.smin ) THEN
832  isub = j
833  smin = d( j )
834  END IF
835  180 CONTINUE
836  IF( isub.NE.n+1-i ) THEN
837 *
838 * Swap singular values and vectors
839 *
840  d( isub ) = d( n+1-i )
841  d( n+1-i ) = smin
842  IF( ncvt.GT.0 )
843  $ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
844  $ ldvt )
845  IF( nru.GT.0 )
846  $ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
847  IF( ncc.GT.0 )
848  $ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
849  END IF
850  190 CONTINUE
851  GO TO 220
852 *
853 * Maximum number of iterations exceeded, failure to converge
854 *
855  200 CONTINUE
856  info = 0
857  DO 210 i = 1, n - 1
858  IF( e( i ).NE.zero )
859  $ info = info + 1
860  210 CONTINUE
861  220 CONTINUE
862  RETURN
863 *
864 * End of SBDSQR
865 *
866  END
slas2
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:109
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
sbdsqr
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:242
slasr
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:201
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:94
slasq1
subroutine slasq1(N, D, E, WORK, INFO)
SLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
Definition: slasq1.f:110
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slasv2
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:140
slartg
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f:99