Annotation of rpl/lapack/lapack/dlalsd.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLALSD
! 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 November 2011
! 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: * =====================================================================
1.1 bertrand 179: SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
180: $ RANK, WORK, IWORK, INFO )
181: *
1.9 ! bertrand 182: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 183: * -- LAPACK is a software package provided by Univ. of Tennessee, --
184: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 185: * November 2011
1.1 bertrand 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