Annotation of rpl/lapack/lapack/dlasdq.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLASDQ
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASDQ + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasdq.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasdq.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasdq.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
! 22: * U, LDU, C, LDC, WORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER UPLO
! 26: * INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
! 30: * $ VT( LDVT, * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLASDQ computes the singular value decomposition (SVD) of a real
! 40: *> (upper or lower) bidiagonal matrix with diagonal D and offdiagonal
! 41: *> E, accumulating the transformations if desired. Letting B denote
! 42: *> the input bidiagonal matrix, the algorithm computes orthogonal
! 43: *> matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose
! 44: *> of P). The singular values S are overwritten on D.
! 45: *>
! 46: *> The input matrix U is changed to U * Q if desired.
! 47: *> The input matrix VT is changed to P**T * VT if desired.
! 48: *> The input matrix C is changed to Q**T * C if desired.
! 49: *>
! 50: *> See "Computing Small Singular Values of Bidiagonal Matrices With
! 51: *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
! 52: *> LAPACK Working Note #3, for a detailed description of the algorithm.
! 53: *> \endverbatim
! 54: *
! 55: * Arguments:
! 56: * ==========
! 57: *
! 58: *> \param[in] UPLO
! 59: *> \verbatim
! 60: *> UPLO is CHARACTER*1
! 61: *> On entry, UPLO specifies whether the input bidiagonal matrix
! 62: *> is upper or lower bidiagonal, and wether it is square are
! 63: *> not.
! 64: *> UPLO = 'U' or 'u' B is upper bidiagonal.
! 65: *> UPLO = 'L' or 'l' B is lower bidiagonal.
! 66: *> \endverbatim
! 67: *>
! 68: *> \param[in] SQRE
! 69: *> \verbatim
! 70: *> SQRE is INTEGER
! 71: *> = 0: then the input matrix is N-by-N.
! 72: *> = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and
! 73: *> (N+1)-by-N if UPLU = 'L'.
! 74: *>
! 75: *> The bidiagonal matrix has
! 76: *> N = NL + NR + 1 rows and
! 77: *> M = N + SQRE >= N columns.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] N
! 81: *> \verbatim
! 82: *> N is INTEGER
! 83: *> On entry, N specifies the number of rows and columns
! 84: *> in the matrix. N must be at least 0.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in] NCVT
! 88: *> \verbatim
! 89: *> NCVT is INTEGER
! 90: *> On entry, NCVT specifies the number of columns of
! 91: *> the matrix VT. NCVT must be at least 0.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] NRU
! 95: *> \verbatim
! 96: *> NRU is INTEGER
! 97: *> On entry, NRU specifies the number of rows of
! 98: *> the matrix U. NRU must be at least 0.
! 99: *> \endverbatim
! 100: *>
! 101: *> \param[in] NCC
! 102: *> \verbatim
! 103: *> NCC is INTEGER
! 104: *> On entry, NCC specifies the number of columns of
! 105: *> the matrix C. NCC must be at least 0.
! 106: *> \endverbatim
! 107: *>
! 108: *> \param[in,out] D
! 109: *> \verbatim
! 110: *> D is DOUBLE PRECISION array, dimension (N)
! 111: *> On entry, D contains the diagonal entries of the
! 112: *> bidiagonal matrix whose SVD is desired. On normal exit,
! 113: *> D contains the singular values in ascending order.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[in,out] E
! 117: *> \verbatim
! 118: *> E is DOUBLE PRECISION array.
! 119: *> dimension is (N-1) if SQRE = 0 and N if SQRE = 1.
! 120: *> On entry, the entries of E contain the offdiagonal entries
! 121: *> of the bidiagonal matrix whose SVD is desired. On normal
! 122: *> exit, E will contain 0. If the algorithm does not converge,
! 123: *> D and E will contain the diagonal and superdiagonal entries
! 124: *> of a bidiagonal matrix orthogonally equivalent to the one
! 125: *> given as input.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in,out] VT
! 129: *> \verbatim
! 130: *> VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
! 131: *> On entry, contains a matrix which on exit has been
! 132: *> premultiplied by P**T, dimension N-by-NCVT if SQRE = 0
! 133: *> and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0).
! 134: *> \endverbatim
! 135: *>
! 136: *> \param[in] LDVT
! 137: *> \verbatim
! 138: *> LDVT is INTEGER
! 139: *> On entry, LDVT specifies the leading dimension of VT as
! 140: *> declared in the calling (sub) program. LDVT must be at
! 141: *> least 1. If NCVT is nonzero LDVT must also be at least N.
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[in,out] U
! 145: *> \verbatim
! 146: *> U is DOUBLE PRECISION array, dimension (LDU, N)
! 147: *> On entry, contains a matrix which on exit has been
! 148: *> postmultiplied by Q, dimension NRU-by-N if SQRE = 0
! 149: *> and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0).
! 150: *> \endverbatim
! 151: *>
! 152: *> \param[in] LDU
! 153: *> \verbatim
! 154: *> LDU is INTEGER
! 155: *> On entry, LDU specifies the leading dimension of U as
! 156: *> declared in the calling (sub) program. LDU must be at
! 157: *> least max( 1, NRU ) .
! 158: *> \endverbatim
! 159: *>
! 160: *> \param[in,out] C
! 161: *> \verbatim
! 162: *> C is DOUBLE PRECISION array, dimension (LDC, NCC)
! 163: *> On entry, contains an N-by-NCC matrix which on exit
! 164: *> has been premultiplied by Q**T dimension N-by-NCC if SQRE = 0
! 165: *> and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0).
! 166: *> \endverbatim
! 167: *>
! 168: *> \param[in] LDC
! 169: *> \verbatim
! 170: *> LDC is INTEGER
! 171: *> On entry, LDC specifies the leading dimension of C as
! 172: *> declared in the calling (sub) program. LDC must be at
! 173: *> least 1. If NCC is nonzero, LDC must also be at least N.
! 174: *> \endverbatim
! 175: *>
! 176: *> \param[out] WORK
! 177: *> \verbatim
! 178: *> WORK is DOUBLE PRECISION array, dimension (4*N)
! 179: *> Workspace. Only referenced if one of NCVT, NRU, or NCC is
! 180: *> nonzero, and if N is at least 2.
! 181: *> \endverbatim
! 182: *>
! 183: *> \param[out] INFO
! 184: *> \verbatim
! 185: *> INFO is INTEGER
! 186: *> On exit, a value of 0 indicates a successful exit.
! 187: *> If INFO < 0, argument number -INFO is illegal.
! 188: *> If INFO > 0, the algorithm did not converge, and INFO
! 189: *> specifies how many superdiagonals did not converge.
! 190: *> \endverbatim
! 191: *
! 192: * Authors:
! 193: * ========
! 194: *
! 195: *> \author Univ. of Tennessee
! 196: *> \author Univ. of California Berkeley
! 197: *> \author Univ. of Colorado Denver
! 198: *> \author NAG Ltd.
! 199: *
! 200: *> \date November 2011
! 201: *
! 202: *> \ingroup auxOTHERauxiliary
! 203: *
! 204: *> \par Contributors:
! 205: * ==================
! 206: *>
! 207: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 208: *> California at Berkeley, USA
! 209: *>
! 210: * =====================================================================
1.1 bertrand 211: SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT,
212: $ U, LDU, C, LDC, WORK, INFO )
213: *
1.9 ! bertrand 214: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 215: * -- LAPACK is a software package provided by Univ. of Tennessee, --
216: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 217: * November 2011
1.1 bertrand 218: *
219: * .. Scalar Arguments ..
220: CHARACTER UPLO
221: INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE
222: * ..
223: * .. Array Arguments ..
224: DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ),
225: $ VT( LDVT, * ), WORK( * )
226: * ..
227: *
228: * =====================================================================
229: *
230: * .. Parameters ..
231: DOUBLE PRECISION ZERO
232: PARAMETER ( ZERO = 0.0D+0 )
233: * ..
234: * .. Local Scalars ..
235: LOGICAL ROTATE
236: INTEGER I, ISUB, IUPLO, J, NP1, SQRE1
237: DOUBLE PRECISION CS, R, SMIN, SN
238: * ..
239: * .. External Subroutines ..
240: EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA
241: * ..
242: * .. External Functions ..
243: LOGICAL LSAME
244: EXTERNAL LSAME
245: * ..
246: * .. Intrinsic Functions ..
247: INTRINSIC MAX
248: * ..
249: * .. Executable Statements ..
250: *
251: * Test the input parameters.
252: *
253: INFO = 0
254: IUPLO = 0
255: IF( LSAME( UPLO, 'U' ) )
256: $ IUPLO = 1
257: IF( LSAME( UPLO, 'L' ) )
258: $ IUPLO = 2
259: IF( IUPLO.EQ.0 ) THEN
260: INFO = -1
261: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
262: INFO = -2
263: ELSE IF( N.LT.0 ) THEN
264: INFO = -3
265: ELSE IF( NCVT.LT.0 ) THEN
266: INFO = -4
267: ELSE IF( NRU.LT.0 ) THEN
268: INFO = -5
269: ELSE IF( NCC.LT.0 ) THEN
270: INFO = -6
271: ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR.
272: $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN
273: INFO = -10
274: ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN
275: INFO = -12
276: ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR.
277: $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN
278: INFO = -14
279: END IF
280: IF( INFO.NE.0 ) THEN
281: CALL XERBLA( 'DLASDQ', -INFO )
282: RETURN
283: END IF
284: IF( N.EQ.0 )
285: $ RETURN
286: *
287: * ROTATE is true if any singular vectors desired, false otherwise
288: *
289: ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 )
290: NP1 = N + 1
291: SQRE1 = SQRE
292: *
293: * If matrix non-square upper bidiagonal, rotate to be lower
294: * bidiagonal. The rotations are on the right.
295: *
296: IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN
297: DO 10 I = 1, N - 1
298: CALL DLARTG( D( I ), E( I ), CS, SN, R )
299: D( I ) = R
300: E( I ) = SN*D( I+1 )
301: D( I+1 ) = CS*D( I+1 )
302: IF( ROTATE ) THEN
303: WORK( I ) = CS
304: WORK( N+I ) = SN
305: END IF
306: 10 CONTINUE
307: CALL DLARTG( D( N ), E( N ), CS, SN, R )
308: D( N ) = R
309: E( N ) = ZERO
310: IF( ROTATE ) THEN
311: WORK( N ) = CS
312: WORK( N+N ) = SN
313: END IF
314: IUPLO = 2
315: SQRE1 = 0
316: *
317: * Update singular vectors if desired.
318: *
319: IF( NCVT.GT.0 )
320: $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ),
321: $ WORK( NP1 ), VT, LDVT )
322: END IF
323: *
324: * If matrix lower bidiagonal, rotate to be upper bidiagonal
325: * by applying Givens rotations on the left.
326: *
327: IF( IUPLO.EQ.2 ) THEN
328: DO 20 I = 1, N - 1
329: CALL DLARTG( D( I ), E( I ), CS, SN, R )
330: D( I ) = R
331: E( I ) = SN*D( I+1 )
332: D( I+1 ) = CS*D( I+1 )
333: IF( ROTATE ) THEN
334: WORK( I ) = CS
335: WORK( N+I ) = SN
336: END IF
337: 20 CONTINUE
338: *
339: * If matrix (N+1)-by-N lower bidiagonal, one additional
340: * rotation is needed.
341: *
342: IF( SQRE1.EQ.1 ) THEN
343: CALL DLARTG( D( N ), E( N ), CS, SN, R )
344: D( N ) = R
345: IF( ROTATE ) THEN
346: WORK( N ) = CS
347: WORK( N+N ) = SN
348: END IF
349: END IF
350: *
351: * Update singular vectors if desired.
352: *
353: IF( NRU.GT.0 ) THEN
354: IF( SQRE1.EQ.0 ) THEN
355: CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ),
356: $ WORK( NP1 ), U, LDU )
357: ELSE
358: CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ),
359: $ WORK( NP1 ), U, LDU )
360: END IF
361: END IF
362: IF( NCC.GT.0 ) THEN
363: IF( SQRE1.EQ.0 ) THEN
364: CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ),
365: $ WORK( NP1 ), C, LDC )
366: ELSE
367: CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ),
368: $ WORK( NP1 ), C, LDC )
369: END IF
370: END IF
371: END IF
372: *
373: * Call DBDSQR to compute the SVD of the reduced real
374: * N-by-N upper bidiagonal matrix.
375: *
376: CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C,
377: $ LDC, WORK, INFO )
378: *
379: * Sort the singular values into ascending order (insertion sort on
380: * singular values, but only one transposition per singular vector)
381: *
382: DO 40 I = 1, N
383: *
384: * Scan for smallest D(I).
385: *
386: ISUB = I
387: SMIN = D( I )
388: DO 30 J = I + 1, N
389: IF( D( J ).LT.SMIN ) THEN
390: ISUB = J
391: SMIN = D( J )
392: END IF
393: 30 CONTINUE
394: IF( ISUB.NE.I ) THEN
395: *
396: * Swap singular values and vectors.
397: *
398: D( ISUB ) = D( I )
399: D( I ) = SMIN
400: IF( NCVT.GT.0 )
401: $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT )
402: IF( NRU.GT.0 )
403: $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 )
404: IF( NCC.GT.0 )
405: $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC )
406: END IF
407: 40 CONTINUE
408: *
409: RETURN
410: *
411: * End of DLASDQ
412: *
413: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack