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