![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLALSD 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLALSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLALSD( 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 * DOUBLE PRECISION RCOND 00028 * .. 00029 * .. Array Arguments .. 00030 * INTEGER IWORK( * ) 00031 * DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 00032 * COMPLEX*16 B( LDB, * ), WORK( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> ZLALSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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*16 array, dimension at least 00140 *> (N * NRHS). 00141 *> \endverbatim 00142 *> 00143 *> \param[out] RWORK 00144 *> \verbatim 00145 *> RWORK is DOUBLE PRECISION array, dimension at least 00146 *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00147 *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ), 00148 *> where 00149 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 00150 *> \endverbatim 00151 *> 00152 *> \param[out] IWORK 00153 *> \verbatim 00154 *> IWORK is INTEGER array, dimension at least 00155 *> (3*N*NLVL + 11*N). 00156 *> \endverbatim 00157 *> 00158 *> \param[out] INFO 00159 *> \verbatim 00160 *> INFO is INTEGER 00161 *> = 0: successful exit. 00162 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00163 *> > 0: The algorithm failed to compute a singular value while 00164 *> working on the submatrix lying in rows and columns 00165 *> INFO/(N+1) through MOD(INFO,N+1). 00166 *> \endverbatim 00167 * 00168 * Authors: 00169 * ======== 00170 * 00171 *> \author Univ. of Tennessee 00172 *> \author Univ. of California Berkeley 00173 *> \author Univ. of Colorado Denver 00174 *> \author NAG Ltd. 00175 * 00176 *> \date November 2011 00177 * 00178 *> \ingroup complex16OTHERcomputational 00179 * 00180 *> \par Contributors: 00181 * ================== 00182 *> 00183 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of 00184 *> California at Berkeley, USA \n 00185 *> Osni Marques, LBNL/NERSC, USA \n 00186 * 00187 * ===================================================================== 00188 SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, 00189 $ RANK, WORK, RWORK, IWORK, INFO ) 00190 * 00191 * -- LAPACK computational routine (version 3.4.0) -- 00192 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00193 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00194 * November 2011 00195 * 00196 * .. Scalar Arguments .. 00197 CHARACTER UPLO 00198 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ 00199 DOUBLE PRECISION RCOND 00200 * .. 00201 * .. Array Arguments .. 00202 INTEGER IWORK( * ) 00203 DOUBLE PRECISION D( * ), E( * ), RWORK( * ) 00204 COMPLEX*16 B( LDB, * ), WORK( * ) 00205 * .. 00206 * 00207 * ===================================================================== 00208 * 00209 * .. Parameters .. 00210 DOUBLE PRECISION ZERO, ONE, TWO 00211 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 00212 COMPLEX*16 CZERO 00213 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) ) 00214 * .. 00215 * .. Local Scalars .. 00216 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 00217 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB, 00218 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG, 00219 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB, 00220 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1, 00221 $ U, VT, Z 00222 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL 00223 * .. 00224 * .. External Functions .. 00225 INTEGER IDAMAX 00226 DOUBLE PRECISION DLAMCH, DLANST 00227 EXTERNAL IDAMAX, DLAMCH, DLANST 00228 * .. 00229 * .. External Subroutines .. 00230 EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET, 00231 $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA, 00232 $ ZLASCL, ZLASET 00233 * .. 00234 * .. Intrinsic Functions .. 00235 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN 00236 * .. 00237 * .. Executable Statements .. 00238 * 00239 * Test the input parameters. 00240 * 00241 INFO = 0 00242 * 00243 IF( N.LT.0 ) THEN 00244 INFO = -3 00245 ELSE IF( NRHS.LT.1 ) THEN 00246 INFO = -4 00247 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 00248 INFO = -8 00249 END IF 00250 IF( INFO.NE.0 ) THEN 00251 CALL XERBLA( 'ZLALSD', -INFO ) 00252 RETURN 00253 END IF 00254 * 00255 EPS = DLAMCH( 'Epsilon' ) 00256 * 00257 * Set up the tolerance. 00258 * 00259 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 00260 RCND = EPS 00261 ELSE 00262 RCND = RCOND 00263 END IF 00264 * 00265 RANK = 0 00266 * 00267 * Quick return if possible. 00268 * 00269 IF( N.EQ.0 ) THEN 00270 RETURN 00271 ELSE IF( N.EQ.1 ) THEN 00272 IF( D( 1 ).EQ.ZERO ) THEN 00273 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB ) 00274 ELSE 00275 RANK = 1 00276 CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 00277 D( 1 ) = ABS( D( 1 ) ) 00278 END IF 00279 RETURN 00280 END IF 00281 * 00282 * Rotate the matrix if it is lower bidiagonal. 00283 * 00284 IF( UPLO.EQ.'L' ) THEN 00285 DO 10 I = 1, N - 1 00286 CALL DLARTG( D( I ), E( I ), CS, SN, R ) 00287 D( I ) = R 00288 E( I ) = SN*D( I+1 ) 00289 D( I+1 ) = CS*D( I+1 ) 00290 IF( NRHS.EQ.1 ) THEN 00291 CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 00292 ELSE 00293 RWORK( I*2-1 ) = CS 00294 RWORK( I*2 ) = SN 00295 END IF 00296 10 CONTINUE 00297 IF( NRHS.GT.1 ) THEN 00298 DO 30 I = 1, NRHS 00299 DO 20 J = 1, N - 1 00300 CS = RWORK( J*2-1 ) 00301 SN = RWORK( J*2 ) 00302 CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 00303 20 CONTINUE 00304 30 CONTINUE 00305 END IF 00306 END IF 00307 * 00308 * Scale. 00309 * 00310 NM1 = N - 1 00311 ORGNRM = DLANST( 'M', N, D, E ) 00312 IF( ORGNRM.EQ.ZERO ) THEN 00313 CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB ) 00314 RETURN 00315 END IF 00316 * 00317 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 00318 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 00319 * 00320 * If N is smaller than the minimum divide size SMLSIZ, then solve 00321 * the problem with another solver. 00322 * 00323 IF( N.LE.SMLSIZ ) THEN 00324 IRWU = 1 00325 IRWVT = IRWU + N*N 00326 IRWWRK = IRWVT + N*N 00327 IRWRB = IRWWRK 00328 IRWIB = IRWRB + N*NRHS 00329 IRWB = IRWIB + N*NRHS 00330 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N ) 00331 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N ) 00332 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N, 00333 $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1, 00334 $ RWORK( IRWWRK ), INFO ) 00335 IF( INFO.NE.0 ) THEN 00336 RETURN 00337 END IF 00338 * 00339 * In the real version, B is passed to DLASDQ and multiplied 00340 * internally by Q**H. Here B is complex and that product is 00341 * computed below in two steps (real and imaginary parts). 00342 * 00343 J = IRWB - 1 00344 DO 50 JCOL = 1, NRHS 00345 DO 40 JROW = 1, N 00346 J = J + 1 00347 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00348 40 CONTINUE 00349 50 CONTINUE 00350 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 00351 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 00352 J = IRWB - 1 00353 DO 70 JCOL = 1, NRHS 00354 DO 60 JROW = 1, N 00355 J = J + 1 00356 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00357 60 CONTINUE 00358 70 CONTINUE 00359 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 00360 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 00361 JREAL = IRWRB - 1 00362 JIMAG = IRWIB - 1 00363 DO 90 JCOL = 1, NRHS 00364 DO 80 JROW = 1, N 00365 JREAL = JREAL + 1 00366 JIMAG = JIMAG + 1 00367 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00368 $ RWORK( JIMAG ) ) 00369 80 CONTINUE 00370 90 CONTINUE 00371 * 00372 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00373 DO 100 I = 1, N 00374 IF( D( I ).LE.TOL ) THEN 00375 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 00376 ELSE 00377 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 00378 $ LDB, INFO ) 00379 RANK = RANK + 1 00380 END IF 00381 100 CONTINUE 00382 * 00383 * Since B is complex, the following call to DGEMM is performed 00384 * in two steps (real and imaginary parts). That is for V * B 00385 * (in the real version of the code V**H is stored in WORK). 00386 * 00387 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 00388 * $ WORK( NWORK ), N ) 00389 * 00390 J = IRWB - 1 00391 DO 120 JCOL = 1, NRHS 00392 DO 110 JROW = 1, N 00393 J = J + 1 00394 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00395 110 CONTINUE 00396 120 CONTINUE 00397 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 00398 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 00399 J = IRWB - 1 00400 DO 140 JCOL = 1, NRHS 00401 DO 130 JROW = 1, N 00402 J = J + 1 00403 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00404 130 CONTINUE 00405 140 CONTINUE 00406 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 00407 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 00408 JREAL = IRWRB - 1 00409 JIMAG = IRWIB - 1 00410 DO 160 JCOL = 1, NRHS 00411 DO 150 JROW = 1, N 00412 JREAL = JREAL + 1 00413 JIMAG = JIMAG + 1 00414 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00415 $ RWORK( JIMAG ) ) 00416 150 CONTINUE 00417 160 CONTINUE 00418 * 00419 * Unscale. 00420 * 00421 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00422 CALL DLASRT( 'D', N, D, INFO ) 00423 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00424 * 00425 RETURN 00426 END IF 00427 * 00428 * Book-keeping and setting up some constants. 00429 * 00430 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 00431 * 00432 SMLSZP = SMLSIZ + 1 00433 * 00434 U = 1 00435 VT = 1 + SMLSIZ*N 00436 DIFL = VT + SMLSZP*N 00437 DIFR = DIFL + NLVL*N 00438 Z = DIFR + NLVL*N*2 00439 C = Z + NLVL*N 00440 S = C + N 00441 POLES = S + N 00442 GIVNUM = POLES + 2*NLVL*N 00443 NRWORK = GIVNUM + 2*NLVL*N 00444 BX = 1 00445 * 00446 IRWRB = NRWORK 00447 IRWIB = IRWRB + SMLSIZ*NRHS 00448 IRWB = IRWIB + SMLSIZ*NRHS 00449 * 00450 SIZEI = 1 + N 00451 K = SIZEI + N 00452 GIVPTR = K + N 00453 PERM = GIVPTR + N 00454 GIVCOL = PERM + NLVL*N 00455 IWK = GIVCOL + NLVL*N*2 00456 * 00457 ST = 1 00458 SQRE = 0 00459 ICMPQ1 = 1 00460 ICMPQ2 = 0 00461 NSUB = 0 00462 * 00463 DO 170 I = 1, N 00464 IF( ABS( D( I ) ).LT.EPS ) THEN 00465 D( I ) = SIGN( EPS, D( I ) ) 00466 END IF 00467 170 CONTINUE 00468 * 00469 DO 240 I = 1, NM1 00470 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 00471 NSUB = NSUB + 1 00472 IWORK( NSUB ) = ST 00473 * 00474 * Subproblem found. First determine its size and then 00475 * apply divide and conquer on it. 00476 * 00477 IF( I.LT.NM1 ) THEN 00478 * 00479 * A subproblem with E(I) small for I < NM1. 00480 * 00481 NSIZE = I - ST + 1 00482 IWORK( SIZEI+NSUB-1 ) = NSIZE 00483 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 00484 * 00485 * A subproblem with E(NM1) not too small but I = NM1. 00486 * 00487 NSIZE = N - ST + 1 00488 IWORK( SIZEI+NSUB-1 ) = NSIZE 00489 ELSE 00490 * 00491 * A subproblem with E(NM1) small. This implies an 00492 * 1-by-1 subproblem at D(N), which is not solved 00493 * explicitly. 00494 * 00495 NSIZE = I - ST + 1 00496 IWORK( SIZEI+NSUB-1 ) = NSIZE 00497 NSUB = NSUB + 1 00498 IWORK( NSUB ) = N 00499 IWORK( SIZEI+NSUB-1 ) = 1 00500 CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 00501 END IF 00502 ST1 = ST - 1 00503 IF( NSIZE.EQ.1 ) THEN 00504 * 00505 * This is a 1-by-1 subproblem and is not solved 00506 * explicitly. 00507 * 00508 CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 00509 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00510 * 00511 * This is a small subproblem and is solved by DLASDQ. 00512 * 00513 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 00514 $ RWORK( VT+ST1 ), N ) 00515 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 00516 $ RWORK( U+ST1 ), N ) 00517 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ), 00518 $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ), 00519 $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ), 00520 $ INFO ) 00521 IF( INFO.NE.0 ) THEN 00522 RETURN 00523 END IF 00524 * 00525 * In the real version, B is passed to DLASDQ and multiplied 00526 * internally by Q**H. Here B is complex and that product is 00527 * computed below in two steps (real and imaginary parts). 00528 * 00529 J = IRWB - 1 00530 DO 190 JCOL = 1, NRHS 00531 DO 180 JROW = ST, ST + NSIZE - 1 00532 J = J + 1 00533 RWORK( J ) = DBLE( B( JROW, JCOL ) ) 00534 180 CONTINUE 00535 190 CONTINUE 00536 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00537 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 00538 $ ZERO, RWORK( IRWRB ), NSIZE ) 00539 J = IRWB - 1 00540 DO 210 JCOL = 1, NRHS 00541 DO 200 JROW = ST, ST + NSIZE - 1 00542 J = J + 1 00543 RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 00544 200 CONTINUE 00545 210 CONTINUE 00546 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00547 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 00548 $ ZERO, RWORK( IRWIB ), NSIZE ) 00549 JREAL = IRWRB - 1 00550 JIMAG = IRWIB - 1 00551 DO 230 JCOL = 1, NRHS 00552 DO 220 JROW = ST, ST + NSIZE - 1 00553 JREAL = JREAL + 1 00554 JIMAG = JIMAG + 1 00555 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00556 $ RWORK( JIMAG ) ) 00557 220 CONTINUE 00558 230 CONTINUE 00559 * 00560 CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 00561 $ WORK( BX+ST1 ), N ) 00562 ELSE 00563 * 00564 * A large problem. Solve it using divide and conquer. 00565 * 00566 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 00567 $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ), 00568 $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ), 00569 $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ), 00570 $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 00571 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 00572 $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ), 00573 $ RWORK( S+ST1 ), RWORK( NRWORK ), 00574 $ IWORK( IWK ), INFO ) 00575 IF( INFO.NE.0 ) THEN 00576 RETURN 00577 END IF 00578 BXST = BX + ST1 00579 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 00580 $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N, 00581 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 00582 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 00583 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 00584 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00585 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 00586 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 00587 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 00588 IF( INFO.NE.0 ) THEN 00589 RETURN 00590 END IF 00591 END IF 00592 ST = I + 1 00593 END IF 00594 240 CONTINUE 00595 * 00596 * Apply the singular values and treat the tiny ones as zero. 00597 * 00598 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 00599 * 00600 DO 250 I = 1, N 00601 * 00602 * Some of the elements in D can be negative because 1-by-1 00603 * subproblems were not solved explicitly. 00604 * 00605 IF( ABS( D( I ) ).LE.TOL ) THEN 00606 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N ) 00607 ELSE 00608 RANK = RANK + 1 00609 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 00610 $ WORK( BX+I-1 ), N, INFO ) 00611 END IF 00612 D( I ) = ABS( D( I ) ) 00613 250 CONTINUE 00614 * 00615 * Now apply back the right singular vectors. 00616 * 00617 ICMPQ2 = 1 00618 DO 320 I = 1, NSUB 00619 ST = IWORK( I ) 00620 ST1 = ST - 1 00621 NSIZE = IWORK( SIZEI+I-1 ) 00622 BXST = BX + ST1 00623 IF( NSIZE.EQ.1 ) THEN 00624 CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 00625 ELSE IF( NSIZE.LE.SMLSIZ ) THEN 00626 * 00627 * Since B and BX are complex, the following call to DGEMM 00628 * is performed in two steps (real and imaginary parts). 00629 * 00630 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00631 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO, 00632 * $ B( ST, 1 ), LDB ) 00633 * 00634 J = BXST - N - 1 00635 JREAL = IRWB - 1 00636 DO 270 JCOL = 1, NRHS 00637 J = J + N 00638 DO 260 JROW = 1, NSIZE 00639 JREAL = JREAL + 1 00640 RWORK( JREAL ) = DBLE( WORK( J+JROW ) ) 00641 260 CONTINUE 00642 270 CONTINUE 00643 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00644 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 00645 $ RWORK( IRWRB ), NSIZE ) 00646 J = BXST - N - 1 00647 JIMAG = IRWB - 1 00648 DO 290 JCOL = 1, NRHS 00649 J = J + N 00650 DO 280 JROW = 1, NSIZE 00651 JIMAG = JIMAG + 1 00652 RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) ) 00653 280 CONTINUE 00654 290 CONTINUE 00655 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 00656 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 00657 $ RWORK( IRWIB ), NSIZE ) 00658 JREAL = IRWRB - 1 00659 JIMAG = IRWIB - 1 00660 DO 310 JCOL = 1, NRHS 00661 DO 300 JROW = ST, ST + NSIZE - 1 00662 JREAL = JREAL + 1 00663 JIMAG = JIMAG + 1 00664 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 00665 $ RWORK( JIMAG ) ) 00666 300 CONTINUE 00667 310 CONTINUE 00668 ELSE 00669 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 00670 $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N, 00671 $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 00672 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 00673 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 00674 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 00675 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 00676 $ RWORK( C+ST1 ), RWORK( S+ST1 ), 00677 $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 00678 IF( INFO.NE.0 ) THEN 00679 RETURN 00680 END IF 00681 END IF 00682 320 CONTINUE 00683 * 00684 * Unscale and sort the singular values. 00685 * 00686 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 00687 CALL DLASRT( 'D', N, D, INFO ) 00688 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 00689 * 00690 RETURN 00691 * 00692 * End of ZLALSD 00693 * 00694 END