KTH framework for Nek5000 toolboxes; testing version  0.0.1
dbdsqr.f
Go to the documentation of this file.
1  SUBROUTINE dbdsqr( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
2  $ LDU, C, LDC, WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * October 31, 1999
8 *
9 * .. Scalar Arguments ..
10  CHARACTER UPLO
11  INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
12 * ..
13 * .. Array Arguments ..
14  DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
15  $ vt( ldvt, * ), work( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DBDSQR computes the singular value decomposition (SVD) of a real
22 * N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P'
23 * denotes the transpose of P), where S is a diagonal matrix with
24 * non-negative diagonal elements (the singular values of B), and Q
25 * and P are orthogonal matrices.
26 *
27 * The routine computes S, and optionally computes U * Q, P' * VT,
28 * or Q' * C, for given real input matrices U, VT, and C.
29 *
30 * See "Computing Small Singular Values of Bidiagonal Matrices With
31 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
32 * LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
33 * no. 5, pp. 873-912, Sept 1990) and
34 * "Accurate singular values and differential qd algorithms," by
35 * B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
36 * Department, University of California at Berkeley, July 1992
37 * for a detailed description of the algorithm.
38 *
39 * Arguments
40 * =========
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': B is upper bidiagonal;
44 * = 'L': B is lower bidiagonal.
45 *
46 * N (input) INTEGER
47 * The order of the matrix B. N >= 0.
48 *
49 * NCVT (input) INTEGER
50 * The number of columns of the matrix VT. NCVT >= 0.
51 *
52 * NRU (input) INTEGER
53 * The number of rows of the matrix U. NRU >= 0.
54 *
55 * NCC (input) INTEGER
56 * The number of columns of the matrix C. NCC >= 0.
57 *
58 * D (input/output) DOUBLE PRECISION array, dimension (N)
59 * On entry, the n diagonal elements of the bidiagonal matrix B.
60 * On exit, if INFO=0, the singular values of B in decreasing
61 * order.
62 *
63 * E (input/output) DOUBLE PRECISION array, dimension (N)
64 * On entry, the elements of E contain the
65 * offdiagonal elements of the bidiagonal matrix whose SVD
66 * is desired. On normal exit (INFO = 0), E is destroyed.
67 * If the algorithm does not converge (INFO > 0), D and E
68 * will contain the diagonal and superdiagonal elements of a
69 * bidiagonal matrix orthogonally equivalent to the one given
70 * as input. E(N) is used for workspace.
71 *
72 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
73 * On entry, an N-by-NCVT matrix VT.
74 * On exit, VT is overwritten by P' * VT.
75 * VT is not referenced if NCVT = 0.
76 *
77 * LDVT (input) INTEGER
78 * The leading dimension of the array VT.
79 * LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
80 *
81 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
82 * On entry, an NRU-by-N matrix U.
83 * On exit, U is overwritten by U * Q.
84 * U is not referenced if NRU = 0.
85 *
86 * LDU (input) INTEGER
87 * The leading dimension of the array U. LDU >= max(1,NRU).
88 *
89 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
90 * On entry, an N-by-NCC matrix C.
91 * On exit, C is overwritten by Q' * C.
92 * C is not referenced if NCC = 0.
93 *
94 * LDC (input) INTEGER
95 * The leading dimension of the array C.
96 * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
97 *
98 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
99 *
100 * INFO (output) INTEGER
101 * = 0: successful exit
102 * < 0: If INFO = -i, the i-th argument had an illegal value
103 * > 0: the algorithm did not converge; D and E contain the
104 * elements of a bidiagonal matrix which is orthogonally
105 * similar to the input matrix B; if INFO = i, i
106 * elements of E have not converged to zero.
107 *
108 * Internal Parameters
109 * ===================
110 *
111 * TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
112 * TOLMUL controls the convergence criterion of the QR loop.
113 * If it is positive, TOLMUL*EPS is the desired relative
114 * precision in the computed singular values.
115 * If it is negative, abs(TOLMUL*EPS*sigma_max) is the
116 * desired absolute accuracy in the computed singular
117 * values (corresponds to relative accuracy
118 * abs(TOLMUL*EPS) in the largest singular value.
119 * abs(TOLMUL) should be between 1 and 1/EPS, and preferably
120 * between 10 (for fast convergence) and .1/EPS
121 * (for there to be some accuracy in the results).
122 * Default is to lose at either one eighth or 2 of the
123 * available decimal digits in each computed singular value
124 * (whichever is smaller).
125 *
126 * MAXITR INTEGER, default = 6
127 * MAXITR controls the maximum number of passes of the
128 * algorithm through its inner loop. The algorithms stops
129 * (and so fails to converge) if the number of passes
130 * through the inner loop exceeds MAXITR*N**2.
131 *
132 * =====================================================================
133 *
134 * .. Parameters ..
135  DOUBLE PRECISION ZERO
136  parameter( zero = 0.0d0 )
137  DOUBLE PRECISION ONE
138  parameter( one = 1.0d0 )
139  DOUBLE PRECISION NEGONE
140  parameter( negone = -1.0d0 )
141  DOUBLE PRECISION HNDRTH
142  parameter( hndrth = 0.01d0 )
143  DOUBLE PRECISION TEN
144  parameter( ten = 10.0d0 )
145  DOUBLE PRECISION HNDRD
146  parameter( hndrd = 100.0d0 )
147  DOUBLE PRECISION MEIGTH
148  parameter( meigth = -0.125d0 )
149  INTEGER MAXITR
150  parameter( maxitr = 6 )
151 * ..
152 * .. Local Scalars ..
153  LOGICAL LOWER, ROTATE
154  INTEGER I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,
155  $ nm12, nm13, oldll, oldm
156  DOUBLE PRECISION ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
157  $ oldcs, oldsn, r, shift, sigmn, sigmx, sinl,
158  $ sinr, sll, smax, smin, sminl, sminlo, sminoa,
159  $ sn, thresh, tol, tolmul, unfl
160 * ..
161 * .. External Functions ..
162  LOGICAL LSAME
163  DOUBLE PRECISION DLAMCH
164  EXTERNAL lsame, dlamch
165 * ..
166 * .. External Subroutines ..
167  EXTERNAL dlartg, dlas2, dlasq1, dlasr, dlasv2, drot,
168  $ dscal, dswap, xerbla
169 * ..
170 * .. Intrinsic Functions ..
171  INTRINSIC abs, dble, max, min, sign, sqrt
172 * ..
173 * .. Executable Statements ..
174 *
175 * Test the input parameters.
176 *
177  info = 0
178  lower = lsame( uplo, 'L' )
179  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
180  info = -1
181  ELSE IF( n.LT.0 ) THEN
182  info = -2
183  ELSE IF( ncvt.LT.0 ) THEN
184  info = -3
185  ELSE IF( nru.LT.0 ) THEN
186  info = -4
187  ELSE IF( ncc.LT.0 ) THEN
188  info = -5
189  ELSE IF( ( ncvt.EQ.0 .AND. ldvt.LT.1 ) .OR.
190  $ ( ncvt.GT.0 .AND. ldvt.LT.max( 1, n ) ) ) THEN
191  info = -9
192  ELSE IF( ldu.LT.max( 1, nru ) ) THEN
193  info = -11
194  ELSE IF( ( ncc.EQ.0 .AND. ldc.LT.1 ) .OR.
195  $ ( ncc.GT.0 .AND. ldc.LT.max( 1, n ) ) ) THEN
196  info = -13
197  END IF
198  IF( info.NE.0 ) THEN
199  CALL xerbla( 'DBDSQR', -info )
200  RETURN
201  END IF
202  IF( n.EQ.0 )
203  $ RETURN
204  IF( n.EQ.1 )
205  $ GO TO 160
206 *
207 * ROTATE is true if any singular vectors desired, false otherwise
208 *
209  rotate = ( ncvt.GT.0 ) .OR. ( nru.GT.0 ) .OR. ( ncc.GT.0 )
210 *
211 * If no singular vectors desired, use qd algorithm
212 *
213  IF( .NOT.rotate ) THEN
214  CALL dlasq1( n, d, e, work, info )
215  RETURN
216  END IF
217 *
218  nm1 = n - 1
219  nm12 = nm1 + nm1
220  nm13 = nm12 + nm1
221  idir = 0
222 *
223 * Get machine constants
224 *
225  eps = dlamch( 'Epsilon' )
226  unfl = dlamch( 'Safe minimum' )
227 *
228 * If matrix lower bidiagonal, rotate to be upper bidiagonal
229 * by applying Givens rotations on the left
230 *
231  IF( lower ) THEN
232  DO 10 i = 1, n - 1
233  CALL dlartg( d( i ), e( i ), cs, sn, r )
234  d( i ) = r
235  e( i ) = sn*d( i+1 )
236  d( i+1 ) = cs*d( i+1 )
237  work( i ) = cs
238  work( nm1+i ) = sn
239  10 CONTINUE
240 *
241 * Update singular vectors if desired
242 *
243  IF( nru.GT.0 )
244  $ CALL dlasr( 'R', 'V', 'F', nru, n, work( 1 ), work( n ), u,
245  $ ldu )
246  IF( ncc.GT.0 )
247  $ CALL dlasr( 'L', 'V', 'F', n, ncc, work( 1 ), work( n ), c,
248  $ ldc )
249  END IF
250 *
251 * Compute singular values to relative accuracy TOL
252 * (By setting TOL to be negative, algorithm will compute
253 * singular values to absolute accuracy ABS(TOL)*norm(input matrix))
254 *
255  tolmul = max( ten, min( hndrd, eps**meigth ) )
256  tol = tolmul*eps
257 *
258 * Compute approximate maximum, minimum singular values
259 *
260  smax = zero
261  DO 20 i = 1, n
262  smax = max( smax, abs( d( i ) ) )
263  20 CONTINUE
264  DO 30 i = 1, n - 1
265  smax = max( smax, abs( e( i ) ) )
266  30 CONTINUE
267  sminl = zero
268  IF( tol.GE.zero ) THEN
269 *
270 * Relative accuracy desired
271 *
272  sminoa = abs( d( 1 ) )
273  IF( sminoa.EQ.zero )
274  $ GO TO 50
275  mu = sminoa
276  DO 40 i = 2, n
277  mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
278  sminoa = min( sminoa, mu )
279  IF( sminoa.EQ.zero )
280  $ GO TO 50
281  40 CONTINUE
282  50 CONTINUE
283  sminoa = sminoa / sqrt( dble( n ) )
284  thresh = max( tol*sminoa, maxitr*n*n*unfl )
285  ELSE
286 *
287 * Absolute accuracy desired
288 *
289  thresh = max( abs( tol )*smax, maxitr*n*n*unfl )
290  END IF
291 *
292 * Prepare for main iteration loop for the singular values
293 * (MAXIT is the maximum number of passes through the inner
294 * loop permitted before nonconvergence signalled.)
295 *
296  maxit = maxitr*n*n
297  iter = 0
298  oldll = -1
299  oldm = -1
300 *
301 * M points to last element of unconverged part of matrix
302 *
303  m = n
304 *
305 * Begin main iteration loop
306 *
307  60 CONTINUE
308 *
309 * Check for convergence or exceeding iteration count
310 *
311  IF( m.LE.1 )
312  $ GO TO 160
313  IF( iter.GT.maxit )
314  $ GO TO 200
315 *
316 * Find diagonal block of matrix to work on
317 *
318  IF( tol.LT.zero .AND. abs( d( m ) ).LE.thresh )
319  $ d( m ) = zero
320  smax = abs( d( m ) )
321  smin = smax
322  DO 70 lll = 1, m - 1
323  ll = m - lll
324  abss = abs( d( ll ) )
325  abse = abs( e( ll ) )
326  IF( tol.LT.zero .AND. abss.LE.thresh )
327  $ d( ll ) = zero
328  IF( abse.LE.thresh )
329  $ GO TO 80
330  smin = min( smin, abss )
331  smax = max( smax, abss, abse )
332  70 CONTINUE
333  ll = 0
334  GO TO 90
335  80 CONTINUE
336  e( ll ) = zero
337 *
338 * Matrix splits since E(LL) = 0
339 *
340  IF( ll.EQ.m-1 ) THEN
341 *
342 * Convergence of bottom singular value, return to top of loop
343 *
344  m = m - 1
345  GO TO 60
346  END IF
347  90 CONTINUE
348  ll = ll + 1
349 *
350 * E(LL) through E(M-1) are nonzero, E(LL-1) is zero
351 *
352  IF( ll.EQ.m-1 ) THEN
353 *
354 * 2 by 2 block, handle separately
355 *
356  CALL dlasv2( d( m-1 ), e( m-1 ), d( m ), sigmn, sigmx, sinr,
357  $ cosr, sinl, cosl )
358  d( m-1 ) = sigmx
359  e( m-1 ) = zero
360  d( m ) = sigmn
361 *
362 * Compute singular vectors, if desired
363 *
364  IF( ncvt.GT.0 )
365  $ CALL drot( ncvt, vt( m-1, 1 ), ldvt, vt( m, 1 ), ldvt, cosr,
366  $ sinr )
367  IF( nru.GT.0 )
368  $ CALL drot( nru, u( 1, m-1 ), 1, u( 1, m ), 1, cosl, sinl )
369  IF( ncc.GT.0 )
370  $ CALL drot( ncc, c( m-1, 1 ), ldc, c( m, 1 ), ldc, cosl,
371  $ sinl )
372  m = m - 2
373  GO TO 60
374  END IF
375 *
376 * If working on new submatrix, choose shift direction
377 * (from larger end diagonal element towards smaller)
378 *
379  IF( ll.GT.oldm .OR. m.LT.oldll ) THEN
380  IF( abs( d( ll ) ).GE.abs( d( m ) ) ) THEN
381 *
382 * Chase bulge from top (big end) to bottom (small end)
383 *
384  idir = 1
385  ELSE
386 *
387 * Chase bulge from bottom (big end) to top (small end)
388 *
389  idir = 2
390  END IF
391  END IF
392 *
393 * Apply convergence tests
394 *
395  IF( idir.EQ.1 ) THEN
396 *
397 * Run convergence test in forward direction
398 * First apply standard test to bottom of matrix
399 *
400  IF( abs( e( m-1 ) ).LE.abs( tol )*abs( d( m ) ) .OR.
401  $ ( tol.LT.zero .AND. abs( e( m-1 ) ).LE.thresh ) ) THEN
402  e( m-1 ) = zero
403  GO TO 60
404  END IF
405 *
406  IF( tol.GE.zero ) THEN
407 *
408 * If relative accuracy desired,
409 * apply convergence criterion forward
410 *
411  mu = abs( d( ll ) )
412  sminl = mu
413  DO 100 lll = ll, m - 1
414  IF( abs( e( lll ) ).LE.tol*mu ) THEN
415  e( lll ) = zero
416  GO TO 60
417  END IF
418  sminlo = sminl
419  mu = abs( d( lll+1 ) )*( mu / ( mu+abs( e( lll ) ) ) )
420  sminl = min( sminl, mu )
421  100 CONTINUE
422  END IF
423 *
424  ELSE
425 *
426 * Run convergence test in backward direction
427 * First apply standard test to top of matrix
428 *
429  IF( abs( e( ll ) ).LE.abs( tol )*abs( d( ll ) ) .OR.
430  $ ( tol.LT.zero .AND. abs( e( ll ) ).LE.thresh ) ) THEN
431  e( ll ) = zero
432  GO TO 60
433  END IF
434 *
435  IF( tol.GE.zero ) THEN
436 *
437 * If relative accuracy desired,
438 * apply convergence criterion backward
439 *
440  mu = abs( d( m ) )
441  sminl = mu
442  DO 110 lll = m - 1, ll, -1
443  IF( abs( e( lll ) ).LE.tol*mu ) THEN
444  e( lll ) = zero
445  GO TO 60
446  END IF
447  sminlo = sminl
448  mu = abs( d( lll ) )*( mu / ( mu+abs( e( lll ) ) ) )
449  sminl = min( sminl, mu )
450  110 CONTINUE
451  END IF
452  END IF
453  oldll = ll
454  oldm = m
455 *
456 * Compute shift. First, test if shifting would ruin relative
457 * accuracy, and if so set the shift to zero.
458 *
459  IF( tol.GE.zero .AND. n*tol*( sminl / smax ).LE.
460  $ max( eps, hndrth*tol ) ) THEN
461 *
462 * Use a zero shift to avoid loss of relative accuracy
463 *
464  shift = zero
465  ELSE
466 *
467 * Compute the shift from 2-by-2 block at end of matrix
468 *
469  IF( idir.EQ.1 ) THEN
470  sll = abs( d( ll ) )
471  CALL dlas2( d( m-1 ), e( m-1 ), d( m ), shift, r )
472  ELSE
473  sll = abs( d( m ) )
474  CALL dlas2( d( ll ), e( ll ), d( ll+1 ), shift, r )
475  END IF
476 *
477 * Test if shift negligible, and if so set to zero
478 *
479  IF( sll.GT.zero ) THEN
480  IF( ( shift / sll )**2.LT.eps )
481  $ shift = zero
482  END IF
483  END IF
484 *
485 * Increment iteration count
486 *
487  iter = iter + m - ll
488 *
489 * If SHIFT = 0, do simplified QR iteration
490 *
491  IF( shift.EQ.zero ) THEN
492  IF( idir.EQ.1 ) THEN
493 *
494 * Chase bulge from top to bottom
495 * Save cosines and sines for later singular vector updates
496 *
497  cs = one
498  oldcs = one
499  DO 120 i = ll, m - 1
500  CALL dlartg( d( i )*cs, e( i ), cs, sn, r )
501  IF( i.GT.ll )
502  $ e( i-1 ) = oldsn*r
503  CALL dlartg( oldcs*r, d( i+1 )*sn, oldcs, oldsn, d( i ) )
504  work( i-ll+1 ) = cs
505  work( i-ll+1+nm1 ) = sn
506  work( i-ll+1+nm12 ) = oldcs
507  work( i-ll+1+nm13 ) = oldsn
508  120 CONTINUE
509  h = d( m )*cs
510  d( m ) = h*oldcs
511  e( m-1 ) = h*oldsn
512 *
513 * Update singular vectors
514 *
515  IF( ncvt.GT.0 )
516  $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
517  $ work( n ), vt( ll, 1 ), ldvt )
518  IF( nru.GT.0 )
519  $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
520  $ work( nm13+1 ), u( 1, ll ), ldu )
521  IF( ncc.GT.0 )
522  $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
523  $ work( nm13+1 ), c( ll, 1 ), ldc )
524 *
525 * Test convergence
526 *
527  IF( abs( e( m-1 ) ).LE.thresh )
528  $ e( m-1 ) = zero
529 *
530  ELSE
531 *
532 * Chase bulge from bottom to top
533 * Save cosines and sines for later singular vector updates
534 *
535  cs = one
536  oldcs = one
537  DO 130 i = m, ll + 1, -1
538  CALL dlartg( d( i )*cs, e( i-1 ), cs, sn, r )
539  IF( i.LT.m )
540  $ e( i ) = oldsn*r
541  CALL dlartg( oldcs*r, d( i-1 )*sn, oldcs, oldsn, d( i ) )
542  work( i-ll ) = cs
543  work( i-ll+nm1 ) = -sn
544  work( i-ll+nm12 ) = oldcs
545  work( i-ll+nm13 ) = -oldsn
546  130 CONTINUE
547  h = d( ll )*cs
548  d( ll ) = h*oldcs
549  e( ll ) = h*oldsn
550 *
551 * Update singular vectors
552 *
553  IF( ncvt.GT.0 )
554  $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
555  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
556  IF( nru.GT.0 )
557  $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
558  $ work( n ), u( 1, ll ), ldu )
559  IF( ncc.GT.0 )
560  $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
561  $ work( n ), c( ll, 1 ), ldc )
562 *
563 * Test convergence
564 *
565  IF( abs( e( ll ) ).LE.thresh )
566  $ e( ll ) = zero
567  END IF
568  ELSE
569 *
570 * Use nonzero shift
571 *
572  IF( idir.EQ.1 ) THEN
573 *
574 * Chase bulge from top to bottom
575 * Save cosines and sines for later singular vector updates
576 *
577  f = ( abs( d( ll ) )-shift )*
578  $ ( sign( one, d( ll ) )+shift / d( ll ) )
579  g = e( ll )
580  DO 140 i = ll, m - 1
581  CALL dlartg( f, g, cosr, sinr, r )
582  IF( i.GT.ll )
583  $ e( i-1 ) = r
584  f = cosr*d( i ) + sinr*e( i )
585  e( i ) = cosr*e( i ) - sinr*d( i )
586  g = sinr*d( i+1 )
587  d( i+1 ) = cosr*d( i+1 )
588  CALL dlartg( f, g, cosl, sinl, r )
589  d( i ) = r
590  f = cosl*e( i ) + sinl*d( i+1 )
591  d( i+1 ) = cosl*d( i+1 ) - sinl*e( i )
592  IF( i.LT.m-1 ) THEN
593  g = sinl*e( i+1 )
594  e( i+1 ) = cosl*e( i+1 )
595  END IF
596  work( i-ll+1 ) = cosr
597  work( i-ll+1+nm1 ) = sinr
598  work( i-ll+1+nm12 ) = cosl
599  work( i-ll+1+nm13 ) = sinl
600  140 CONTINUE
601  e( m-1 ) = f
602 *
603 * Update singular vectors
604 *
605  IF( ncvt.GT.0 )
606  $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncvt, work( 1 ),
607  $ work( n ), vt( ll, 1 ), ldvt )
608  IF( nru.GT.0 )
609  $ CALL dlasr( 'R', 'V', 'F', nru, m-ll+1, work( nm12+1 ),
610  $ work( nm13+1 ), u( 1, ll ), ldu )
611  IF( ncc.GT.0 )
612  $ CALL dlasr( 'L', 'V', 'F', m-ll+1, ncc, work( nm12+1 ),
613  $ work( nm13+1 ), c( ll, 1 ), ldc )
614 *
615 * Test convergence
616 *
617  IF( abs( e( m-1 ) ).LE.thresh )
618  $ e( m-1 ) = zero
619 *
620  ELSE
621 *
622 * Chase bulge from bottom to top
623 * Save cosines and sines for later singular vector updates
624 *
625  f = ( abs( d( m ) )-shift )*( sign( one, d( m ) )+shift /
626  $ d( m ) )
627  g = e( m-1 )
628  DO 150 i = m, ll + 1, -1
629  CALL dlartg( f, g, cosr, sinr, r )
630  IF( i.LT.m )
631  $ e( i ) = r
632  f = cosr*d( i ) + sinr*e( i-1 )
633  e( i-1 ) = cosr*e( i-1 ) - sinr*d( i )
634  g = sinr*d( i-1 )
635  d( i-1 ) = cosr*d( i-1 )
636  CALL dlartg( f, g, cosl, sinl, r )
637  d( i ) = r
638  f = cosl*e( i-1 ) + sinl*d( i-1 )
639  d( i-1 ) = cosl*d( i-1 ) - sinl*e( i-1 )
640  IF( i.GT.ll+1 ) THEN
641  g = sinl*e( i-2 )
642  e( i-2 ) = cosl*e( i-2 )
643  END IF
644  work( i-ll ) = cosr
645  work( i-ll+nm1 ) = -sinr
646  work( i-ll+nm12 ) = cosl
647  work( i-ll+nm13 ) = -sinl
648  150 CONTINUE
649  e( ll ) = f
650 *
651 * Test convergence
652 *
653  IF( abs( e( ll ) ).LE.thresh )
654  $ e( ll ) = zero
655 *
656 * Update singular vectors if desired
657 *
658  IF( ncvt.GT.0 )
659  $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncvt, work( nm12+1 ),
660  $ work( nm13+1 ), vt( ll, 1 ), ldvt )
661  IF( nru.GT.0 )
662  $ CALL dlasr( 'R', 'V', 'B', nru, m-ll+1, work( 1 ),
663  $ work( n ), u( 1, ll ), ldu )
664  IF( ncc.GT.0 )
665  $ CALL dlasr( 'L', 'V', 'B', m-ll+1, ncc, work( 1 ),
666  $ work( n ), c( ll, 1 ), ldc )
667  END IF
668  END IF
669 *
670 * QR iteration finished, go back and check convergence
671 *
672  GO TO 60
673 *
674 * All singular values converged, so make them positive
675 *
676  160 CONTINUE
677  DO 170 i = 1, n
678  IF( d( i ).LT.zero ) THEN
679  d( i ) = -d( i )
680 *
681 * Change sign of singular vectors, if desired
682 *
683  IF( ncvt.GT.0 )
684  $ CALL dscal( ncvt, negone, vt( i, 1 ), ldvt )
685  END IF
686  170 CONTINUE
687 *
688 * Sort the singular values into decreasing order (insertion sort on
689 * singular values, but only one transposition per singular vector)
690 *
691  DO 190 i = 1, n - 1
692 *
693 * Scan for smallest D(I)
694 *
695  isub = 1
696  smin = d( 1 )
697  DO 180 j = 2, n + 1 - i
698  IF( d( j ).LE.smin ) THEN
699  isub = j
700  smin = d( j )
701  END IF
702  180 CONTINUE
703  IF( isub.NE.n+1-i ) THEN
704 *
705 * Swap singular values and vectors
706 *
707  d( isub ) = d( n+1-i )
708  d( n+1-i ) = smin
709  IF( ncvt.GT.0 )
710  $ CALL dswap( ncvt, vt( isub, 1 ), ldvt, vt( n+1-i, 1 ),
711  $ ldvt )
712  IF( nru.GT.0 )
713  $ CALL dswap( nru, u( 1, isub ), 1, u( 1, n+1-i ), 1 )
714  IF( ncc.GT.0 )
715  $ CALL dswap( ncc, c( isub, 1 ), ldc, c( n+1-i, 1 ), ldc )
716  END IF
717  190 CONTINUE
718  GO TO 220
719 *
720 * Maximum number of iterations exceeded, failure to converge
721 *
722  200 CONTINUE
723  info = 0
724  DO 210 i = 1, n - 1
725  IF( e( i ).NE.zero )
726  $ info = info + 1
727  210 CONTINUE
728  220 CONTINUE
729  RETURN
730 *
731 * End of DBDSQR
732 *
733  END
dbdsqr
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
Definition: dbdsqr.f:3
dlartg
subroutine dlartg(F, G, CS, SN, R)
Definition: dlartg.f:2
dlas2
subroutine dlas2(F, G, H, SSMIN, SSMAX)
Definition: dlas2.f:2
dlasq1
subroutine dlasq1(N, D, E, WORK, INFO)
Definition: dlasq1.f:2
dlasr
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
Definition: dlasr.f:2
dlasv2
subroutine dlasv2(F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL)
Definition: dlasv2.f:2
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
Definition: drot.f:2
dscal
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2
dswap
subroutine dswap(n, dx, incx, dy, incy)
Definition: dswap.f:2
xerbla
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2