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: *> \date September 2012
168: *
169: *> \ingroup doubleOTHERcomputational
170: *
171: *> \par Contributors:
172: * ==================
173: *>
174: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
175: *> California at Berkeley, USA \n
176: *> Osni Marques, LBNL/NERSC, USA \n
177: *
178: * =====================================================================
179: SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
180: $ RANK, WORK, IWORK, INFO )
181: *
182: * -- LAPACK computational routine (version 3.4.2) --
183: * -- LAPACK is a software package provided by Univ. of Tennessee, --
184: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
185: * September 2012
186: *
187: * .. Scalar Arguments ..
188: CHARACTER UPLO
189: INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
190: DOUBLE PRECISION RCOND
191: * ..
192: * .. Array Arguments ..
193: INTEGER IWORK( * )
194: DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * )
195: * ..
196: *
197: * =====================================================================
198: *
199: * .. Parameters ..
200: DOUBLE PRECISION ZERO, ONE, TWO
201: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
202: * ..
203: * .. Local Scalars ..
204: INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
205: $ GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
206: $ NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
207: $ SMLSZP, SQRE, ST, ST1, U, VT, Z
208: DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL
209: * ..
210: * .. External Functions ..
211: INTEGER IDAMAX
212: DOUBLE PRECISION DLAMCH, DLANST
213: EXTERNAL IDAMAX, DLAMCH, DLANST
214: * ..
215: * .. External Subroutines ..
216: EXTERNAL DCOPY, DGEMM, DLACPY, DLALSA, DLARTG, DLASCL,
217: $ DLASDA, DLASDQ, DLASET, DLASRT, DROT, XERBLA
218: * ..
219: * .. Intrinsic Functions ..
220: INTRINSIC ABS, DBLE, INT, LOG, SIGN
221: * ..
222: * .. Executable Statements ..
223: *
224: * Test the input parameters.
225: *
226: INFO = 0
227: *
228: IF( N.LT.0 ) THEN
229: INFO = -3
230: ELSE IF( NRHS.LT.1 ) THEN
231: INFO = -4
232: ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
233: INFO = -8
234: END IF
235: IF( INFO.NE.0 ) THEN
236: CALL XERBLA( 'DLALSD', -INFO )
237: RETURN
238: END IF
239: *
240: EPS = DLAMCH( 'Epsilon' )
241: *
242: * Set up the tolerance.
243: *
244: IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
245: RCND = EPS
246: ELSE
247: RCND = RCOND
248: END IF
249: *
250: RANK = 0
251: *
252: * Quick return if possible.
253: *
254: IF( N.EQ.0 ) THEN
255: RETURN
256: ELSE IF( N.EQ.1 ) THEN
257: IF( D( 1 ).EQ.ZERO ) THEN
258: CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
259: ELSE
260: RANK = 1
261: CALL DLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
262: D( 1 ) = ABS( D( 1 ) )
263: END IF
264: RETURN
265: END IF
266: *
267: * Rotate the matrix if it is lower bidiagonal.
268: *
269: IF( UPLO.EQ.'L' ) THEN
270: DO 10 I = 1, N - 1
271: CALL DLARTG( D( I ), E( I ), CS, SN, R )
272: D( I ) = R
273: E( I ) = SN*D( I+1 )
274: D( I+1 ) = CS*D( I+1 )
275: IF( NRHS.EQ.1 ) THEN
276: CALL DROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
277: ELSE
278: WORK( I*2-1 ) = CS
279: WORK( I*2 ) = SN
280: END IF
281: 10 CONTINUE
282: IF( NRHS.GT.1 ) THEN
283: DO 30 I = 1, NRHS
284: DO 20 J = 1, N - 1
285: CS = WORK( J*2-1 )
286: SN = WORK( J*2 )
287: CALL DROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
288: 20 CONTINUE
289: 30 CONTINUE
290: END IF
291: END IF
292: *
293: * Scale.
294: *
295: NM1 = N - 1
296: ORGNRM = DLANST( 'M', N, D, E )
297: IF( ORGNRM.EQ.ZERO ) THEN
298: CALL DLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
299: RETURN
300: END IF
301: *
302: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
303: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
304: *
305: * If N is smaller than the minimum divide size SMLSIZ, then solve
306: * the problem with another solver.
307: *
308: IF( N.LE.SMLSIZ ) THEN
309: NWORK = 1 + N*N
310: CALL DLASET( 'A', N, N, ZERO, ONE, WORK, N )
311: CALL DLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
312: $ LDB, WORK( NWORK ), INFO )
313: IF( INFO.NE.0 ) THEN
314: RETURN
315: END IF
316: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
317: DO 40 I = 1, N
318: IF( D( I ).LE.TOL ) THEN
319: CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
320: ELSE
321: CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
322: $ LDB, INFO )
323: RANK = RANK + 1
324: END IF
325: 40 CONTINUE
326: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
327: $ WORK( NWORK ), N )
328: CALL DLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
329: *
330: * Unscale.
331: *
332: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
333: CALL DLASRT( 'D', N, D, INFO )
334: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
335: *
336: RETURN
337: END IF
338: *
339: * Book-keeping and setting up some constants.
340: *
341: NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
342: *
343: SMLSZP = SMLSIZ + 1
344: *
345: U = 1
346: VT = 1 + SMLSIZ*N
347: DIFL = VT + SMLSZP*N
348: DIFR = DIFL + NLVL*N
349: Z = DIFR + NLVL*N*2
350: C = Z + NLVL*N
351: S = C + N
352: POLES = S + N
353: GIVNUM = POLES + 2*NLVL*N
354: BX = GIVNUM + 2*NLVL*N
355: NWORK = BX + N*NRHS
356: *
357: SIZEI = 1 + N
358: K = SIZEI + N
359: GIVPTR = K + N
360: PERM = GIVPTR + N
361: GIVCOL = PERM + NLVL*N
362: IWK = GIVCOL + NLVL*N*2
363: *
364: ST = 1
365: SQRE = 0
366: ICMPQ1 = 1
367: ICMPQ2 = 0
368: NSUB = 0
369: *
370: DO 50 I = 1, N
371: IF( ABS( D( I ) ).LT.EPS ) THEN
372: D( I ) = SIGN( EPS, D( I ) )
373: END IF
374: 50 CONTINUE
375: *
376: DO 60 I = 1, NM1
377: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
378: NSUB = NSUB + 1
379: IWORK( NSUB ) = ST
380: *
381: * Subproblem found. First determine its size and then
382: * apply divide and conquer on it.
383: *
384: IF( I.LT.NM1 ) THEN
385: *
386: * A subproblem with E(I) small for I < NM1.
387: *
388: NSIZE = I - ST + 1
389: IWORK( SIZEI+NSUB-1 ) = NSIZE
390: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
391: *
392: * A subproblem with E(NM1) not too small but I = NM1.
393: *
394: NSIZE = N - ST + 1
395: IWORK( SIZEI+NSUB-1 ) = NSIZE
396: ELSE
397: *
398: * A subproblem with E(NM1) small. This implies an
399: * 1-by-1 subproblem at D(N), which is not solved
400: * explicitly.
401: *
402: NSIZE = I - ST + 1
403: IWORK( SIZEI+NSUB-1 ) = NSIZE
404: NSUB = NSUB + 1
405: IWORK( NSUB ) = N
406: IWORK( SIZEI+NSUB-1 ) = 1
407: CALL DCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
408: END IF
409: ST1 = ST - 1
410: IF( NSIZE.EQ.1 ) THEN
411: *
412: * This is a 1-by-1 subproblem and is not solved
413: * explicitly.
414: *
415: CALL DCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
416: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
417: *
418: * This is a small subproblem and is solved by DLASDQ.
419: *
420: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
421: $ WORK( VT+ST1 ), N )
422: CALL DLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
423: $ E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
424: $ N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
425: IF( INFO.NE.0 ) THEN
426: RETURN
427: END IF
428: CALL DLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
429: $ WORK( BX+ST1 ), N )
430: ELSE
431: *
432: * A large problem. Solve it using divide and conquer.
433: *
434: CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
435: $ E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
436: $ IWORK( K+ST1 ), WORK( DIFL+ST1 ),
437: $ WORK( DIFR+ST1 ), WORK( Z+ST1 ),
438: $ WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
439: $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
440: $ WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
441: $ WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
442: $ INFO )
443: IF( INFO.NE.0 ) THEN
444: RETURN
445: END IF
446: BXST = BX + ST1
447: CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
448: $ LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
449: $ WORK( VT+ST1 ), IWORK( K+ST1 ),
450: $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
451: $ WORK( Z+ST1 ), WORK( POLES+ST1 ),
452: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
453: $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
454: $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
455: $ IWORK( IWK ), INFO )
456: IF( INFO.NE.0 ) THEN
457: RETURN
458: END IF
459: END IF
460: ST = I + 1
461: END IF
462: 60 CONTINUE
463: *
464: * Apply the singular values and treat the tiny ones as zero.
465: *
466: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
467: *
468: DO 70 I = 1, N
469: *
470: * Some of the elements in D can be negative because 1-by-1
471: * subproblems were not solved explicitly.
472: *
473: IF( ABS( D( I ) ).LE.TOL ) THEN
474: CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
475: ELSE
476: RANK = RANK + 1
477: CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
478: $ WORK( BX+I-1 ), N, INFO )
479: END IF
480: D( I ) = ABS( D( I ) )
481: 70 CONTINUE
482: *
483: * Now apply back the right singular vectors.
484: *
485: ICMPQ2 = 1
486: DO 80 I = 1, NSUB
487: ST = IWORK( I )
488: ST1 = ST - 1
489: NSIZE = IWORK( SIZEI+I-1 )
490: BXST = BX + ST1
491: IF( NSIZE.EQ.1 ) THEN
492: CALL DCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
493: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
494: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
495: $ WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
496: $ B( ST, 1 ), LDB )
497: ELSE
498: CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
499: $ B( ST, 1 ), LDB, WORK( U+ST1 ), N,
500: $ WORK( VT+ST1 ), IWORK( K+ST1 ),
501: $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
502: $ WORK( Z+ST1 ), WORK( POLES+ST1 ),
503: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
504: $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
505: $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
506: $ IWORK( IWK ), INFO )
507: IF( INFO.NE.0 ) THEN
508: RETURN
509: END IF
510: END IF
511: 80 CONTINUE
512: *
513: * Unscale and sort the singular values.
514: *
515: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
516: CALL DLASRT( 'D', N, D, INFO )
517: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
518: *
519: RETURN
520: *
521: * End of DLALSD
522: *
523: 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