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