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