1: SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
2: $ WORK, IWORK, INFO )
3: *
4: * -- LAPACK routine (version 3.2.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * June 2010
8: *
9: * .. Scalar Arguments ..
10: CHARACTER COMPQ, UPLO
11: INTEGER INFO, LDU, LDVT, N
12: * ..
13: * .. Array Arguments ..
14: INTEGER IQ( * ), IWORK( * )
15: DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ),
16: $ VT( LDVT, * ), WORK( * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * DBDSDC computes the singular value decomposition (SVD) of a real
23: * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT,
24: * using a divide and conquer method, where S is a diagonal matrix
25: * with non-negative diagonal elements (the singular values of B), and
26: * U and VT are orthogonal matrices of left and right singular vectors,
27: * respectively. DBDSDC can be used to compute all singular values,
28: * and optionally, singular vectors or singular vectors in compact form.
29: *
30: * This code makes very mild assumptions about floating point
31: * arithmetic. It will work on machines with a guard digit in
32: * add/subtract, or on those binary machines without guard digits
33: * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
34: * It could conceivably fail on hexadecimal or decimal machines
35: * without guard digits, but we know of none. See DLASD3 for details.
36: *
37: * The code currently calls DLASDQ if singular values only are desired.
38: * However, it can be slightly modified to compute singular values
39: * using the divide and conquer method.
40: *
41: * Arguments
42: * =========
43: *
44: * UPLO (input) CHARACTER*1
45: * = 'U': B is upper bidiagonal.
46: * = 'L': B is lower bidiagonal.
47: *
48: * COMPQ (input) CHARACTER*1
49: * Specifies whether singular vectors are to be computed
50: * as follows:
51: * = 'N': Compute singular values only;
52: * = 'P': Compute singular values and compute singular
53: * vectors in compact form;
54: * = 'I': Compute singular values and singular vectors.
55: *
56: * N (input) INTEGER
57: * The order of the matrix B. N >= 0.
58: *
59: * D (input/output) DOUBLE PRECISION array, dimension (N)
60: * On entry, the n diagonal elements of the bidiagonal matrix B.
61: * On exit, if INFO=0, the singular values of B.
62: *
63: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
64: * On entry, the elements of E contain the offdiagonal
65: * elements of the bidiagonal matrix whose SVD is desired.
66: * On exit, E has been destroyed.
67: *
68: * U (output) DOUBLE PRECISION array, dimension (LDU,N)
69: * If COMPQ = 'I', then:
70: * On exit, if INFO = 0, U contains the left singular vectors
71: * of the bidiagonal matrix.
72: * For other values of COMPQ, U is not referenced.
73: *
74: * LDU (input) INTEGER
75: * The leading dimension of the array U. LDU >= 1.
76: * If singular vectors are desired, then LDU >= max( 1, N ).
77: *
78: * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
79: * If COMPQ = 'I', then:
80: * On exit, if INFO = 0, VT' contains the right singular
81: * vectors of the bidiagonal matrix.
82: * For other values of COMPQ, VT is not referenced.
83: *
84: * LDVT (input) INTEGER
85: * The leading dimension of the array VT. LDVT >= 1.
86: * If singular vectors are desired, then LDVT >= max( 1, N ).
87: *
88: * Q (output) DOUBLE PRECISION array, dimension (LDQ)
89: * If COMPQ = 'P', then:
90: * On exit, if INFO = 0, Q and IQ contain the left
91: * and right singular vectors in a compact form,
92: * requiring O(N log N) space instead of 2*N**2.
93: * In particular, Q contains all the DOUBLE PRECISION data in
94: * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
95: * words of memory, where SMLSIZ is returned by ILAENV and
96: * is equal to the maximum size of the subproblems at the
97: * bottom of the computation tree (usually about 25).
98: * For other values of COMPQ, Q is not referenced.
99: *
100: * IQ (output) INTEGER array, dimension (LDIQ)
101: * If COMPQ = 'P', then:
102: * On exit, if INFO = 0, Q and IQ contain the left
103: * and right singular vectors in a compact form,
104: * requiring O(N log N) space instead of 2*N**2.
105: * In particular, IQ contains all INTEGER data in
106: * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
107: * words of memory, where SMLSIZ is returned by ILAENV and
108: * is equal to the maximum size of the subproblems at the
109: * bottom of the computation tree (usually about 25).
110: * For other values of COMPQ, IQ is not referenced.
111: *
112: * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
113: * If COMPQ = 'N' then LWORK >= (4 * N).
114: * If COMPQ = 'P' then LWORK >= (6 * N).
115: * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
116: *
117: * IWORK (workspace) INTEGER array, dimension (8*N)
118: *
119: * INFO (output) INTEGER
120: * = 0: successful exit.
121: * < 0: if INFO = -i, the i-th argument had an illegal value.
122: * > 0: The algorithm failed to compute a singular value.
123: * The update process of divide and conquer failed.
124: *
125: * Further Details
126: * ===============
127: *
128: * Based on contributions by
129: * Ming Gu and Huan Ren, Computer Science Division, University of
130: * California at Berkeley, USA
131: *
132: * =====================================================================
133: * Changed dimension statement in comment describing E from (N) to
134: * (N-1). Sven, 17 Feb 05.
135: * =====================================================================
136: *
137: * .. Parameters ..
138: DOUBLE PRECISION ZERO, ONE, TWO
139: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
140: * ..
141: * .. Local Scalars ..
142: INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
143: $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
144: $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
145: $ SMLSZP, SQRE, START, WSTART, Z
146: DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN
147: * ..
148: * .. External Functions ..
149: LOGICAL LSAME
150: INTEGER ILAENV
151: DOUBLE PRECISION DLAMCH, DLANST
152: EXTERNAL LSAME, ILAENV, DLAMCH, DLANST
153: * ..
154: * .. External Subroutines ..
155: EXTERNAL DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
156: $ DLASET, DLASR, DSWAP, XERBLA
157: * ..
158: * .. Intrinsic Functions ..
159: INTRINSIC ABS, DBLE, INT, LOG, SIGN
160: * ..
161: * .. Executable Statements ..
162: *
163: * Test the input parameters.
164: *
165: INFO = 0
166: *
167: IUPLO = 0
168: IF( LSAME( UPLO, 'U' ) )
169: $ IUPLO = 1
170: IF( LSAME( UPLO, 'L' ) )
171: $ IUPLO = 2
172: IF( LSAME( COMPQ, 'N' ) ) THEN
173: ICOMPQ = 0
174: ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
175: ICOMPQ = 1
176: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
177: ICOMPQ = 2
178: ELSE
179: ICOMPQ = -1
180: END IF
181: IF( IUPLO.EQ.0 ) THEN
182: INFO = -1
183: ELSE IF( ICOMPQ.LT.0 ) THEN
184: INFO = -2
185: ELSE IF( N.LT.0 ) THEN
186: INFO = -3
187: ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
188: $ N ) ) ) THEN
189: INFO = -7
190: ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
191: $ N ) ) ) THEN
192: INFO = -9
193: END IF
194: IF( INFO.NE.0 ) THEN
195: CALL XERBLA( 'DBDSDC', -INFO )
196: RETURN
197: END IF
198: *
199: * Quick return if possible
200: *
201: IF( N.EQ.0 )
202: $ RETURN
203: SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
204: IF( N.EQ.1 ) THEN
205: IF( ICOMPQ.EQ.1 ) THEN
206: Q( 1 ) = SIGN( ONE, D( 1 ) )
207: Q( 1+SMLSIZ*N ) = ONE
208: ELSE IF( ICOMPQ.EQ.2 ) THEN
209: U( 1, 1 ) = SIGN( ONE, D( 1 ) )
210: VT( 1, 1 ) = ONE
211: END IF
212: D( 1 ) = ABS( D( 1 ) )
213: RETURN
214: END IF
215: NM1 = N - 1
216: *
217: * If matrix lower bidiagonal, rotate to be upper bidiagonal
218: * by applying Givens rotations on the left
219: *
220: WSTART = 1
221: QSTART = 3
222: IF( ICOMPQ.EQ.1 ) THEN
223: CALL DCOPY( N, D, 1, Q( 1 ), 1 )
224: CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
225: END IF
226: IF( IUPLO.EQ.2 ) THEN
227: QSTART = 5
228: WSTART = 2*N - 1
229: DO 10 I = 1, N - 1
230: CALL DLARTG( D( I ), E( I ), CS, SN, R )
231: D( I ) = R
232: E( I ) = SN*D( I+1 )
233: D( I+1 ) = CS*D( I+1 )
234: IF( ICOMPQ.EQ.1 ) THEN
235: Q( I+2*N ) = CS
236: Q( I+3*N ) = SN
237: ELSE IF( ICOMPQ.EQ.2 ) THEN
238: WORK( I ) = CS
239: WORK( NM1+I ) = -SN
240: END IF
241: 10 CONTINUE
242: END IF
243: *
244: * If ICOMPQ = 0, use DLASDQ to compute the singular values.
245: *
246: IF( ICOMPQ.EQ.0 ) THEN
247: CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
248: $ LDU, WORK( WSTART ), INFO )
249: GO TO 40
250: END IF
251: *
252: * If N is smaller than the minimum divide size SMLSIZ, then solve
253: * the problem with another solver.
254: *
255: IF( N.LE.SMLSIZ ) THEN
256: IF( ICOMPQ.EQ.2 ) THEN
257: CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
258: CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
259: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
260: $ LDU, WORK( WSTART ), INFO )
261: ELSE IF( ICOMPQ.EQ.1 ) THEN
262: IU = 1
263: IVT = IU + N
264: CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
265: $ N )
266: CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
267: $ N )
268: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
269: $ Q( IVT+( QSTART-1 )*N ), N,
270: $ Q( IU+( QSTART-1 )*N ), N,
271: $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
272: $ INFO )
273: END IF
274: GO TO 40
275: END IF
276: *
277: IF( ICOMPQ.EQ.2 ) THEN
278: CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
279: CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
280: END IF
281: *
282: * Scale.
283: *
284: ORGNRM = DLANST( 'M', N, D, E )
285: IF( ORGNRM.EQ.ZERO )
286: $ RETURN
287: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
288: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
289: *
290: EPS = DLAMCH( 'Epsilon' )
291: *
292: MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
293: SMLSZP = SMLSIZ + 1
294: *
295: IF( ICOMPQ.EQ.1 ) THEN
296: IU = 1
297: IVT = 1 + SMLSIZ
298: DIFL = IVT + SMLSZP
299: DIFR = DIFL + MLVL
300: Z = DIFR + MLVL*2
301: IC = Z + MLVL
302: IS = IC + 1
303: POLES = IS + 1
304: GIVNUM = POLES + 2*MLVL
305: *
306: K = 1
307: GIVPTR = 2
308: PERM = 3
309: GIVCOL = PERM + MLVL
310: END IF
311: *
312: DO 20 I = 1, N
313: IF( ABS( D( I ) ).LT.EPS ) THEN
314: D( I ) = SIGN( EPS, D( I ) )
315: END IF
316: 20 CONTINUE
317: *
318: START = 1
319: SQRE = 0
320: *
321: DO 30 I = 1, NM1
322: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
323: *
324: * Subproblem found. First determine its size and then
325: * apply divide and conquer on it.
326: *
327: IF( I.LT.NM1 ) THEN
328: *
329: * A subproblem with E(I) small for I < NM1.
330: *
331: NSIZE = I - START + 1
332: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
333: *
334: * A subproblem with E(NM1) not too small but I = NM1.
335: *
336: NSIZE = N - START + 1
337: ELSE
338: *
339: * A subproblem with E(NM1) small. This implies an
340: * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem
341: * first.
342: *
343: NSIZE = I - START + 1
344: IF( ICOMPQ.EQ.2 ) THEN
345: U( N, N ) = SIGN( ONE, D( N ) )
346: VT( N, N ) = ONE
347: ELSE IF( ICOMPQ.EQ.1 ) THEN
348: Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
349: Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
350: END IF
351: D( N ) = ABS( D( N ) )
352: END IF
353: IF( ICOMPQ.EQ.2 ) THEN
354: CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
355: $ U( START, START ), LDU, VT( START, START ),
356: $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
357: ELSE
358: CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
359: $ E( START ), Q( START+( IU+QSTART-2 )*N ), N,
360: $ Q( START+( IVT+QSTART-2 )*N ),
361: $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
362: $ N ), Q( START+( DIFR+QSTART-2 )*N ),
363: $ Q( START+( Z+QSTART-2 )*N ),
364: $ Q( START+( POLES+QSTART-2 )*N ),
365: $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
366: $ N, IQ( START+PERM*N ),
367: $ Q( START+( GIVNUM+QSTART-2 )*N ),
368: $ Q( START+( IC+QSTART-2 )*N ),
369: $ Q( START+( IS+QSTART-2 )*N ),
370: $ WORK( WSTART ), IWORK, INFO )
371: END IF
372: IF( INFO.NE.0 ) THEN
373: RETURN
374: END IF
375: START = I + 1
376: END IF
377: 30 CONTINUE
378: *
379: * Unscale
380: *
381: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
382: 40 CONTINUE
383: *
384: * Use Selection Sort to minimize swaps of singular vectors
385: *
386: DO 60 II = 2, N
387: I = II - 1
388: KK = I
389: P = D( I )
390: DO 50 J = II, N
391: IF( D( J ).GT.P ) THEN
392: KK = J
393: P = D( J )
394: END IF
395: 50 CONTINUE
396: IF( KK.NE.I ) THEN
397: D( KK ) = D( I )
398: D( I ) = P
399: IF( ICOMPQ.EQ.1 ) THEN
400: IQ( I ) = KK
401: ELSE IF( ICOMPQ.EQ.2 ) THEN
402: CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
403: CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
404: END IF
405: ELSE IF( ICOMPQ.EQ.1 ) THEN
406: IQ( I ) = I
407: END IF
408: 60 CONTINUE
409: *
410: * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
411: *
412: IF( ICOMPQ.EQ.1 ) THEN
413: IF( IUPLO.EQ.1 ) THEN
414: IQ( N ) = 1
415: ELSE
416: IQ( N ) = 0
417: END IF
418: END IF
419: *
420: * If B is lower bidiagonal, update U by those Givens rotations
421: * which rotated B to be upper bidiagonal
422: *
423: IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
424: $ CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
425: *
426: RETURN
427: *
428: * End of DBDSDC
429: *
430: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dbdsdc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack