Mercurial > octave-nkf

comparison libcruft/lapack/dlalsd.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 DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
2 $ RANK, WORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
11 DOUBLE PRECISION RCOND
12 * ..
13 * .. Array Arguments ..
14 INTEGER IWORK( * )
15 DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLALSD uses the singular value decomposition of A to solve the least
22 * squares problem of finding X to minimize the Euclidean norm of each
23 * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
24 * are N-by-NRHS. The solution X overwrites B.
25 *
26 * The singular values of A smaller than RCOND times the largest
27 * singular value are treated as zero in solving the least squares
28 * problem; in this case a minimum norm solution is returned.
29 * The actual singular values are returned in D in ascending order.
30 *
31 * This code makes very mild assumptions about floating point
32 * arithmetic. It will work on machines with a guard digit in
33 * add/subtract, or on those binary machines without guard digits
34 * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
35 * It could conceivably fail on hexadecimal or decimal machines
36 * without guard digits, but we know of none.
37 *
38 * Arguments
39 * =========
40 *
41 * UPLO (input) CHARACTER*1
42 * = 'U': D and E define an upper bidiagonal matrix.
43 * = 'L': D and E define a lower bidiagonal matrix.
44 *
45 * SMLSIZ (input) INTEGER
46 * The maximum size of the subproblems at the bottom of the
47 * computation tree.
48 *
49 * N (input) INTEGER
50 * The dimension of the bidiagonal matrix. N >= 0.
51 *
52 * NRHS (input) INTEGER
53 * The number of columns of B. NRHS must be at least 1.
54 *
55 * D (input/output) DOUBLE PRECISION array, dimension (N)
56 * On entry D contains the main diagonal of the bidiagonal
57 * matrix. On exit, if INFO = 0, D contains its singular values.
58 *
59 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
60 * Contains the super-diagonal entries of the bidiagonal matrix.
61 * On exit, E has been destroyed.
62 *
63 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
64 * On input, B contains the right hand sides of the least
65 * squares problem. On output, B contains the solution X.
66 *
67 * LDB (input) INTEGER
68 * The leading dimension of B in the calling subprogram.
69 * LDB must be at least max(1,N).
70 *
71 * RCOND (input) DOUBLE PRECISION
72 * The singular values of A less than or equal to RCOND times
73 * the largest singular value are treated as zero in solving
74 * the least squares problem. If RCOND is negative,
75 * machine precision is used instead.
76 * For example, if diag(S)*X=B were the least squares problem,
77 * where diag(S) is a diagonal matrix of singular values, the
78 * solution would be X(i) = B(i) / S(i) if S(i) is greater than
79 * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
80 * RCOND*max(S).
81 *
82 * RANK (output) INTEGER
83 * The number of singular values of A greater than RCOND times
84 * the largest singular value.
85 *
86 * WORK (workspace) DOUBLE PRECISION array, dimension at least
87 * (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
88 * where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
89 *
90 * IWORK (workspace) INTEGER array, dimension at least
91 * (3*N*NLVL + 11*N)
92 *
93 * INFO (output) INTEGER
94 * = 0: successful exit.
95 * < 0: if INFO = -i, the i-th argument had an illegal value.
96 * > 0: The algorithm failed to compute an singular value while
97 * working on the submatrix lying in rows and columns
98 * INFO/(N+1) through MOD(INFO,N+1).
99 *
100 * Further Details
101 * ===============
102 *
103 * Based on contributions by
104 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
105 * California at Berkeley, USA
106 * Osni Marques, LBNL/NERSC, USA
107 *
108 * =====================================================================
109 *
110 * .. Parameters ..
111 DOUBLE PRECISION ZERO, ONE, TWO
112 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
113 * ..
114 * .. Local Scalars ..
115 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
116 $ GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
117 $ NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
118 $ SMLSZP, SQRE, ST, ST1, U, VT, Z
119 DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
120 * ..
121 * .. External Functions ..
122 INTEGER IDAMAX
123 DOUBLE PRECISION DLAMCH, DLANST
124 EXTERNAL IDAMAX, DLAMCH, DLANST
125 * ..
126 * .. External Subroutines ..
127 EXTERNAL DCOPY, DGEMM, DLACPY, DLALSA, DLARTG, DLASCL,
128 $ DLASDA, DLASDQ, DLASET, DLASRT, DROT, XERBLA
129 * ..
130 * .. Intrinsic Functions ..
131 INTRINSIC ABS, DBLE, INT, LOG, SIGN
132 * ..
133 * .. Executable Statements ..
134 *
135 * Test the input parameters.
136 *
137 INFO = 0
138 *
139 IF( N.LT.0 ) THEN
140 INFO = -3
141 ELSE IF( NRHS.LT.1 ) THEN
142 INFO = -4
143 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
144 INFO = -8
145 END IF
146 IF( INFO.NE.0 ) THEN
147 CALL XERBLA( 'DLALSD', -INFO )
148 RETURN
149 END IF
150 *
151 EPS = DLAMCH( 'Epsilon' )
152 *
153 * Set up the tolerance.
154 *
155 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
156 RCND = EPS
157 ELSE
158 RCND = RCOND
159 END IF
160 *
161 RANK = 0
162 *
163 * Quick return if possible.
164 *
165 IF( N.EQ.0 ) THEN
166 RETURN
167 ELSE IF( N.EQ.1 ) THEN
168 IF( D( 1 ).EQ.ZERO ) THEN
169 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
170 ELSE
171 RANK = 1
172 CALL DLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
173 D( 1 ) = ABS( D( 1 ) )
174 END IF
175 RETURN
176 END IF
177 *
178 * Rotate the matrix if it is lower bidiagonal.
179 *
180 IF( UPLO.EQ.'L' ) THEN
181 DO 10 I = 1, N - 1
182 CALL DLARTG( D( I ), E( I ), CS, SN, R )
183 D( I ) = R
184 E( I ) = SN*D( I+1 )
185 D( I+1 ) = CS*D( I+1 )
186 IF( NRHS.EQ.1 ) THEN
187 CALL DROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
188 ELSE
189 WORK( I*2-1 ) = CS
190 WORK( I*2 ) = SN
191 END IF
192 10 CONTINUE
193 IF( NRHS.GT.1 ) THEN
194 DO 30 I = 1, NRHS
195 DO 20 J = 1, N - 1
196 CS = WORK( J*2-1 )
197 SN = WORK( J*2 )
198 CALL DROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
199 20 CONTINUE
200 30 CONTINUE
201 END IF
202 END IF
203 *
204 * Scale.
205 *
206 NM1 = N - 1
207 ORGNRM = DLANST( 'M', N, D, E )
208 IF( ORGNRM.EQ.ZERO ) THEN
209 CALL DLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
210 RETURN
211 END IF
212 *
213 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
214 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
215 *
216 * If N is smaller than the minimum divide size SMLSIZ, then solve
217 * the problem with another solver.
218 *
219 IF( N.LE.SMLSIZ ) THEN
220 NWORK = 1 + N*N
221 CALL DLASET( 'A', N, N, ZERO, ONE, WORK, N )
222 CALL DLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
223 $ LDB, WORK( NWORK ), INFO )
224 IF( INFO.NE.0 ) THEN
225 RETURN
226 END IF
227 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
228 DO 40 I = 1, N
229 IF( D( I ).LE.TOL ) THEN
230 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
231 ELSE
232 CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
233 $ LDB, INFO )
234 RANK = RANK + 1
235 END IF
236 40 CONTINUE
237 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
238 $ WORK( NWORK ), N )
239 CALL DLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
240 *
241 * Unscale.
242 *
243 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
244 CALL DLASRT( 'D', N, D, INFO )
245 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
246 *
247 RETURN
248 END IF
249 *
250 * Book-keeping and setting up some constants.
251 *
252 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
253 *
254 SMLSZP = SMLSIZ + 1
255 *
256 U = 1
257 VT = 1 + SMLSIZ*N
258 DIFL = VT + SMLSZP*N
259 DIFR = DIFL + NLVL*N
260 Z = DIFR + NLVL*N*2
261 C = Z + NLVL*N
262 S = C + N
263 POLES = S + N
264 GIVNUM = POLES + 2*NLVL*N
265 BX = GIVNUM + 2*NLVL*N
266 NWORK = BX + N*NRHS
267 *
268 SIZEI = 1 + N
269 K = SIZEI + N
270 GIVPTR = K + N
271 PERM = GIVPTR + N
272 GIVCOL = PERM + NLVL*N
273 IWK = GIVCOL + NLVL*N*2
274 *
275 ST = 1
276 SQRE = 0
277 ICMPQ1 = 1
278 ICMPQ2 = 0
279 NSUB = 0
280 *
281 DO 50 I = 1, N
282 IF( ABS( D( I ) ).LT.EPS ) THEN
283 D( I ) = SIGN( EPS, D( I ) )
284 END IF
285 50 CONTINUE
286 *
287 DO 60 I = 1, NM1
288 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
289 NSUB = NSUB + 1
290 IWORK( NSUB ) = ST
291 *
292 * Subproblem found. First determine its size and then
293 * apply divide and conquer on it.
294 *
295 IF( I.LT.NM1 ) THEN
296 *
297 * A subproblem with E(I) small for I < NM1.
298 *
299 NSIZE = I - ST + 1
300 IWORK( SIZEI+NSUB-1 ) = NSIZE
301 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
302 *
303 * A subproblem with E(NM1) not too small but I = NM1.
304 *
305 NSIZE = N - ST + 1
306 IWORK( SIZEI+NSUB-1 ) = NSIZE
307 ELSE
308 *
309 * A subproblem with E(NM1) small. This implies an
310 * 1-by-1 subproblem at D(N), which is not solved
311 * explicitly.
312 *
313 NSIZE = I - ST + 1
314 IWORK( SIZEI+NSUB-1 ) = NSIZE
315 NSUB = NSUB + 1
316 IWORK( NSUB ) = N
317 IWORK( SIZEI+NSUB-1 ) = 1
318 CALL DCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
319 END IF
320 ST1 = ST - 1
321 IF( NSIZE.EQ.1 ) THEN
322 *
323 * This is a 1-by-1 subproblem and is not solved
324 * explicitly.
325 *
326 CALL DCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
327 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
328 *
329 * This is a small subproblem and is solved by DLASDQ.
330 *
331 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
332 $ WORK( VT+ST1 ), N )
333 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
334 $ E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
335 $ N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
336 IF( INFO.NE.0 ) THEN
337 RETURN
338 END IF
339 CALL DLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
340 $ WORK( BX+ST1 ), N )
341 ELSE
342 *
343 * A large problem. Solve it using divide and conquer.
344 *
345 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
346 $ E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
347 $ IWORK( K+ST1 ), WORK( DIFL+ST1 ),
348 $ WORK( DIFR+ST1 ), WORK( Z+ST1 ),
349 $ WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
350 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
351 $ WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
352 $ WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
353 $ INFO )
354 IF( INFO.NE.0 ) THEN
355 RETURN
356 END IF
357 BXST = BX + ST1
358 CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
359 $ LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
360 $ WORK( VT+ST1 ), IWORK( K+ST1 ),
361 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
362 $ WORK( Z+ST1 ), WORK( POLES+ST1 ),
363 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
364 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
365 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
366 $ IWORK( IWK ), INFO )
367 IF( INFO.NE.0 ) THEN
368 RETURN
369 END IF
370 END IF
371 ST = I + 1
372 END IF
373 60 CONTINUE
374 *
375 * Apply the singular values and treat the tiny ones as zero.
376 *
377 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
378 *
379 DO 70 I = 1, N
380 *
381 * Some of the elements in D can be negative because 1-by-1
382 * subproblems were not solved explicitly.
383 *
384 IF( ABS( D( I ) ).LE.TOL ) THEN
385 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
386 ELSE
387 RANK = RANK + 1
388 CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
389 $ WORK( BX+I-1 ), N, INFO )
390 END IF
391 D( I ) = ABS( D( I ) )
392 70 CONTINUE
393 *
394 * Now apply back the right singular vectors.
395 *
396 ICMPQ2 = 1
397 DO 80 I = 1, NSUB
398 ST = IWORK( I )
399 ST1 = ST - 1
400 NSIZE = IWORK( SIZEI+I-1 )
401 BXST = BX + ST1
402 IF( NSIZE.EQ.1 ) THEN
403 CALL DCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
404 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
405 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
406 $ WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
407 $ B( ST, 1 ), LDB )
408 ELSE
409 CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
410 $ B( ST, 1 ), LDB, WORK( U+ST1 ), N,
411 $ WORK( VT+ST1 ), IWORK( K+ST1 ),
412 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
413 $ WORK( Z+ST1 ), WORK( POLES+ST1 ),
414 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
415 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
416 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
417 $ IWORK( IWK ), INFO )
418 IF( INFO.NE.0 ) THEN
419 RETURN
420 END IF
421 END IF
422 80 CONTINUE
423 *
424 * Unscale and sort the singular values.
425 *
426 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
427 CALL DLASRT( 'D', N, D, INFO )
428 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
429 *
430 RETURN
431 *
432 * End of DLALSD
433 *
434 END