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