 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, 00002 $ WORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER COMPQ, UPLO 00011 INTEGER INFO, LDU, LDVT, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IQ( * ), IWORK( * ) 00015 DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ), 00016 $ VT( LDVT, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DBDSDC computes the singular value decomposition (SVD) of a real 00023 * N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, 00024 * using a divide and conquer method, where S is a diagonal matrix 00025 * with non-negative diagonal elements (the singular values of B), and 00026 * U and VT are orthogonal matrices of left and right singular vectors, 00027 * respectively. DBDSDC can be used to compute all singular values, 00028 * and optionally, singular vectors or singular vectors in compact form. 00029 * 00030 * This code makes very mild assumptions about floating point 00031 * arithmetic. It will work on machines with a guard digit in 00032 * add/subtract, or on those binary machines without guard digits 00033 * which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. 00034 * It could conceivably fail on hexadecimal or decimal machines 00035 * without guard digits, but we know of none. See DLASD3 for details. 00036 * 00037 * The code currently calls DLASDQ if singular values only are desired. 00038 * However, it can be slightly modified to compute singular values 00039 * using the divide and conquer method. 00040 * 00041 * Arguments 00042 * ========= 00043 * 00044 * UPLO (input) CHARACTER*1 00045 * = 'U': B is upper bidiagonal. 00046 * = 'L': B is lower bidiagonal. 00047 * 00048 * COMPQ (input) CHARACTER*1 00049 * Specifies whether singular vectors are to be computed 00050 * as follows: 00051 * = 'N': Compute singular values only; 00052 * = 'P': Compute singular values and compute singular 00053 * vectors in compact form; 00054 * = 'I': Compute singular values and singular vectors. 00055 * 00056 * N (input) INTEGER 00057 * The order of the matrix B. N >= 0. 00058 * 00059 * D (input/output) DOUBLE PRECISION array, dimension (N) 00060 * On entry, the n diagonal elements of the bidiagonal matrix B. 00061 * On exit, if INFO=0, the singular values of B. 00062 * 00063 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00064 * On entry, the elements of E contain the offdiagonal 00065 * elements of the bidiagonal matrix whose SVD is desired. 00066 * On exit, E has been destroyed. 00067 * 00068 * U (output) DOUBLE PRECISION array, dimension (LDU,N) 00069 * If COMPQ = 'I', then: 00070 * On exit, if INFO = 0, U contains the left singular vectors 00071 * of the bidiagonal matrix. 00072 * For other values of COMPQ, U is not referenced. 00073 * 00074 * LDU (input) INTEGER 00075 * The leading dimension of the array U. LDU >= 1. 00076 * If singular vectors are desired, then LDU >= max( 1, N ). 00077 * 00078 * VT (output) DOUBLE PRECISION array, dimension (LDVT,N) 00079 * If COMPQ = 'I', then: 00080 * On exit, if INFO = 0, VT**T contains the right singular 00081 * vectors of the bidiagonal matrix. 00082 * For other values of COMPQ, VT is not referenced. 00083 * 00084 * LDVT (input) INTEGER 00085 * The leading dimension of the array VT. LDVT >= 1. 00086 * If singular vectors are desired, then LDVT >= max( 1, N ). 00087 * 00088 * Q (output) DOUBLE PRECISION array, dimension (LDQ) 00089 * If COMPQ = 'P', then: 00090 * On exit, if INFO = 0, Q and IQ contain the left 00091 * and right singular vectors in a compact form, 00092 * requiring O(N log N) space instead of 2*N**2. 00093 * In particular, Q contains all the DOUBLE PRECISION data in 00094 * LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) 00095 * words of memory, where SMLSIZ is returned by ILAENV and 00096 * is equal to the maximum size of the subproblems at the 00097 * bottom of the computation tree (usually about 25). 00098 * For other values of COMPQ, Q is not referenced. 00099 * 00100 * IQ (output) INTEGER array, dimension (LDIQ) 00101 * If COMPQ = 'P', then: 00102 * On exit, if INFO = 0, Q and IQ contain the left 00103 * and right singular vectors in a compact form, 00104 * requiring O(N log N) space instead of 2*N**2. 00105 * In particular, IQ contains all INTEGER data in 00106 * LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) 00107 * words of memory, where SMLSIZ is returned by ILAENV and 00108 * is equal to the maximum size of the subproblems at the 00109 * bottom of the computation tree (usually about 25). 00110 * For other values of COMPQ, IQ is not referenced. 00111 * 00112 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00113 * If COMPQ = 'N' then LWORK >= (4 * N). 00114 * If COMPQ = 'P' then LWORK >= (6 * N). 00115 * If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). 00116 * 00117 * IWORK (workspace) INTEGER array, dimension (8*N) 00118 * 00119 * INFO (output) INTEGER 00120 * = 0: successful exit. 00121 * < 0: if INFO = -i, the i-th argument had an illegal value. 00122 * > 0: The algorithm failed to compute a singular value. 00123 * The update process of divide and conquer failed. 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 * Changed dimension statement in comment describing E from (N) to 00134 * (N-1). Sven, 17 Feb 05. 00135 * ===================================================================== 00136 * 00137 * .. Parameters .. 00138 DOUBLE PRECISION ZERO, ONE, TWO 00139 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00140 * .. 00141 * .. Local Scalars .. 00142 INTEGER DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC, 00143 $ ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK, 00144 $ MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ, 00145 $ SMLSZP, SQRE, START, WSTART, Z 00146 DOUBLE PRECISION CS, EPS, ORGNRM, P, R, SN 00147 * .. 00148 * .. External Functions .. 00149 LOGICAL LSAME 00150 INTEGER ILAENV 00151 DOUBLE PRECISION DLAMCH, DLANST 00152 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST 00153 * .. 00154 * .. External Subroutines .. 00155 EXTERNAL DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ, 00156 $ DLASET, DLASR, DSWAP, XERBLA 00157 * .. 00158 * .. Intrinsic Functions .. 00159 INTRINSIC ABS, DBLE, INT, LOG, SIGN 00160 * .. 00161 * .. Executable Statements .. 00162 * 00163 * Test the input parameters. 00164 * 00165 INFO = 0 00166 * 00167 IUPLO = 0 00168 IF( LSAME( UPLO, 'U' ) ) 00169 $ IUPLO = 1 00170 IF( LSAME( UPLO, 'L' ) ) 00171 $ IUPLO = 2 00172 IF( LSAME( COMPQ, 'N' ) ) THEN 00173 ICOMPQ = 0 00174 ELSE IF( LSAME( COMPQ, 'P' ) ) THEN 00175 ICOMPQ = 1 00176 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN 00177 ICOMPQ = 2 00178 ELSE 00179 ICOMPQ = -1 00180 END IF 00181 IF( IUPLO.EQ.0 ) THEN 00182 INFO = -1 00183 ELSE IF( ICOMPQ.LT.0 ) THEN 00184 INFO = -2 00185 ELSE IF( N.LT.0 ) THEN 00186 INFO = -3 00187 ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT. 00188 $ N ) ) ) THEN 00189 INFO = -7 00190 ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT. 00191 $ N ) ) ) THEN 00192 INFO = -9 00193 END IF 00194 IF( INFO.NE.0 ) THEN 00195 CALL XERBLA( 'DBDSDC', -INFO ) 00196 RETURN 00197 END IF 00198 * 00199 * Quick return if possible 00200 * 00201 IF( N.EQ.0 ) 00202 $ RETURN 00203 SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 ) 00204 IF( N.EQ.1 ) THEN 00205 IF( ICOMPQ.EQ.1 ) THEN 00206 Q( 1 ) = SIGN( ONE, D( 1 ) ) 00207 Q( 1+SMLSIZ*N ) = ONE 00208 ELSE IF( ICOMPQ.EQ.2 ) THEN 00209 U( 1, 1 ) = SIGN( ONE, D( 1 ) ) 00210 VT( 1, 1 ) = ONE 00211 END IF 00212 D( 1 ) = ABS( D( 1 ) ) 00213 RETURN 00214 END IF 00215 NM1 = N - 1 00216 * 00217 * If matrix lower bidiagonal, rotate to be upper bidiagonal 00218 * by applying Givens rotations on the left 00219 * 00220 WSTART = 1 00221 QSTART = 3 00222 IF( ICOMPQ.EQ.1 ) THEN 00223 CALL DCOPY( N, D, 1, Q( 1 ), 1 ) 00224 CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 ) 00225 END IF 00226 IF( IUPLO.EQ.2 ) THEN 00227 QSTART = 5 00228 WSTART = 2*N - 1 00229 DO 10 I = 1, N - 1 00230 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00231 D( I ) = R 00232 E( I ) = SN*D( I+1 ) 00233 D( I+1 ) = CS*D( I+1 ) 00234 IF( ICOMPQ.EQ.1 ) THEN 00235 Q( I+2*N ) = CS 00236 Q( I+3*N ) = SN 00237 ELSE IF( ICOMPQ.EQ.2 ) THEN 00238 WORK( I ) = CS 00239 WORK( NM1+I ) = -SN 00240 END IF 00241 10 CONTINUE 00242 END IF 00243 * 00244 * If ICOMPQ = 0, use DLASDQ to compute the singular values. 00245 * 00246 IF( ICOMPQ.EQ.0 ) THEN 00247 CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U, 00248 $ LDU, WORK( WSTART ), INFO ) 00249 GO TO 40 00250 END IF 00251 * 00252 * If N is smaller than the minimum divide size SMLSIZ, then solve 00253 * the problem with another solver. 00254 * 00255 IF( N.LE.SMLSIZ ) THEN 00256 IF( ICOMPQ.EQ.2 ) THEN 00257 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU ) 00258 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 00259 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U, 00260 $ LDU, WORK( WSTART ), INFO ) 00261 ELSE IF( ICOMPQ.EQ.1 ) THEN 00262 IU = 1 00263 IVT = IU + N 00264 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ), 00265 $ N ) 00266 CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ), 00267 $ N ) 00268 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, 00269 $ Q( IVT+( QSTART-1 )*N ), N, 00270 $ Q( IU+( QSTART-1 )*N ), N, 00271 $ Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ), 00272 $ INFO ) 00273 END IF 00274 GO TO 40 00275 END IF 00276 * 00277 IF( ICOMPQ.EQ.2 ) THEN 00278 CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU ) 00279 CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT ) 00280 END IF 00281 * 00282 * Scale. 00283 * 00284 ORGNRM = DLANST( 'M', N, D, E ) 00285 IF( ORGNRM.EQ.ZERO ) 00286 $ RETURN 00287 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR ) 00288 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR ) 00289 * 00290 EPS = (0.9D+0)*DLAMCH( 'Epsilon' ) 00291 * 00292 MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 00293 SMLSZP = SMLSIZ + 1 00294 * 00295 IF( ICOMPQ.EQ.1 ) THEN 00296 IU = 1 00297 IVT = 1 + SMLSIZ 00298 DIFL = IVT + SMLSZP 00299 DIFR = DIFL + MLVL 00300 Z = DIFR + MLVL*2 00301 IC = Z + MLVL 00302 IS = IC + 1 00303 POLES = IS + 1 00304 GIVNUM = POLES + 2*MLVL 00305 * 00306 K = 1 00307 GIVPTR = 2 00308 PERM = 3 00309 GIVCOL = PERM + MLVL 00310 END IF 00311 * 00312 DO 20 I = 1, N 00313 IF( ABS( D( I ) ).LT.EPS ) THEN 00314 D( I ) = SIGN( EPS, D( I ) ) 00315 END IF 00316 20 CONTINUE 00317 * 00318 START = 1 00319 SQRE = 0 00320 * 00321 DO 30 I = 1, NM1 00322 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 00323 * 00324 * Subproblem found. First determine its size and then 00325 * apply divide and conquer on it. 00326 * 00327 IF( I.LT.NM1 ) THEN 00328 * 00329 * A subproblem with E(I) small for I < NM1. 00330 * 00331 NSIZE = I - START + 1 00332 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 00333 * 00334 * A subproblem with E(NM1) not too small but I = NM1. 00335 * 00336 NSIZE = N - START + 1 00337 ELSE 00338 * 00339 * A subproblem with E(NM1) small. This implies an 00340 * 1-by-1 subproblem at D(N). Solve this 1-by-1 problem 00341 * first. 00342 * 00343 NSIZE = I - START + 1 00344 IF( ICOMPQ.EQ.2 ) THEN 00345 U( N, N ) = SIGN( ONE, D( N ) ) 00346 VT( N, N ) = ONE 00347 ELSE IF( ICOMPQ.EQ.1 ) THEN 00348 Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) ) 00349 Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE 00350 END IF 00351 D( N ) = ABS( D( N ) ) 00352 END IF 00353 IF( ICOMPQ.EQ.2 ) THEN 00354 CALL DLASD0( NSIZE, SQRE, D( START ), E( START ), 00355 $ U( START, START ), LDU, VT( START, START ), 00356 $ LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO ) 00357 ELSE 00358 CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ), 00359 $ E( START ), Q( START+( IU+QSTART-2 )*N ), N, 00360 $ Q( START+( IVT+QSTART-2 )*N ), 00361 $ IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )* 00362 $ N ), Q( START+( DIFR+QSTART-2 )*N ), 00363 $ Q( START+( Z+QSTART-2 )*N ), 00364 $ Q( START+( POLES+QSTART-2 )*N ), 00365 $ IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ), 00366 $ N, IQ( START+PERM*N ), 00367 $ Q( START+( GIVNUM+QSTART-2 )*N ), 00368 $ Q( START+( IC+QSTART-2 )*N ), 00369 $ Q( START+( IS+QSTART-2 )*N ), 00370 $ WORK( WSTART ), IWORK, INFO ) 00371 END IF 00372 IF( INFO.NE.0 ) THEN 00373 RETURN 00374 END IF 00375 START = I + 1 00376 END IF 00377 30 CONTINUE 00378 * 00379 * Unscale 00380 * 00381 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR ) 00382 40 CONTINUE 00383 * 00384 * Use Selection Sort to minimize swaps of singular vectors 00385 * 00386 DO 60 II = 2, N 00387 I = II - 1 00388 KK = I 00389 P = D( I ) 00390 DO 50 J = II, N 00391 IF( D( J ).GT.P ) THEN 00392 KK = J 00393 P = D( J ) 00394 END IF 00395 50 CONTINUE 00396 IF( KK.NE.I ) THEN 00397 D( KK ) = D( I ) 00398 D( I ) = P 00399 IF( ICOMPQ.EQ.1 ) THEN 00400 IQ( I ) = KK 00401 ELSE IF( ICOMPQ.EQ.2 ) THEN 00402 CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 ) 00403 CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT ) 00404 END IF 00405 ELSE IF( ICOMPQ.EQ.1 ) THEN 00406 IQ( I ) = I 00407 END IF 00408 60 CONTINUE 00409 * 00410 * If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO 00411 * 00412 IF( ICOMPQ.EQ.1 ) THEN 00413 IF( IUPLO.EQ.1 ) THEN 00414 IQ( N ) = 1 00415 ELSE 00416 IQ( N ) = 0 00417 END IF 00418 END IF 00419 * 00420 * If B is lower bidiagonal, update U by those Givens rotations 00421 * which rotated B to be upper bidiagonal 00422 * 00423 IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) ) 00424 $ CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU ) 00425 * 00426 RETURN 00427 * 00428 * End of DBDSDC 00429 * 00430 END
1.7.3 
