 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASDQ( UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, 00002 $ U, LDU, C, LDC, WORK, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO 00011 INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU, SQRE 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), 00015 $ VT( LDVT, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLASDQ computes the singular value decomposition (SVD) of a real 00022 * (upper or lower) bidiagonal matrix with diagonal D and offdiagonal 00023 * E, accumulating the transformations if desired. Letting B denote 00024 * the input bidiagonal matrix, the algorithm computes orthogonal 00025 * matrices Q and P such that B = Q * S * P**T (P**T denotes the transpose 00026 * of P). The singular values S are overwritten on D. 00027 * 00028 * The input matrix U is changed to U * Q if desired. 00029 * The input matrix VT is changed to P**T * VT if desired. 00030 * The input matrix C is changed to Q**T * C if desired. 00031 * 00032 * See "Computing Small Singular Values of Bidiagonal Matrices With 00033 * Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, 00034 * LAPACK Working Note #3, for a detailed description of the algorithm. 00035 * 00036 * Arguments 00037 * ========= 00038 * 00039 * UPLO (input) CHARACTER*1 00040 * On entry, UPLO specifies whether the input bidiagonal matrix 00041 * is upper or lower bidiagonal, and wether it is square are 00042 * not. 00043 * UPLO = 'U' or 'u' B is upper bidiagonal. 00044 * UPLO = 'L' or 'l' B is lower bidiagonal. 00045 * 00046 * SQRE (input) INTEGER 00047 * = 0: then the input matrix is N-by-N. 00048 * = 1: then the input matrix is N-by-(N+1) if UPLU = 'U' and 00049 * (N+1)-by-N if UPLU = 'L'. 00050 * 00051 * The bidiagonal matrix has 00052 * N = NL + NR + 1 rows and 00053 * M = N + SQRE >= N columns. 00054 * 00055 * N (input) INTEGER 00056 * On entry, N specifies the number of rows and columns 00057 * in the matrix. N must be at least 0. 00058 * 00059 * NCVT (input) INTEGER 00060 * On entry, NCVT specifies the number of columns of 00061 * the matrix VT. NCVT must be at least 0. 00062 * 00063 * NRU (input) INTEGER 00064 * On entry, NRU specifies the number of rows of 00065 * the matrix U. NRU must be at least 0. 00066 * 00067 * NCC (input) INTEGER 00068 * On entry, NCC specifies the number of columns of 00069 * the matrix C. NCC must be at least 0. 00070 * 00071 * D (input/output) DOUBLE PRECISION array, dimension (N) 00072 * On entry, D contains the diagonal entries of the 00073 * bidiagonal matrix whose SVD is desired. On normal exit, 00074 * D contains the singular values in ascending order. 00075 * 00076 * E (input/output) DOUBLE PRECISION array. 00077 * dimension is (N-1) if SQRE = 0 and N if SQRE = 1. 00078 * On entry, the entries of E contain the offdiagonal entries 00079 * of the bidiagonal matrix whose SVD is desired. On normal 00080 * exit, E will contain 0. If the algorithm does not converge, 00081 * D and E will contain the diagonal and superdiagonal entries 00082 * of a bidiagonal matrix orthogonally equivalent to the one 00083 * given as input. 00084 * 00085 * VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) 00086 * On entry, contains a matrix which on exit has been 00087 * premultiplied by P**T, dimension N-by-NCVT if SQRE = 0 00088 * and (N+1)-by-NCVT if SQRE = 1 (not referenced if NCVT=0). 00089 * 00090 * LDVT (input) INTEGER 00091 * On entry, LDVT specifies the leading dimension of VT as 00092 * declared in the calling (sub) program. LDVT must be at 00093 * least 1. If NCVT is nonzero LDVT must also be at least N. 00094 * 00095 * U (input/output) DOUBLE PRECISION array, dimension (LDU, N) 00096 * On entry, contains a matrix which on exit has been 00097 * postmultiplied by Q, dimension NRU-by-N if SQRE = 0 00098 * and NRU-by-(N+1) if SQRE = 1 (not referenced if NRU=0). 00099 * 00100 * LDU (input) INTEGER 00101 * On entry, LDU specifies the leading dimension of U as 00102 * declared in the calling (sub) program. LDU must be at 00103 * least max( 1, NRU ) . 00104 * 00105 * C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) 00106 * On entry, contains an N-by-NCC matrix which on exit 00107 * has been premultiplied by Q**T dimension N-by-NCC if SQRE = 0 00108 * and (N+1)-by-NCC if SQRE = 1 (not referenced if NCC=0). 00109 * 00110 * LDC (input) INTEGER 00111 * On entry, LDC specifies the leading dimension of C as 00112 * declared in the calling (sub) program. LDC must be at 00113 * least 1. If NCC is nonzero, LDC must also be at least N. 00114 * 00115 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) 00116 * Workspace. Only referenced if one of NCVT, NRU, or NCC is 00117 * nonzero, and if N is at least 2. 00118 * 00119 * INFO (output) INTEGER 00120 * On exit, a value of 0 indicates a successful exit. 00121 * If INFO < 0, argument number -INFO is illegal. 00122 * If INFO > 0, the algorithm did not converge, and INFO 00123 * specifies how many superdiagonals did not converge. 00124 * 00125 * Further Details 00126 * =============== 00127 * 00128 * Based on contributions by 00129 * Ming Gu and Huan Ren, Computer Science Division, University of 00130 * California at Berkeley, USA 00131 * 00132 * ===================================================================== 00133 * 00134 * .. Parameters .. 00135 DOUBLE PRECISION ZERO 00136 PARAMETER ( ZERO = 0.0D+0 ) 00137 * .. 00138 * .. Local Scalars .. 00139 LOGICAL ROTATE 00140 INTEGER I, ISUB, IUPLO, J, NP1, SQRE1 00141 DOUBLE PRECISION CS, R, SMIN, SN 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL DBDSQR, DLARTG, DLASR, DSWAP, XERBLA 00145 * .. 00146 * .. External Functions .. 00147 LOGICAL LSAME 00148 EXTERNAL LSAME 00149 * .. 00150 * .. Intrinsic Functions .. 00151 INTRINSIC MAX 00152 * .. 00153 * .. Executable Statements .. 00154 * 00155 * Test the input parameters. 00156 * 00157 INFO = 0 00158 IUPLO = 0 00159 IF( LSAME( UPLO, 'U' ) ) 00160 $ IUPLO = 1 00161 IF( LSAME( UPLO, 'L' ) ) 00162 $ IUPLO = 2 00163 IF( IUPLO.EQ.0 ) THEN 00164 INFO = -1 00165 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00166 INFO = -2 00167 ELSE IF( N.LT.0 ) THEN 00168 INFO = -3 00169 ELSE IF( NCVT.LT.0 ) THEN 00170 INFO = -4 00171 ELSE IF( NRU.LT.0 ) THEN 00172 INFO = -5 00173 ELSE IF( NCC.LT.0 ) THEN 00174 INFO = -6 00175 ELSE IF( ( NCVT.EQ.0 .AND. LDVT.LT.1 ) .OR. 00176 $ ( NCVT.GT.0 .AND. LDVT.LT.MAX( 1, N ) ) ) THEN 00177 INFO = -10 00178 ELSE IF( LDU.LT.MAX( 1, NRU ) ) THEN 00179 INFO = -12 00180 ELSE IF( ( NCC.EQ.0 .AND. LDC.LT.1 ) .OR. 00181 $ ( NCC.GT.0 .AND. LDC.LT.MAX( 1, N ) ) ) THEN 00182 INFO = -14 00183 END IF 00184 IF( INFO.NE.0 ) THEN 00185 CALL XERBLA( 'DLASDQ', -INFO ) 00186 RETURN 00187 END IF 00188 IF( N.EQ.0 ) 00189 $ RETURN 00190 * 00191 * ROTATE is true if any singular vectors desired, false otherwise 00192 * 00193 ROTATE = ( NCVT.GT.0 ) .OR. ( NRU.GT.0 ) .OR. ( NCC.GT.0 ) 00194 NP1 = N + 1 00195 SQRE1 = SQRE 00196 * 00197 * If matrix non-square upper bidiagonal, rotate to be lower 00198 * bidiagonal. The rotations are on the right. 00199 * 00200 IF( ( IUPLO.EQ.1 ) .AND. ( SQRE1.EQ.1 ) ) THEN 00201 DO 10 I = 1, N - 1 00202 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00203 D( I ) = R 00204 E( I ) = SN*D( I+1 ) 00205 D( I+1 ) = CS*D( I+1 ) 00206 IF( ROTATE ) THEN 00207 WORK( I ) = CS 00208 WORK( N+I ) = SN 00209 END IF 00210 10 CONTINUE 00211 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 00212 D( N ) = R 00213 E( N ) = ZERO 00214 IF( ROTATE ) THEN 00215 WORK( N ) = CS 00216 WORK( N+N ) = SN 00217 END IF 00218 IUPLO = 2 00219 SQRE1 = 0 00220 * 00221 * Update singular vectors if desired. 00222 * 00223 IF( NCVT.GT.0 ) 00224 $ CALL DLASR( 'L', 'V', 'F', NP1, NCVT, WORK( 1 ), 00225 $ WORK( NP1 ), VT, LDVT ) 00226 END IF 00227 * 00228 * If matrix lower bidiagonal, rotate to be upper bidiagonal 00229 * by applying Givens rotations on the left. 00230 * 00231 IF( IUPLO.EQ.2 ) THEN 00232 DO 20 I = 1, N - 1 00233 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00234 D( I ) = R 00235 E( I ) = SN*D( I+1 ) 00236 D( I+1 ) = CS*D( I+1 ) 00237 IF( ROTATE ) THEN 00238 WORK( I ) = CS 00239 WORK( N+I ) = SN 00240 END IF 00241 20 CONTINUE 00242 * 00243 * If matrix (N+1)-by-N lower bidiagonal, one additional 00244 * rotation is needed. 00245 * 00246 IF( SQRE1.EQ.1 ) THEN 00247 CALL DLARTG( D( N ), E( N ), CS, SN, R ) 00248 D( N ) = R 00249 IF( ROTATE ) THEN 00250 WORK( N ) = CS 00251 WORK( N+N ) = SN 00252 END IF 00253 END IF 00254 * 00255 * Update singular vectors if desired. 00256 * 00257 IF( NRU.GT.0 ) THEN 00258 IF( SQRE1.EQ.0 ) THEN 00259 CALL DLASR( 'R', 'V', 'F', NRU, N, WORK( 1 ), 00260 $ WORK( NP1 ), U, LDU ) 00261 ELSE 00262 CALL DLASR( 'R', 'V', 'F', NRU, NP1, WORK( 1 ), 00263 $ WORK( NP1 ), U, LDU ) 00264 END IF 00265 END IF 00266 IF( NCC.GT.0 ) THEN 00267 IF( SQRE1.EQ.0 ) THEN 00268 CALL DLASR( 'L', 'V', 'F', N, NCC, WORK( 1 ), 00269 $ WORK( NP1 ), C, LDC ) 00270 ELSE 00271 CALL DLASR( 'L', 'V', 'F', NP1, NCC, WORK( 1 ), 00272 $ WORK( NP1 ), C, LDC ) 00273 END IF 00274 END IF 00275 END IF 00276 * 00277 * Call DBDSQR to compute the SVD of the reduced real 00278 * N-by-N upper bidiagonal matrix. 00279 * 00280 CALL DBDSQR( 'U', N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, 00281 $ LDC, WORK, INFO ) 00282 * 00283 * Sort the singular values into ascending order (insertion sort on 00284 * singular values, but only one transposition per singular vector) 00285 * 00286 DO 40 I = 1, N 00287 * 00288 * Scan for smallest D(I). 00289 * 00290 ISUB = I 00291 SMIN = D( I ) 00292 DO 30 J = I + 1, N 00293 IF( D( J ).LT.SMIN ) THEN 00294 ISUB = J 00295 SMIN = D( J ) 00296 END IF 00297 30 CONTINUE 00298 IF( ISUB.NE.I ) THEN 00299 * 00300 * Swap singular values and vectors. 00301 * 00302 D( ISUB ) = D( I ) 00303 D( I ) = SMIN 00304 IF( NCVT.GT.0 ) 00305 $ CALL DSWAP( NCVT, VT( ISUB, 1 ), LDVT, VT( I, 1 ), LDVT ) 00306 IF( NRU.GT.0 ) 00307 $ CALL DSWAP( NRU, U( 1, ISUB ), 1, U( 1, I ), 1 ) 00308 IF( NCC.GT.0 ) 00309 $ CALL DSWAP( NCC, C( ISUB, 1 ), LDC, C( I, 1 ), LDC ) 00310 END IF 00311 40 CONTINUE 00312 * 00313 RETURN 00314 * 00315 * End of DLASDQ 00316 * 00317 END
1.7.3 
