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