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