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