 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 00002 $ RANK, WORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2.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 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER UPLO 00011 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ 00012 DOUBLE PRECISION RCOND 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ) 00016 DOUBLE PRECISION B( LDB, * ), D( * ), E( * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DLALSD uses the singular value decomposition of A to solve the least 00023 * squares problem of finding X to minimize the Euclidean norm of each 00024 * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B 00025 * are N-by-NRHS. The solution X overwrites B. 00026 * 00027 * The singular values of A smaller than RCOND times the largest 00028 * singular value are treated as zero in solving the least squares 00029 * problem; in this case a minimum norm solution is returned. 00030 * The actual singular values are returned in D in ascending order. 00031 * 00032 * This code makes very mild assumptions about floating point 00033 * arithmetic. It will work on machines with a guard digit in 00034 * add/subtract, or on those binary machines without guard digits 00035 * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. 00036 * It could conceivably fail on hexadecimal or decimal machines 00037 * without guard digits, but we know of none. 00038 * 00039 * Arguments 00040 * ========= 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': D and E define an upper bidiagonal matrix. 00044 * = 'L': D and E define a lower bidiagonal matrix. 00045 * 00046 * SMLSIZ (input) INTEGER 00047 * The maximum size of the subproblems at the bottom of the 00048 * computation tree. 00049 * 00050 * N (input) INTEGER 00051 * The dimension of the bidiagonal matrix. N >= 0. 00052 * 00053 * NRHS (input) INTEGER 00054 * The number of columns of B. NRHS must be at least 1. 00055 * 00056 * D (input/output) DOUBLE PRECISION array, dimension (N) 00057 * On entry D contains the main diagonal of the bidiagonal 00058 * matrix. On exit, if INFO = 0, D contains its singular values. 00059 * 00060 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00061 * Contains the super-diagonal entries of the bidiagonal matrix. 00062 * On exit, E has been destroyed. 00063 * 00064 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) 00065 * On input, B contains the right hand sides of the least 00066 * squares problem. On output, B contains the solution X. 00067 * 00068 * LDB (input) INTEGER 00069 * The leading dimension of B in the calling subprogram. 00070 * LDB must be at least max(1,N). 00071 * 00072 * RCOND (input) DOUBLE PRECISION 00073 * The singular values of A less than or equal to RCOND times 00074 * the largest singular value are treated as zero in solving 00075 * the least squares problem. If RCOND is negative, 00076 * machine precision is used instead. 00077 * For example, if diag(S)*X=B were the least squares problem, 00078 * where diag(S) is a diagonal matrix of singular values, the 00079 * solution would be X(i) = B(i) / S(i) if S(i) is greater than 00080 * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to 00081 * RCOND*max(S). 00082 * 00083 * RANK (output) INTEGER 00084 * The number of singular values of A greater than RCOND times 00085 * the largest singular value. 00086 * 00087 * WORK (workspace) DOUBLE PRECISION array, dimension at least 00088 * (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2), 00089 * where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1). 00090 * 00091 * IWORK (workspace) INTEGER array, dimension at least 00092 * (3*N*NLVL + 11*N) 00093 * 00094 * INFO (output) INTEGER 00095 * = 0: successful exit. 00096 * < 0: if INFO = -i, the i-th argument had an illegal value. 00097 * > 0: The algorithm failed to compute a singular value while 00098 * working on the submatrix lying in rows and columns 00099 * INFO/(N+1) through MOD(INFO,N+1). 00100 * 00101 * Further Details 00102 * =============== 00103 * 00104 * Based on contributions by 00105 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00106 * California at Berkeley, USA 00107 * Osni Marques, LBNL/NERSC, USA 00108 * 00109 * ===================================================================== 00110 * 00111 * .. Parameters .. 00112 DOUBLE PRECISION ZERO, ONE, TWO 00113 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00114 * .. 00115 * .. Local Scalars .. 00116 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 00117 $ GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL, 00118 $ NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI, 00119 $ SMLSZP, SQRE, ST, ST1, U, VT, Z 00120 DOUBLE PRECISION CS, EPS, ORGNRM, R, RCND, SN, TOL 00121 * .. 00122 * .. External Functions .. 00123 INTEGER IDAMAX 00124 DOUBLE PRECISION DLAMCH, DLANST 00125 EXTERNAL IDAMAX, DLAMCH, DLANST 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL DCOPY, DGEMM, DLACPY, DLALSA, DLARTG, DLASCL, 00129 $ DLASDA, DLASDQ, DLASET, DLASRT, DROT, XERBLA 00130 * .. 00131 * .. Intrinsic Functions .. 00132 INTRINSIC ABS, DBLE, INT, LOG, SIGN 00133 * .. 00134 * .. Executable Statements .. 00135 * 00136 * Test the input parameters. 00137 * 00138 INFO = 0 00139 * 00140 IF( N.LT.0 ) THEN 00141 INFO = -3 00142 ELSE IF( NRHS.LT.1 ) THEN 00143 INFO = -4 00144 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 00145 INFO = -8 00146 END IF 00147 IF( INFO.NE.0 ) THEN 00148 CALL XERBLA( 'DLALSD', -INFO ) 00149 RETURN 00150 END IF 00151 * 00152 EPS = DLAMCH( 'Epsilon' ) 00153 * 00154 * Set up the tolerance. 00155 * 00156 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 00157 RCND = EPS 00158 ELSE 00159 RCND = RCOND 00160 END IF 00161 * 00162 RANK = 0 00163 * 00164 * Quick return if possible. 00165 * 00166 IF( N.EQ.0 ) THEN 00167 RETURN 00168 ELSE IF( N.EQ.1 ) THEN 00169 IF( D( 1 ).EQ.ZERO ) THEN 00170 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB ) 00171 ELSE 00172 RANK = 1 00173 CALL DLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 00174 D( 1 ) = ABS( D( 1 ) ) 00175 END IF 00176 RETURN 00177 END IF 00178 * 00179 * Rotate the matrix if it is lower bidiagonal. 00180 * 00181 IF( UPLO.EQ.'L' ) THEN 00182 DO 10 I = 1, N - 1 00183 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00184 D( I ) = R 00185 E( I ) = SN*D( I+1 ) 00186 D( I+1 ) = CS*D( I+1 ) 00187 IF( NRHS.EQ.1 ) THEN 00188 CALL DROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 00189 ELSE 00190 WORK( I*2-1 ) = CS 00191 WORK( I*2 ) = SN 00192 END IF 00193 10 CONTINUE 00194 IF( NRHS.GT.1 ) THEN 00195 DO 30 I = 1, NRHS 00196 DO 20 J = 1, N - 1 00197 CS = WORK( J*2-1 ) 00198 SN = WORK( J*2 ) 00199 CALL DROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 00200 20 CONTINUE 00201 30 CONTINUE 00202 END IF 00203 END IF 00204 * 00205 * Scale. 00206 * 00207 NM1 = N - 1 00208 ORGNRM = DLANST( 'M', N, D, E ) 00209 IF( ORGNRM.EQ.ZERO ) THEN 00210 CALL DLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB ) 00211 RETURN 00212 END IF 00213 * 00214 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00215 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 00216 * 00217 * If N is smaller than the minimum divide size SMLSIZ, then solve 00218 * the problem with another solver. 00219 * 00220 IF( N.LE.SMLSIZ ) THEN 00221 NWORK = 1 + N*N 00222 CALL DLASET( 'A', N, N, ZERO, ONE, WORK, N ) 00223 CALL DLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B, 00224 $ LDB, WORK( NWORK ), INFO ) 00225 IF( INFO.NE.0 ) THEN 00226 RETURN 00227 END IF 00228 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00229 DO 40 I = 1, N 00230 IF( D( I ).LE.TOL ) THEN 00231 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB ) 00232 ELSE 00233 CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 00234 $ LDB, INFO ) 00235 RANK = RANK + 1 00236 END IF 00237 40 CONTINUE 00238 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 00239 $ WORK( NWORK ), N ) 00240 CALL DLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB ) 00241 * 00242 * Unscale. 00243 * 00244 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00245 CALL DLASRT( 'D', N, D, INFO ) 00246 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00247 * 00248 RETURN 00249 END IF 00250 * 00251 * Book-keeping and setting up some constants. 00252 * 00253 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 00254 * 00255 SMLSZP = SMLSIZ + 1 00256 * 00257 U = 1 00258 VT = 1 + SMLSIZ*N 00259 DIFL = VT + SMLSZP*N 00260 DIFR = DIFL + NLVL*N 00261 Z = DIFR + NLVL*N*2 00262 C = Z + NLVL*N 00263 S = C + N 00264 POLES = S + N 00265 GIVNUM = POLES + 2*NLVL*N 00266 BX = GIVNUM + 2*NLVL*N 00267 NWORK = BX + N*NRHS 00268 * 00269 SIZEI = 1 + N 00270 K = SIZEI + N 00271 GIVPTR = K + N 00272 PERM = GIVPTR + N 00273 GIVCOL = PERM + NLVL*N 00274 IWK = GIVCOL + NLVL*N*2 00275 * 00276 ST = 1 00277 SQRE = 0 00278 ICMPQ1 = 1 00279 ICMPQ2 = 0 00280 NSUB = 0 00281 * 00282 DO 50 I = 1, N 00283 IF( ABS( D( I ) ).LT.EPS ) THEN 00284 D( I ) = SIGN( EPS, D( I ) ) 00285 END IF 00286 50 CONTINUE 00287 * 00288 DO 60 I = 1, NM1 00289 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 00290 NSUB = NSUB + 1 00291 IWORK( NSUB ) = ST 00292 * 00293 * Subproblem found. First determine its size and then 00294 * apply divide and conquer on it. 00295 * 00296 IF( I.LT.NM1 ) THEN 00297 * 00298 * A subproblem with E(I) small for I < NM1. 00299 * 00300 NSIZE = I - ST + 1 00301 IWORK( SIZEI+NSUB-1 ) = NSIZE 00302 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 00303 * 00304 * A subproblem with E(NM1) not too small but I = NM1. 00305 * 00306 NSIZE = N - ST + 1 00307 IWORK( SIZEI+NSUB-1 ) = NSIZE 00308 ELSE 00309 * 00310 * A subproblem with E(NM1) small. This implies an 00311 * 1-by-1 subproblem at D(N), which is not solved 00312 * explicitly. 00313 * 00314 NSIZE = I - ST + 1 00315 IWORK( SIZEI+NSUB-1 ) = NSIZE 00316 NSUB = NSUB + 1 00317 IWORK( NSUB ) = N 00318 IWORK( SIZEI+NSUB-1 ) = 1 00319 CALL DCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 00320 END IF 00321 ST1 = ST - 1 00322 IF( NSIZE.EQ.1 ) THEN 00323 * 00324 * This is a 1-by-1 subproblem and is not solved 00325 * explicitly. 00326 * 00327 CALL DCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 00328 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00329 * 00330 * This is a small subproblem and is solved by DLASDQ. 00331 * 00332 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 00333 $ WORK( VT+ST1 ), N ) 00334 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ), 00335 $ E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ), 00336 $ N, B( ST, 1 ), LDB, WORK( NWORK ), INFO ) 00337 IF( INFO.NE.0 ) THEN 00338 RETURN 00339 END IF 00340 CALL DLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 00341 $ WORK( BX+ST1 ), N ) 00342 ELSE 00343 * 00344 * A large problem. Solve it using divide and conquer. 00345 * 00346 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 00347 $ E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ), 00348 $ IWORK( K+ST1 ), WORK( DIFL+ST1 ), 00349 $ WORK( DIFR+ST1 ), WORK( Z+ST1 ), 00350 $ WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 00351 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 00352 $ WORK( GIVNUM+ST1 ), WORK( C+ST1 ), 00353 $ WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ), 00354 $ INFO ) 00355 IF( INFO.NE.0 ) THEN 00356 RETURN 00357 END IF 00358 BXST = BX + ST1 00359 CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 00360 $ LDB, WORK( BXST ), N, WORK( U+ST1 ), N, 00361 $ WORK( VT+ST1 ), IWORK( K+ST1 ), 00362 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ), 00363 $ WORK( Z+ST1 ), WORK( POLES+ST1 ), 00364 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00365 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ), 00366 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ), 00367 $ IWORK( IWK ), INFO ) 00368 IF( INFO.NE.0 ) THEN 00369 RETURN 00370 END IF 00371 END IF 00372 ST = I + 1 00373 END IF 00374 60 CONTINUE 00375 * 00376 * Apply the singular values and treat the tiny ones as zero. 00377 * 00378 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00379 * 00380 DO 70 I = 1, N 00381 * 00382 * Some of the elements in D can be negative because 1-by-1 00383 * subproblems were not solved explicitly. 00384 * 00385 IF( ABS( D( I ) ).LE.TOL ) THEN 00386 CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N ) 00387 ELSE 00388 RANK = RANK + 1 00389 CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 00390 $ WORK( BX+I-1 ), N, INFO ) 00391 END IF 00392 D( I ) = ABS( D( I ) ) 00393 70 CONTINUE 00394 * 00395 * Now apply back the right singular vectors. 00396 * 00397 ICMPQ2 = 1 00398 DO 80 I = 1, NSUB 00399 ST = IWORK( I ) 00400 ST1 = ST - 1 00401 NSIZE = IWORK( SIZEI+I-1 ) 00402 BXST = BX + ST1 00403 IF( NSIZE.EQ.1 ) THEN 00404 CALL DCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 00405 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00406 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00407 $ WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO, 00408 $ B( ST, 1 ), LDB ) 00409 ELSE 00410 CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 00411 $ B( ST, 1 ), LDB, WORK( U+ST1 ), N, 00412 $ WORK( VT+ST1 ), IWORK( K+ST1 ), 00413 $ WORK( DIFL+ST1 ), WORK( DIFR+ST1 ), 00414 $ WORK( Z+ST1 ), WORK( POLES+ST1 ), 00415 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00416 $ IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ), 00417 $ WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ), 00418 $ IWORK( IWK ), INFO ) 00419 IF( INFO.NE.0 ) THEN 00420 RETURN 00421 END IF 00422 END IF 00423 80 CONTINUE 00424 * 00425 * Unscale and sort the singular values. 00426 * 00427 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00428 CALL DLASRT( 'D', N, D, INFO ) 00429 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00430 * 00431 RETURN 00432 * 00433 * End of DLALSD 00434 * 00435 END
1.7.3 
