LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
sbdsqr.f
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, sminl, 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  sminl = 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  smin = smax
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  smin = min( smin, abss )
463  smax = max( smax, abss, abse )
464  70 CONTINUE
465  ll = 0
466  GO TO 90
467  80 CONTINUE
468  e( ll ) = zero
469 *
470 * Matrix splits since E(LL) = 0
471 *
472  IF( ll.EQ.m-1 ) THEN
473 *
474 * Convergence of bottom singular value, return to top of loop
475 *
476  m = m - 1
477  GO TO 60
478  END IF
479  90 CONTINUE
480  ll = ll + 1
481 *
482 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
483 *
484  IF( ll.EQ.m-1 ) THEN
485 *
486 * 2 by 2 block, handle separately
487 *
488  CALL slasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
489  $ cosr, sinl, cosl )
490  d( m-1 ) = sigmx
491  e( m-1 ) = zero
492  d( m ) = sigmn
493 *
494 * Compute singular vectors, if desired
495 *
496  IF( ncvt.GT.0 )
497  $ CALL srot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
498  $ sinr )
499  IF( nru.GT.0 )
500  $ CALL srot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
501  IF( ncc.GT.0 )
502  $ CALL srot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
503  $ sinl )
504  m = m - 2
505  GO TO 60
506  END IF
507 *
508 * If working on new submatrix, choose shift direction
509 * (from larger end diagonal element towards smaller)
510 *
511  IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
512  IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
513 *
514 * Chase bulge from top (big end) to bottom (small end)
515 *
516  idir = 1
517  ELSE
518 *
519 * Chase bulge from bottom (big end) to top (small end)
520 *
521  idir = 2
522  END IF
523  END IF
524 *
525 * Apply convergence tests
526 *
527  IF( idir.EQ.1 ) THEN
528 *
529 * Run convergence test in forward direction
530 * First apply standard test to bottom of matrix
531 *
532  IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
533  $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
534  e( m-1 ) = zero
535  GO TO 60
536  END IF
537 *
538  IF( tol.GE.zero ) THEN
539 *
540 * If relative accuracy desired,
541 * apply convergence criterion forward
542 *
543  mu = abs( d( ll ) )
544  sminl = mu
545  DO 100 lll = ll, m - 1
546  IF( abs( e( lll ) ).LE.tol*mu ) THEN
547  e( lll ) = zero
548  GO TO 60
549  END IF
550  mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
551  sminl = min( sminl, mu )
552  100 CONTINUE
553  END IF
554 *
555  ELSE
556 *
557 * Run convergence test in backward direction
558 * First apply standard test to top of matrix
559 *
560  IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
561  $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
562  e( ll ) = zero
563  GO TO 60
564  END IF
565 *
566  IF( tol.GE.zero ) THEN
567 *
568 * If relative accuracy desired,
569 * apply convergence criterion backward
570 *
571  mu = abs( d( m ) )
572  sminl = mu
573  DO 110 lll = m - 1, ll, -1
574  IF( abs( e( lll ) ).LE.tol*mu ) THEN
575  e( lll ) = zero
576  GO TO 60
577  END IF
578  mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
579  sminl = min( sminl, mu )
580  110 CONTINUE
581  END IF
582  END IF
583  oldll = ll
584  oldm = m
585 *
586 * Compute shift. First, test if shifting would ruin relative
587 * accuracy, and if so set the shift to zero.
588 *
589  IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
590  $ max( eps, hndrth*tol ) ) THEN
591 *
592 * Use a zero shift to avoid loss of relative accuracy
593 *
594  shift = zero
595  ELSE
596 *
597 * Compute the shift from 2-by-2 block at end of matrix
598 *
599  IF( idir.EQ.1 ) THEN
600  sll = abs( d( ll ) )
601  CALL slas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
602  ELSE
603  sll = abs( d( m ) )
604  CALL slas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
605  END IF
606 *
607 * Test if shift negligible, and if so set to zero
608 *
609  IF( sll.GT.zero ) THEN
610  IF( ( shift / sll )**2.LT.eps )
611  $ shift = zero
612  END IF
613  END IF
614 *
615 * Increment iteration count
616 *
617  iter = iter + m - ll
618 *
619 * If SHIFT = 0, do simplified QR iteration
620 *
621  IF( shift.EQ.zero ) THEN
622  IF( idir.EQ.1 ) THEN
623 *
624 * Chase bulge from top to bottom
625 * Save cosines and sines for later singular vector updates
626 *
627  cs = one
628  oldcs = one
629  DO 120 i = ll, m - 1
630  CALL slartg( d( i )*cs, e( i ), cs, sn, r )
631  IF( i.GT.ll )
632  $ e( i-1 ) = oldsn*r
633  CALL slartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
634  work( i-ll+1 ) = cs
635  work( i-ll+1+nm1 ) = sn
636  work( i-ll+1+nm12 ) = oldcs
637  work( i-ll+1+nm13 ) = oldsn
638  120 CONTINUE
639  h = d( m )*cs
640  d( m ) = h*oldcs
641  e( m-1 ) = h*oldsn
642 *
643 * Update singular vectors
644 *
645  IF( ncvt.GT.0 )
646  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
647  $ work( n ), vt( ll, 1 ), ldvt )
648  IF( nru.GT.0 )
649  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
650  $ work( nm13+1 ), u( 1, ll ), ldu )
651  IF( ncc.GT.0 )
652  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
653  $ work( nm13+1 ), c( ll, 1 ), ldc )
654 *
655 * Test convergence
656 *
657  IF( abs( e( m-1 ) ).LE.thresh )
658  $ e( m-1 ) = zero
659 *
660  ELSE
661 *
662 * Chase bulge from bottom to top
663 * Save cosines and sines for later singular vector updates
664 *
665  cs = one
666  oldcs = one
667  DO 130 i = m, ll + 1, -1
668  CALL slartg( d( i )*cs, e( i-1 ), cs, sn, r )
669  IF( i.LT.m )
670  $ e( i ) = oldsn*r
671  CALL slartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
672  work( i-ll ) = cs
673  work( i-ll+nm1 ) = -sn
674  work( i-ll+nm12 ) = oldcs
675  work( i-ll+nm13 ) = -oldsn
676  130 CONTINUE
677  h = d( ll )*cs
678  d( ll ) = h*oldcs
679  e( ll ) = h*oldsn
680 *
681 * Update singular vectors
682 *
683  IF( ncvt.GT.0 )
684  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
685  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
686  IF( nru.GT.0 )
687  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
688  $ work( n ), u( 1, ll ), ldu )
689  IF( ncc.GT.0 )
690  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
691  $ work( n ), c( ll, 1 ), ldc )
692 *
693 * Test convergence
694 *
695  IF( abs( e( ll ) ).LE.thresh )
696  $ e( ll ) = zero
697  END IF
698  ELSE
699 *
700 * Use nonzero shift
701 *
702  IF( idir.EQ.1 ) THEN
703 *
704 * Chase bulge from top to bottom
705 * Save cosines and sines for later singular vector updates
706 *
707  f = ( abs( d( ll ) )-shift )*
708  $ ( sign( one, d( ll ) )+shift / d( ll ) )
709  g = e( ll )
710  DO 140 i = ll, m - 1
711  CALL slartg( f, g, cosr, sinr, r )
712  IF( i.GT.ll )
713  $ e( i-1 ) = r
714  f = cosr*d( i ) + sinr*e( i )
715  e( i ) = cosr*e( i ) - sinr*d( i )
716  g = sinr*d( i+1 )
717  d( i+1 ) = cosr*d( i+1 )
718  CALL slartg( f, g, cosl, sinl, r )
719  d( i ) = r
720  f = cosl*e( i ) + sinl*d( i+1 )
721  d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
722  IF( i.LT.m-1 ) THEN
723  g = sinl*e( i+1 )
724  e( i+1 ) = cosl*e( i+1 )
725  END IF
726  work( i-ll+1 ) = cosr
727  work( i-ll+1+nm1 ) = sinr
728  work( i-ll+1+nm12 ) = cosl
729  work( i-ll+1+nm13 ) = sinl
730  140 CONTINUE
731  e( m-1 ) = f
732 *
733 * Update singular vectors
734 *
735  IF( ncvt.GT.0 )
736  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
737  $ work( n ), vt( ll, 1 ), ldvt )
738  IF( nru.GT.0 )
739  $ CALL slasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
740  $ work( nm13+1 ), u( 1, ll ), ldu )
741  IF( ncc.GT.0 )
742  $ CALL slasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
743  $ work( nm13+1 ), c( ll, 1 ), ldc )
744 *
745 * Test convergence
746 *
747  IF( abs( e( m-1 ) ).LE.thresh )
748  $ e( m-1 ) = zero
749 *
750  ELSE
751 *
752 * Chase bulge from bottom to top
753 * Save cosines and sines for later singular vector updates
754 *
755  f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
756  $ d( m ) )
757  g = e( m-1 )
758  DO 150 i = m, ll + 1, -1
759  CALL slartg( f, g, cosr, sinr, r )
760  IF( i.LT.m )
761  $ e( i ) = r
762  f = cosr*d( i ) + sinr*e( i-1 )
763  e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
764  g = sinr*d( i-1 )
765  d( i-1 ) = cosr*d( i-1 )
766  CALL slartg( f, g, cosl, sinl, r )
767  d( i ) = r
768  f = cosl*e( i-1 ) + sinl*d( i-1 )
769  d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
770  IF( i.GT.ll+1 ) THEN
771  g = sinl*e( i-2 )
772  e( i-2 ) = cosl*e( i-2 )
773  END IF
774  work( i-ll ) = cosr
775  work( i-ll+nm1 ) = -sinr
776  work( i-ll+nm12 ) = cosl
777  work( i-ll+nm13 ) = -sinl
778  150 CONTINUE
779  e( ll ) = f
780 *
781 * Test convergence
782 *
783  IF( abs( e( ll ) ).LE.thresh )
784  $ e( ll ) = zero
785 *
786 * Update singular vectors if desired
787 *
788  IF( ncvt.GT.0 )
789  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
790  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
791  IF( nru.GT.0 )
792  $ CALL slasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
793  $ work( n ), u( 1, ll ), ldu )
794  IF( ncc.GT.0 )
795  $ CALL slasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
796  $ work( n ), c( ll, 1 ), ldc )
797  END IF
798  END IF
799 *
800 * QR iteration finished, go back and check convergence
801 *
802  GO TO 60
803 *
804 * All singular values converged, so make them positive
805 *
806  160 CONTINUE
807  DO 170 i = 1, n
808  IF( d( i ).LT.zero ) THEN
809  d( i ) = -d( i )
810 *
811 * Change sign of singular vectors, if desired
812 *
813  IF( ncvt.GT.0 )
814  $ CALL sscal( ncvt, negone, vt( i, 1 ), ldvt )
815  END IF
816  170 CONTINUE
817 *
818 * Sort the singular values into decreasing order (insertion sort on
819 * singular values, but only one transposition per singular vector)
820 *
821  DO 190 i = 1, n - 1
822 *
823 * Scan for smallest D(I)
824 *
825  isub = 1
826  smin = d( 1 )
827  DO 180 j = 2, n + 1 - i
828  IF( d( j ).LE.smin ) THEN
829  isub = j
830  smin = d( j )
831  END IF
832  180 CONTINUE
833  IF( isub.NE.n+1-i ) THEN
834 *
835 * Swap singular values and vectors
836 *
837  d( isub ) = d( n+1-i )
838  d( n+1-i ) = smin
839  IF( ncvt.GT.0 )
840  $ CALL sswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
841  $ ldvt )
842  IF( nru.GT.0 )
843  $ CALL sswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
844  IF( ncc.GT.0 )
845  $ CALL sswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
846  END IF
847  190 CONTINUE
848  GO TO 220
849 *
850 * Maximum number of iterations exceeded, failure to converge
851 *
852  200 CONTINUE
853  info = 0
854  DO 210 i = 1, n - 1
855  IF( e( i ).NE.zero )
856  $ info = info + 1
857  210 CONTINUE
858  220 CONTINUE
859  RETURN
860 *
861 * End of SBDSQR
862 *
863  END
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:138
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
sbdsqr
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:240
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:108
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
srot
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
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:199
slas2
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:107