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