Mercurial > octave-nkf

comparison libcruft/lapack/dlasda.f @ 7072:b48d486f641d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-26 15:52:57 by jwe]
author jwe
date Fri, 26 Oct 2007 15:52:58 +0000
parents
children
comparison
equal deleted inserted replaced
7071:c3b479e753dd 7072:b48d486f641d
1 SUBROUTINE DLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
2 $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
3 $ PERM, GIVNUM, C, S, WORK, IWORK, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.1) --
6 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
11 * ..
12 * .. Array Arguments ..
13 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
14 $ K( * ), PERM( LDGCOL, * )
15 DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
16 $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
17 $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
18 $ Z( LDU, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * Using a divide and conquer approach, DLASDA computes the singular
25 * value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
26 * B with diagonal D and offdiagonal E, where M = N + SQRE. The
27 * algorithm computes the singular values in the SVD B = U * S * VT.
28 * The orthogonal matrices U and VT are optionally computed in
29 * compact form.
30 *
31 * A related subroutine, DLASD0, computes the singular values and
32 * the singular vectors in explicit form.
33 *
34 * Arguments
35 * =========
36 *
37 * ICOMPQ (input) INTEGER
38 * Specifies whether singular vectors are to be computed
39 * in compact form, as follows
40 * = 0: Compute singular values only.
41 * = 1: Compute singular vectors of upper bidiagonal
42 * matrix in compact form.
43 *
44 * SMLSIZ (input) INTEGER
45 * The maximum size of the subproblems at the bottom of the
46 * computation tree.
47 *
48 * N (input) INTEGER
49 * The row dimension of the upper bidiagonal matrix. This is
50 * also the dimension of the main diagonal array D.
51 *
52 * SQRE (input) INTEGER
53 * Specifies the column dimension of the bidiagonal matrix.
54 * = 0: The bidiagonal matrix has column dimension M = N;
55 * = 1: The bidiagonal matrix has column dimension M = N + 1.
56 *
57 * D (input/output) DOUBLE PRECISION array, dimension ( N )
58 * On entry D contains the main diagonal of the bidiagonal
59 * matrix. On exit D, if INFO = 0, contains its singular values.
60 *
61 * E (input) DOUBLE PRECISION array, dimension ( M-1 )
62 * Contains the subdiagonal entries of the bidiagonal matrix.
63 * On exit, E has been destroyed.
64 *
65 * U (output) DOUBLE PRECISION array,
66 * dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
67 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
68 * singular vector matrices of all subproblems at the bottom
69 * level.
70 *
71 * LDU (input) INTEGER, LDU = > N.
72 * The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
73 * GIVNUM, and Z.
74 *
75 * VT (output) DOUBLE PRECISION array,
76 * dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
77 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT' contains the right
78 * singular vector matrices of all subproblems at the bottom
79 * level.
80 *
81 * K (output) INTEGER array,
82 * dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
83 * If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
84 * secular equation on the computation tree.
85 *
86 * DIFL (output) DOUBLE PRECISION array, dimension ( LDU, NLVL ),
87 * where NLVL = floor(log_2 (N/SMLSIZ))).
88 *
89 * DIFR (output) DOUBLE PRECISION array,
90 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
91 * dimension ( N ) if ICOMPQ = 0.
92 * If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
93 * record distances between singular values on the I-th
94 * level and singular values on the (I -1)-th level, and
95 * DIFR(1:N, 2 * I ) contains the normalizing factors for
96 * the right singular vector matrix. See DLASD8 for details.
97 *
98 * Z (output) DOUBLE PRECISION array,
99 * dimension ( LDU, NLVL ) if ICOMPQ = 1 and
100 * dimension ( N ) if ICOMPQ = 0.
101 * The first K elements of Z(1, I) contain the components of
102 * the deflation-adjusted updating row vector for subproblems
103 * on the I-th level.
104 *
105 * POLES (output) DOUBLE PRECISION array,
106 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
107 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
108 * POLES(1, 2*I) contain the new and old singular values
109 * involved in the secular equations on the I-th level.
110 *
111 * GIVPTR (output) INTEGER array,
112 * dimension ( N ) if ICOMPQ = 1, and not referenced if
113 * ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
114 * the number of Givens rotations performed on the I-th
115 * problem on the computation tree.
116 *
117 * GIVCOL (output) INTEGER array,
118 * dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
119 * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
120 * GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
121 * of Givens rotations performed on the I-th level on the
122 * computation tree.
123 *
124 * LDGCOL (input) INTEGER, LDGCOL = > N.
125 * The leading dimension of arrays GIVCOL and PERM.
126 *
127 * PERM (output) INTEGER array,
128 * dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced
129 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
130 * permutations done on the I-th level of the computation tree.
131 *
132 * GIVNUM (output) DOUBLE PRECISION array,
133 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not
134 * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
135 * GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
136 * values of Givens rotations performed on the I-th level on
137 * the computation tree.
138 *
139 * C (output) DOUBLE PRECISION array,
140 * dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
141 * If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
142 * C( I ) contains the C-value of a Givens rotation related to
143 * the right null space of the I-th subproblem.
144 *
145 * S (output) DOUBLE PRECISION array, dimension ( N ) if
146 * ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
147 * and the I-th subproblem is not square, on exit, S( I )
148 * contains the S-value of a Givens rotation related to
149 * the right null space of the I-th subproblem.
150 *
151 * WORK (workspace) DOUBLE PRECISION array, dimension
152 * (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
153 *
154 * IWORK (workspace) INTEGER array.
155 * Dimension must be at least (7 * N).
156 *
157 * INFO (output) INTEGER
158 * = 0: successful exit.
159 * < 0: if INFO = -i, the i-th argument had an illegal value.
160 * > 0: if INFO = 1, an singular value did not converge
161 *
162 * Further Details
163 * ===============
164 *
165 * Based on contributions by
166 * Ming Gu and Huan Ren, Computer Science Division, University of
167 * California at Berkeley, USA
168 *
169 * =====================================================================
170 *
171 * .. Parameters ..
172 DOUBLE PRECISION ZERO, ONE
173 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
174 * ..
175 * .. Local Scalars ..
176 INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
177 $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
178 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU,
179 $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI
180 DOUBLE PRECISION ALPHA, BETA
181 * ..
182 * .. External Subroutines ..
183 EXTERNAL DCOPY, DLASD6, DLASDQ, DLASDT, DLASET, XERBLA
184 * ..
185 * .. Executable Statements ..
186 *
187 * Test the input parameters.
188 *
189 INFO = 0
190 *
191 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
192 INFO = -1
193 ELSE IF( SMLSIZ.LT.3 ) THEN
194 INFO = -2
195 ELSE IF( N.LT.0 ) THEN
196 INFO = -3
197 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
198 INFO = -4
199 ELSE IF( LDU.LT.( N+SQRE ) ) THEN
200 INFO = -8
201 ELSE IF( LDGCOL.LT.N ) THEN
202 INFO = -17
203 END IF
204 IF( INFO.NE.0 ) THEN
205 CALL XERBLA( 'DLASDA', -INFO )
206 RETURN
207 END IF
208 *
209 M = N + SQRE
210 *
211 * If the input matrix is too small, call DLASDQ to find the SVD.
212 *
213 IF( N.LE.SMLSIZ ) THEN
214 IF( ICOMPQ.EQ.0 ) THEN
215 CALL DLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU,
216 $ U, LDU, WORK, INFO )
217 ELSE
218 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU,
219 $ U, LDU, WORK, INFO )
220 END IF
221 RETURN
222 END IF
223 *
224 * Book-keeping and set up the computation tree.
225 *
226 INODE = 1
227 NDIML = INODE + N
228 NDIMR = NDIML + N
229 IDXQ = NDIMR + N
230 IWK = IDXQ + N
231 *
232 NCC = 0
233 NRU = 0
234 *
235 SMLSZP = SMLSIZ + 1
236 VF = 1
237 VL = VF + M
238 NWORK1 = VL + M
239 NWORK2 = NWORK1 + SMLSZP*SMLSZP
240 *
241 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
242 $ IWORK( NDIMR ), SMLSIZ )
243 *
244 * for the nodes on bottom level of the tree, solve
245 * their subproblems by DLASDQ.
246 *
247 NDB1 = ( ND+1 ) / 2
248 DO 30 I = NDB1, ND
249 *
250 * IC : center row of each node
251 * NL : number of rows of left subproblem
252 * NR : number of rows of right subproblem
253 * NLF: starting row of the left subproblem
254 * NRF: starting row of the right subproblem
255 *
256 I1 = I - 1
257 IC = IWORK( INODE+I1 )
258 NL = IWORK( NDIML+I1 )
259 NLP1 = NL + 1
260 NR = IWORK( NDIMR+I1 )
261 NLF = IC - NL
262 NRF = IC + 1
263 IDXQI = IDXQ + NLF - 2
264 VFI = VF + NLF - 1
265 VLI = VL + NLF - 1
266 SQREI = 1
267 IF( ICOMPQ.EQ.0 ) THEN
268 CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ),
269 $ SMLSZP )
270 CALL DLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ),
271 $ E( NLF ), WORK( NWORK1 ), SMLSZP,
272 $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL,
273 $ WORK( NWORK2 ), INFO )
274 ITEMP = NWORK1 + NL*SMLSZP
275 CALL DCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
276 CALL DCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
277 ELSE
278 CALL DLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU )
279 CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU )
280 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ),
281 $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU,
282 $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO )
283 CALL DCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 )
284 CALL DCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 )
285 END IF
286 IF( INFO.NE.0 ) THEN
287 RETURN
288 END IF
289 DO 10 J = 1, NL
290 IWORK( IDXQI+J ) = J
291 10 CONTINUE
292 IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN
293 SQREI = 0
294 ELSE
295 SQREI = 1
296 END IF
297 IDXQI = IDXQI + NLP1
298 VFI = VFI + NLP1
299 VLI = VLI + NLP1
300 NRP1 = NR + SQREI
301 IF( ICOMPQ.EQ.0 ) THEN
302 CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ),
303 $ SMLSZP )
304 CALL DLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ),
305 $ E( NRF ), WORK( NWORK1 ), SMLSZP,
306 $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR,
307 $ WORK( NWORK2 ), INFO )
308 ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP
309 CALL DCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
310 CALL DCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
311 ELSE
312 CALL DLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU )
313 CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU )
314 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ),
315 $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU,
316 $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO )
317 CALL DCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 )
318 CALL DCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 )
319 END IF
320 IF( INFO.NE.0 ) THEN
321 RETURN
322 END IF
323 DO 20 J = 1, NR
324 IWORK( IDXQI+J ) = J
325 20 CONTINUE
326 30 CONTINUE
327 *
328 * Now conquer each subproblem bottom-up.
329 *
330 J = 2**NLVL
331 DO 50 LVL = NLVL, 1, -1
332 LVL2 = LVL*2 - 1
333 *
334 * Find the first node LF and last node LL on
335 * the current level LVL.
336 *
337 IF( LVL.EQ.1 ) THEN
338 LF = 1
339 LL = 1
340 ELSE
341 LF = 2**( LVL-1 )
342 LL = 2*LF - 1
343 END IF
344 DO 40 I = LF, LL
345 IM1 = I - 1
346 IC = IWORK( INODE+IM1 )
347 NL = IWORK( NDIML+IM1 )
348 NR = IWORK( NDIMR+IM1 )
349 NLF = IC - NL
350 NRF = IC + 1
351 IF( I.EQ.LL ) THEN
352 SQREI = SQRE
353 ELSE
354 SQREI = 1
355 END IF
356 VFI = VF + NLF - 1
357 VLI = VL + NLF - 1
358 IDXQI = IDXQ + NLF - 1
359 ALPHA = D( IC )
360 BETA = E( IC )
361 IF( ICOMPQ.EQ.0 ) THEN
362 CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
363 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
364 $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL,
365 $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z,
366 $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ),
367 $ IWORK( IWK ), INFO )
368 ELSE
369 J = J - 1
370 CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
371 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
372 $ IWORK( IDXQI ), PERM( NLF, LVL ),
373 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
374 $ GIVNUM( NLF, LVL2 ), LDU,
375 $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ),
376 $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ),
377 $ C( J ), S( J ), WORK( NWORK1 ),
378 $ IWORK( IWK ), INFO )
379 END IF
380 IF( INFO.NE.0 ) THEN
381 RETURN
382 END IF
383 40 CONTINUE
384 50 CONTINUE
385 *
386 RETURN
387 *
388 * End of DLASDA
389 *
390 END