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