Annotation of rpl/lapack/lapack/dlasdq.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
2: $ U, LDU, C, LDC, WORK, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.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: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER UPLO
11: INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
15: $ VT( LDVT, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLASDQ computes the singular value decomposition (SVD) of a real
22: * (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
23: * E, accumulating the transformations if desired. Letting B denote
24: * the input bidiagonal matrix, the algorithm computes orthogonal
25: * matrices Q and P such that B = Q * S * P' (P' denotes the transpose
26: * of P). The singular values S are overwritten on D.
27: *
28: * The input matrix U is changed to U * Q if desired.
29: * The input matrix VT is changed to P' * VT if desired.
30: * The input matrix C is changed to Q' * C if desired.
31: *
32: * See "Computing Small Singular Values of Bidiagonal Matrices With
33: * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
34: * LAPACK Working Note #3, for a detailed description of the algorithm.
35: *
36: * Arguments
37: * =========
38: *
39: * UPLO (input) CHARACTER*1
40: * On entry, UPLO specifies whether the input bidiagonal matrix
41: * is upper or lower bidiagonal, and wether it is square are
42: * not.
43: * UPLO = 'U' or 'u' B is upper bidiagonal.
44: * UPLO = 'L' or 'l' B is lower bidiagonal.
45: *
46: * SQRE (input) INTEGER
47: * = 0: then the input matrix is N-by-N.
48: * = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
49: * (N+1)-by-N if UPLU = 'L'.
50: *
51: * The bidiagonal matrix has
52: * N = NL + NR + 1 rows and
53: * M = N + SQRE >= N columns.
54: *
55: * N (input) INTEGER
56: * On entry, N specifies the number of rows and columns
57: * in the matrix. N must be at least 0.
58: *
59: * NCVT (input) INTEGER
60: * On entry, NCVT specifies the number of columns of
61: * the matrix VT. NCVT must be at least 0.
62: *
63: * NRU (input) INTEGER
64: * On entry, NRU specifies the number of rows of
65: * the matrix U. NRU must be at least 0.
66: *
67: * NCC (input) INTEGER
68: * On entry, NCC specifies the number of columns of
69: * the matrix C. NCC must be at least 0.
70: *
71: * D (input/output) DOUBLE PRECISION array, dimension (N)
72: * On entry, D contains the diagonal entries of the
73: * bidiagonal matrix whose SVD is desired. On normal exit,
74: * D contains the singular values in ascending order.
75: *
76: * E (input/output) DOUBLE PRECISION array.
77: * dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
78: * On entry, the entries of E contain the offdiagonal entries
79: * of the bidiagonal matrix whose SVD is desired. On normal
80: * exit, E will contain 0. If the algorithm does not converge,
81: * D and E will contain the diagonal and superdiagonal entries
82: * of a bidiagonal matrix orthogonally equivalent to the one
83: * given as input.
84: *
85: * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)
86: * On entry, contains a matrix which on exit has been
87: * premultiplied by P', dimension N-by-NCVT if SQRE = 0
88: * and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
89: *
90: * LDVT (input) INTEGER
91: * On entry, LDVT specifies the leading dimension of VT as
92: * declared in the calling (sub) program. LDVT must be at
93: * least 1. If NCVT is nonzero LDVT must also be at least N.
94: *
95: * U (input/output) DOUBLE PRECISION array, dimension (LDU, N)
96: * On entry, contains a matrix which on exit has been
97: * postmultiplied by Q, dimension NRU-by-N if SQRE = 0
98: * and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
99: *
100: * LDU (input) INTEGER
101: * On entry, LDU specifies the leading dimension of U as
102: * declared in the calling (sub) program. LDU must be at
103: * least max( 1, NRU ) .
104: *
105: * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)
106: * On entry, contains an N-by-NCC matrix which on exit
107: * has been premultiplied by Q' dimension N-by-NCC if SQRE = 0
108: * and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
109: *
110: * LDC (input) INTEGER
111: * On entry, LDC specifies the leading dimension of C as
112: * declared in the calling (sub) program. LDC must be at
113: * least 1. If NCC is nonzero, LDC must also be at least N.
114: *
115: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
116: * Workspace. Only referenced if one of NCVT, NRU, or NCC is
117: * nonzero, and if N is at least 2.
118: *
119: * INFO (output) INTEGER
120: * On exit, a value of 0 indicates a successful exit.
121: * If INFO < 0, argument number -INFO is illegal.
122: * If INFO > 0, the algorithm did not converge, and INFO
123: * specifies how many superdiagonals did not converge.
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: *
134: * .. Parameters ..
135: DOUBLE PRECISION ZERO
136: PARAMETER ( ZERO = 0.0D+0 )
137: * ..
138: * .. Local Scalars ..
139: LOGICAL ROTATE
140: INTEGER I, ISUB, IUPLO, J, NP1, SQRE1
141: DOUBLE PRECISION CS, R, SMIN, SN
142: * ..
143: * .. External Subroutines ..
144: EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
145: * ..
146: * .. External Functions ..
147: LOGICAL LSAME
148: EXTERNAL LSAME
149: * ..
150: * .. Intrinsic Functions ..
151: INTRINSIC MAX
152: * ..
153: * .. Executable Statements ..
154: *
155: * Test the input parameters.
156: *
157: INFO = 0
158: IUPLO = 0
159: IF( LSAME( UPLO, 'U' ) )
160: $ IUPLO = 1
161: IF( LSAME( UPLO, 'L' ) )
162: $ IUPLO = 2
163: IF( IUPLO.EQ.0 ) THEN
164: INFO = -1
165: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
166: INFO = -2
167: ELSE IF( N.LT.0 ) THEN
168: INFO = -3
169: ELSE IF( NCVT.LT.0 ) THEN
170: INFO = -4
171: ELSE IF( NRU.LT.0 ) THEN
172: INFO = -5
173: ELSE IF( NCC.LT.0 ) THEN
174: INFO = -6
175: ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
176: $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
177: INFO = -10
178: ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
179: INFO = -12
180: ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
181: $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
182: INFO = -14
183: END IF
184: IF( INFO.NE.0 ) THEN
185: CALL XERBLA( 'DLASDQ', -INFO )
186: RETURN
187: END IF
188: IF( N.EQ.0 )
189: $ RETURN
190: *
191: * ROTATE is true if any singular vectors desired, false otherwise
192: *
193: ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
194: NP1 = N + 1
195: SQRE1 = SQRE
196: *
197: * If matrix non-square upper bidiagonal, rotate to be lower
198: * bidiagonal. The rotations are on the right.
199: *
200: IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
201: DO 10 I = 1, N - 1
202: CALL DLARTG( D( I ), E( I ), CS, SN, R )
203: D( I ) = R
204: E( I ) = SN*D( I+1 )
205: D( I+1 ) = CS*D( I+1 )
206: IF( ROTATE ) THEN
207: WORK( I ) = CS
208: WORK( N+I ) = SN
209: END IF
210: 10 CONTINUE
211: CALL DLARTG( D( N ), E( N ), CS, SN, R )
212: D( N ) = R
213: E( N ) = ZERO
214: IF( ROTATE ) THEN
215: WORK( N ) = CS
216: WORK( N+N ) = SN
217: END IF
218: IUPLO = 2
219: SQRE1 = 0
220: *
221: * Update singular vectors if desired.
222: *
223: IF( NCVT.GT.0 )
224: $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
225: $ WORK( NP1 ), VT, LDVT )
226: END IF
227: *
228: * If matrix lower bidiagonal, rotate to be upper bidiagonal
229: * by applying Givens rotations on the left.
230: *
231: IF( IUPLO.EQ.2 ) THEN
232: DO 20 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: IF( ROTATE ) THEN
238: WORK( I ) = CS
239: WORK( N+I ) = SN
240: END IF
241: 20 CONTINUE
242: *
243: * If matrix (N+1)-by-N lower bidiagonal, one additional
244: * rotation is needed.
245: *
246: IF( SQRE1.EQ.1 ) THEN
247: CALL DLARTG( D( N ), E( N ), CS, SN, R )
248: D( N ) = R
249: IF( ROTATE ) THEN
250: WORK( N ) = CS
251: WORK( N+N ) = SN
252: END IF
253: END IF
254: *
255: * Update singular vectors if desired.
256: *
257: IF( NRU.GT.0 ) THEN
258: IF( SQRE1.EQ.0 ) THEN
259: CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
260: $ WORK( NP1 ), U, LDU )
261: ELSE
262: CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
263: $ WORK( NP1 ), U, LDU )
264: END IF
265: END IF
266: IF( NCC.GT.0 ) THEN
267: IF( SQRE1.EQ.0 ) THEN
268: CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
269: $ WORK( NP1 ), C, LDC )
270: ELSE
271: CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
272: $ WORK( NP1 ), C, LDC )
273: END IF
274: END IF
275: END IF
276: *
277: * Call DBDSQR to compute the SVD of the reduced real
278: * N-by-N upper bidiagonal matrix.
279: *
280: CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
281: $ LDC, WORK, INFO )
282: *
283: * Sort the singular values into ascending order (insertion sort on
284: * singular values, but only one transposition per singular vector)
285: *
286: DO 40 I = 1, N
287: *
288: * Scan for smallest D(I).
289: *
290: ISUB = I
291: SMIN = D( I )
292: DO 30 J = I + 1, N
293: IF( D( J ).LT.SMIN ) THEN
294: ISUB = J
295: SMIN = D( J )
296: END IF
297: 30 CONTINUE
298: IF( ISUB.NE.I ) THEN
299: *
300: * Swap singular values and vectors.
301: *
302: D( ISUB ) = D( I )
303: D( I ) = SMIN
304: IF( NCVT.GT.0 )
305: $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
306: IF( NRU.GT.0 )
307: $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
308: IF( NCC.GT.0 )
309: $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
310: END IF
311: 40 CONTINUE
312: *
313: RETURN
314: *
315: * End of DLASDQ
316: *
317: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack