 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 00002 $ RANK, WORK, RWORK, 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 UPLO 00011 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ 00012 DOUBLE PRECISION RCOND 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ) 00016 DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 00017 COMPLEX*16 B( LDB, * ), WORK( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * ZLALSD uses the singular value decomposition of A to solve the least 00024 * squares problem of finding X to minimize the Euclidean norm of each 00025 * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B 00026 * are N-by-NRHS. The solution X overwrites B. 00027 * 00028 * The singular values of A smaller than RCOND times the largest 00029 * singular value are treated as zero in solving the least squares 00030 * problem; in this case a minimum norm solution is returned. 00031 * The actual singular values are returned in D in ascending order. 00032 * 00033 * This code makes very mild assumptions about floating point 00034 * arithmetic. It will work on machines with a guard digit in 00035 * add/subtract, or on those binary machines without guard digits 00036 * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. 00037 * It could conceivably fail on hexadecimal or decimal machines 00038 * without guard digits, but we know of none. 00039 * 00040 * Arguments 00041 * ========= 00042 * 00043 * UPLO (input) CHARACTER*1 00044 * = 'U': D and E define an upper bidiagonal matrix. 00045 * = 'L': D and E define a lower bidiagonal matrix. 00046 * 00047 * SMLSIZ (input) INTEGER 00048 * The maximum size of the subproblems at the bottom of the 00049 * computation tree. 00050 * 00051 * N (input) INTEGER 00052 * The dimension of the bidiagonal matrix. N >= 0. 00053 * 00054 * NRHS (input) INTEGER 00055 * The number of columns of B. NRHS must be at least 1. 00056 * 00057 * D (input/output) DOUBLE PRECISION array, dimension (N) 00058 * On entry D contains the main diagonal of the bidiagonal 00059 * matrix. On exit, if INFO = 0, D contains its singular values. 00060 * 00061 * E (input/output) DOUBLE PRECISION array, dimension (N-1) 00062 * Contains the super-diagonal entries of the bidiagonal matrix. 00063 * On exit, E has been destroyed. 00064 * 00065 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 00066 * On input, B contains the right hand sides of the least 00067 * squares problem. On output, B contains the solution X. 00068 * 00069 * LDB (input) INTEGER 00070 * The leading dimension of B in the calling subprogram. 00071 * LDB must be at least max(1,N). 00072 * 00073 * RCOND (input) DOUBLE PRECISION 00074 * The singular values of A less than or equal to RCOND times 00075 * the largest singular value are treated as zero in solving 00076 * the least squares problem. If RCOND is negative, 00077 * machine precision is used instead. 00078 * For example, if diag(S)*X=B were the least squares problem, 00079 * where diag(S) is a diagonal matrix of singular values, the 00080 * solution would be X(i) = B(i) / S(i) if S(i) is greater than 00081 * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to 00082 * RCOND*max(S). 00083 * 00084 * RANK (output) INTEGER 00085 * The number of singular values of A greater than RCOND times 00086 * the largest singular value. 00087 * 00088 * WORK (workspace) COMPLEX*16 array, dimension at least 00089 * (N * NRHS). 00090 * 00091 * RWORK (workspace) DOUBLE PRECISION array, dimension at least 00092 * (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00093 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ), 00094 * where 00095 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 00096 * 00097 * IWORK (workspace) INTEGER array, dimension at least 00098 * (3*N*NLVL + 11*N). 00099 * 00100 * INFO (output) INTEGER 00101 * = 0: successful exit. 00102 * < 0: if INFO = -i, the i-th argument had an illegal value. 00103 * > 0: The algorithm failed to compute a singular value while 00104 * working on the submatrix lying in rows and columns 00105 * INFO/(N+1) through MOD(INFO,N+1). 00106 * 00107 * Further Details 00108 * =============== 00109 * 00110 * Based on contributions by 00111 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00112 * California at Berkeley, USA 00113 * Osni Marques, LBNL/NERSC, USA 00114 * 00115 * ===================================================================== 00116 * 00117 * .. Parameters .. 00118 DOUBLE PRECISION ZERO, ONE, TWO 00119 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00120 COMPLEX*16 CZERO 00121 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) ) 00122 * .. 00123 * .. Local Scalars .. 00124 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 00125 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB, 00126 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG, 00127 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB, 00128 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1, 00129 $ U, VT, Z 00130 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL 00131 * .. 00132 * .. External Functions .. 00133 INTEGER IDAMAX 00134 DOUBLE PRECISION DLAMCH, DLANST 00135 EXTERNAL IDAMAX, DLAMCH, DLANST 00136 * .. 00137 * .. External Subroutines .. 00138 EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET, 00139 $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA, 00140 $ ZLASCL, ZLASET 00141 * .. 00142 * .. Intrinsic Functions .. 00143 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN 00144 * .. 00145 * .. Executable Statements .. 00146 * 00147 * Test the input parameters. 00148 * 00149 INFO = 0 00150 * 00151 IF( N.LT.0 ) THEN 00152 INFO = -3 00153 ELSE IF( NRHS.LT.1 ) THEN 00154 INFO = -4 00155 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 00156 INFO = -8 00157 END IF 00158 IF( INFO.NE.0 ) THEN 00159 CALL XERBLA( 'ZLALSD', -INFO ) 00160 RETURN 00161 END IF 00162 * 00163 EPS = DLAMCH( 'Epsilon' ) 00164 * 00165 * Set up the tolerance. 00166 * 00167 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 00168 RCND = EPS 00169 ELSE 00170 RCND = RCOND 00171 END IF 00172 * 00173 RANK = 0 00174 * 00175 * Quick return if possible. 00176 * 00177 IF( N.EQ.0 ) THEN 00178 RETURN 00179 ELSE IF( N.EQ.1 ) THEN 00180 IF( D( 1 ).EQ.ZERO ) THEN 00181 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB ) 00182 ELSE 00183 RANK = 1 00184 CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 00185 D( 1 ) = ABS( D( 1 ) ) 00186 END IF 00187 RETURN 00188 END IF 00189 * 00190 * Rotate the matrix if it is lower bidiagonal. 00191 * 00192 IF( UPLO.EQ.'L' ) THEN 00193 DO 10 I = 1, N - 1 00194 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00195 D( I ) = R 00196 E( I ) = SN*D( I+1 ) 00197 D( I+1 ) = CS*D( I+1 ) 00198 IF( NRHS.EQ.1 ) THEN 00199 CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 00200 ELSE 00201 RWORK( I*2-1 ) = CS 00202 RWORK( I*2 ) = SN 00203 END IF 00204 10 CONTINUE 00205 IF( NRHS.GT.1 ) THEN 00206 DO 30 I = 1, NRHS 00207 DO 20 J = 1, N - 1 00208 CS = RWORK( J*2-1 ) 00209 SN = RWORK( J*2 ) 00210 CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 00211 20 CONTINUE 00212 30 CONTINUE 00213 END IF 00214 END IF 00215 * 00216 * Scale. 00217 * 00218 NM1 = N - 1 00219 ORGNRM = DLANST( 'M', N, D, E ) 00220 IF( ORGNRM.EQ.ZERO ) THEN 00221 CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB ) 00222 RETURN 00223 END IF 00224 * 00225 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00226 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 00227 * 00228 * If N is smaller than the minimum divide size SMLSIZ, then solve 00229 * the problem with another solver. 00230 * 00231 IF( N.LE.SMLSIZ ) THEN 00232 IRWU = 1 00233 IRWVT = IRWU + N*N 00234 IRWWRK = IRWVT + N*N 00235 IRWRB = IRWWRK 00236 IRWIB = IRWRB + N*NRHS 00237 IRWB = IRWIB + N*NRHS 00238 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N ) 00239 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N ) 00240 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N, 00241 $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1, 00242 $ RWORK( IRWWRK ), INFO ) 00243 IF( INFO.NE.0 ) THEN 00244 RETURN 00245 END IF 00246 * 00247 * In the real version, B is passed to DLASDQ and multiplied 00248 * internally by Q**H. Here B is complex and that product is 00249 * computed below in two steps (real and imaginary parts). 00250 * 00251 J = IRWB - 1 00252 DO 50 JCOL = 1, NRHS 00253 DO 40 JROW = 1, N 00254 J = J + 1 00255 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00256 40 CONTINUE 00257 50 CONTINUE 00258 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 00259 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 00260 J = IRWB - 1 00261 DO 70 JCOL = 1, NRHS 00262 DO 60 JROW = 1, N 00263 J = J + 1 00264 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00265 60 CONTINUE 00266 70 CONTINUE 00267 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 00268 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 00269 JREAL = IRWRB - 1 00270 JIMAG = IRWIB - 1 00271 DO 90 JCOL = 1, NRHS 00272 DO 80 JROW = 1, N 00273 JREAL = JREAL + 1 00274 JIMAG = JIMAG + 1 00275 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00276 $ RWORK( JIMAG ) ) 00277 80 CONTINUE 00278 90 CONTINUE 00279 * 00280 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00281 DO 100 I = 1, N 00282 IF( D( I ).LE.TOL ) THEN 00283 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 00284 ELSE 00285 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 00286 $ LDB, INFO ) 00287 RANK = RANK + 1 00288 END IF 00289 100 CONTINUE 00290 * 00291 * Since B is complex, the following call to DGEMM is performed 00292 * in two steps (real and imaginary parts). That is for V * B 00293 * (in the real version of the code V**H is stored in WORK). 00294 * 00295 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 00296 * $ WORK( NWORK ), N ) 00297 * 00298 J = IRWB - 1 00299 DO 120 JCOL = 1, NRHS 00300 DO 110 JROW = 1, N 00301 J = J + 1 00302 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00303 110 CONTINUE 00304 120 CONTINUE 00305 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 00306 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 00307 J = IRWB - 1 00308 DO 140 JCOL = 1, NRHS 00309 DO 130 JROW = 1, N 00310 J = J + 1 00311 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00312 130 CONTINUE 00313 140 CONTINUE 00314 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 00315 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 00316 JREAL = IRWRB - 1 00317 JIMAG = IRWIB - 1 00318 DO 160 JCOL = 1, NRHS 00319 DO 150 JROW = 1, N 00320 JREAL = JREAL + 1 00321 JIMAG = JIMAG + 1 00322 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00323 $ RWORK( JIMAG ) ) 00324 150 CONTINUE 00325 160 CONTINUE 00326 * 00327 * Unscale. 00328 * 00329 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00330 CALL DLASRT( 'D', N, D, INFO ) 00331 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00332 * 00333 RETURN 00334 END IF 00335 * 00336 * Book-keeping and setting up some constants. 00337 * 00338 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 00339 * 00340 SMLSZP = SMLSIZ + 1 00341 * 00342 U = 1 00343 VT = 1 + SMLSIZ*N 00344 DIFL = VT + SMLSZP*N 00345 DIFR = DIFL + NLVL*N 00346 Z = DIFR + NLVL*N*2 00347 C = Z + NLVL*N 00348 S = C + N 00349 POLES = S + N 00350 GIVNUM = POLES + 2*NLVL*N 00351 NRWORK = GIVNUM + 2*NLVL*N 00352 BX = 1 00353 * 00354 IRWRB = NRWORK 00355 IRWIB = IRWRB + SMLSIZ*NRHS 00356 IRWB = IRWIB + SMLSIZ*NRHS 00357 * 00358 SIZEI = 1 + N 00359 K = SIZEI + N 00360 GIVPTR = K + N 00361 PERM = GIVPTR + N 00362 GIVCOL = PERM + NLVL*N 00363 IWK = GIVCOL + NLVL*N*2 00364 * 00365 ST = 1 00366 SQRE = 0 00367 ICMPQ1 = 1 00368 ICMPQ2 = 0 00369 NSUB = 0 00370 * 00371 DO 170 I = 1, N 00372 IF( ABS( D( I ) ).LT.EPS ) THEN 00373 D( I ) = SIGN( EPS, D( I ) ) 00374 END IF 00375 170 CONTINUE 00376 * 00377 DO 240 I = 1, NM1 00378 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 00379 NSUB = NSUB + 1 00380 IWORK( NSUB ) = ST 00381 * 00382 * Subproblem found. First determine its size and then 00383 * apply divide and conquer on it. 00384 * 00385 IF( I.LT.NM1 ) THEN 00386 * 00387 * A subproblem with E(I) small for I < NM1. 00388 * 00389 NSIZE = I - ST + 1 00390 IWORK( SIZEI+NSUB-1 ) = NSIZE 00391 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 00392 * 00393 * A subproblem with E(NM1) not too small but I = NM1. 00394 * 00395 NSIZE = N - ST + 1 00396 IWORK( SIZEI+NSUB-1 ) = NSIZE 00397 ELSE 00398 * 00399 * A subproblem with E(NM1) small. This implies an 00400 * 1-by-1 subproblem at D(N), which is not solved 00401 * explicitly. 00402 * 00403 NSIZE = I - ST + 1 00404 IWORK( SIZEI+NSUB-1 ) = NSIZE 00405 NSUB = NSUB + 1 00406 IWORK( NSUB ) = N 00407 IWORK( SIZEI+NSUB-1 ) = 1 00408 CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 00409 END IF 00410 ST1 = ST - 1 00411 IF( NSIZE.EQ.1 ) THEN 00412 * 00413 * This is a 1-by-1 subproblem and is not solved 00414 * explicitly. 00415 * 00416 CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 00417 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00418 * 00419 * This is a small subproblem and is solved by DLASDQ. 00420 * 00421 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 00422 $ RWORK( VT+ST1 ), N ) 00423 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 00424 $ RWORK( U+ST1 ), N ) 00425 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ), 00426 $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ), 00427 $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ), 00428 $ INFO ) 00429 IF( INFO.NE.0 ) THEN 00430 RETURN 00431 END IF 00432 * 00433 * In the real version, B is passed to DLASDQ and multiplied 00434 * internally by Q**H. Here B is complex and that product is 00435 * computed below in two steps (real and imaginary parts). 00436 * 00437 J = IRWB - 1 00438 DO 190 JCOL = 1, NRHS 00439 DO 180 JROW = ST, ST + NSIZE - 1 00440 J = J + 1 00441 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00442 180 CONTINUE 00443 190 CONTINUE 00444 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00445 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 00446 $ ZERO, RWORK( IRWRB ), NSIZE ) 00447 J = IRWB - 1 00448 DO 210 JCOL = 1, NRHS 00449 DO 200 JROW = ST, ST + NSIZE - 1 00450 J = J + 1 00451 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00452 200 CONTINUE 00453 210 CONTINUE 00454 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00455 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 00456 $ ZERO, RWORK( IRWIB ), NSIZE ) 00457 JREAL = IRWRB - 1 00458 JIMAG = IRWIB - 1 00459 DO 230 JCOL = 1, NRHS 00460 DO 220 JROW = ST, ST + NSIZE - 1 00461 JREAL = JREAL + 1 00462 JIMAG = JIMAG + 1 00463 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00464 $ RWORK( JIMAG ) ) 00465 220 CONTINUE 00466 230 CONTINUE 00467 * 00468 CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 00469 $ WORK( BX+ST1 ), N ) 00470 ELSE 00471 * 00472 * A large problem. Solve it using divide and conquer. 00473 * 00474 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 00475 $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ), 00476 $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ), 00477 $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ), 00478 $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 00479 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 00480 $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ), 00481 $ RWORK( S+ST1 ), RWORK( NRWORK ), 00482 $ IWORK( IWK ), INFO ) 00483 IF( INFO.NE.0 ) THEN 00484 RETURN 00485 END IF 00486 BXST = BX + ST1 00487 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 00488 $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N, 00489 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 00490 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 00491 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 00492 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00493 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 00494 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 00495 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 00496 IF( INFO.NE.0 ) THEN 00497 RETURN 00498 END IF 00499 END IF 00500 ST = I + 1 00501 END IF 00502 240 CONTINUE 00503 * 00504 * Apply the singular values and treat the tiny ones as zero. 00505 * 00506 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00507 * 00508 DO 250 I = 1, N 00509 * 00510 * Some of the elements in D can be negative because 1-by-1 00511 * subproblems were not solved explicitly. 00512 * 00513 IF( ABS( D( I ) ).LE.TOL ) THEN 00514 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N ) 00515 ELSE 00516 RANK = RANK + 1 00517 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 00518 $ WORK( BX+I-1 ), N, INFO ) 00519 END IF 00520 D( I ) = ABS( D( I ) ) 00521 250 CONTINUE 00522 * 00523 * Now apply back the right singular vectors. 00524 * 00525 ICMPQ2 = 1 00526 DO 320 I = 1, NSUB 00527 ST = IWORK( I ) 00528 ST1 = ST - 1 00529 NSIZE = IWORK( SIZEI+I-1 ) 00530 BXST = BX + ST1 00531 IF( NSIZE.EQ.1 ) THEN 00532 CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 00533 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00534 * 00535 * Since B and BX are complex, the following call to DGEMM 00536 * is performed in two steps (real and imaginary parts). 00537 * 00538 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00539 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO, 00540 * $ B( ST, 1 ), LDB ) 00541 * 00542 J = BXST - N - 1 00543 JREAL = IRWB - 1 00544 DO 270 JCOL = 1, NRHS 00545 J = J + N 00546 DO 260 JROW = 1, NSIZE 00547 JREAL = JREAL + 1 00548 RWORK( JREAL ) = DBLE( WORK( J+JROW ) ) 00549 260 CONTINUE 00550 270 CONTINUE 00551 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00552 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 00553 $ RWORK( IRWRB ), NSIZE ) 00554 J = BXST - N - 1 00555 JIMAG = IRWB - 1 00556 DO 290 JCOL = 1, NRHS 00557 J = J + N 00558 DO 280 JROW = 1, NSIZE 00559 JIMAG = JIMAG + 1 00560 RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) ) 00561 280 CONTINUE 00562 290 CONTINUE 00563 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00564 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 00565 $ RWORK( IRWIB ), NSIZE ) 00566 JREAL = IRWRB - 1 00567 JIMAG = IRWIB - 1 00568 DO 310 JCOL = 1, NRHS 00569 DO 300 JROW = ST, ST + NSIZE - 1 00570 JREAL = JREAL + 1 00571 JIMAG = JIMAG + 1 00572 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00573 $ RWORK( JIMAG ) ) 00574 300 CONTINUE 00575 310 CONTINUE 00576 ELSE 00577 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 00578 $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N, 00579 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 00580 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 00581 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 00582 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00583 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 00584 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 00585 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 00586 IF( INFO.NE.0 ) THEN 00587 RETURN 00588 END IF 00589 END IF 00590 320 CONTINUE 00591 * 00592 * Unscale and sort the singular values. 00593 * 00594 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00595 CALL DLASRT( 'D', N, D, INFO ) 00596 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00597 * 00598 RETURN 00599 * 00600 * End of ZLALSD 00601 * 00602 END
1.7.3 
