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