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