 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00002 $ WORK, LWORK, RWORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK 00011 DOUBLE PRECISION RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION RWORK( * ), S( * ) 00016 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZGELSD computes the minimum-norm solution to a real linear least 00023 * squares problem: 00024 * minimize 2-norm(| b - A*x |) 00025 * using the singular value decomposition (SVD) of A. A is an M-by-N 00026 * matrix which may be rank-deficient. 00027 * 00028 * Several right hand side vectors b and solution vectors x can be 00029 * handled in a single call; they are stored as the columns of the 00030 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 00031 * matrix X. 00032 * 00033 * The problem is solved in three steps: 00034 * (1) Reduce the coefficient matrix A to bidiagonal form with 00035 * Householder tranformations, reducing the original problem 00036 * into a "bidiagonal least squares problem" (BLS) 00037 * (2) Solve the BLS using a divide and conquer approach. 00038 * (3) Apply back all the Householder tranformations to solve 00039 * the original least squares problem. 00040 * 00041 * The effective rank of A is determined by treating as zero those 00042 * singular values which are less than RCOND times the largest singular 00043 * value. 00044 * 00045 * The divide and conquer algorithm makes very mild assumptions about 00046 * floating point arithmetic. It will work on machines with a guard 00047 * digit in add/subtract, or on those binary machines without guard 00048 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00049 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00050 * without guard digits, but we know of none. 00051 * 00052 * Arguments 00053 * ========= 00054 * 00055 * M (input) INTEGER 00056 * The number of rows of the matrix A. M >= 0. 00057 * 00058 * N (input) INTEGER 00059 * The number of columns of the matrix A. N >= 0. 00060 * 00061 * NRHS (input) INTEGER 00062 * The number of right hand sides, i.e., the number of columns 00063 * of the matrices B and X. NRHS >= 0. 00064 * 00065 * A (input) COMPLEX*16 array, dimension (LDA,N) 00066 * On entry, the M-by-N matrix A. 00067 * On exit, A has been destroyed. 00068 * 00069 * LDA (input) INTEGER 00070 * The leading dimension of the array A. LDA >= max(1,M). 00071 * 00072 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 00073 * On entry, the M-by-NRHS right hand side matrix B. 00074 * On exit, B is overwritten by the N-by-NRHS solution matrix X. 00075 * If m >= n and RANK = n, the residual sum-of-squares for 00076 * the solution in the i-th column is given by the sum of 00077 * squares of the modulus of elements n+1:m in that column. 00078 * 00079 * LDB (input) INTEGER 00080 * The leading dimension of the array B. LDB >= max(1,M,N). 00081 * 00082 * S (output) DOUBLE PRECISION array, dimension (min(M,N)) 00083 * The singular values of A in decreasing order. 00084 * The condition number of A in the 2-norm = S(1)/S(min(m,n)). 00085 * 00086 * RCOND (input) DOUBLE PRECISION 00087 * RCOND is used to determine the effective rank of A. 00088 * Singular values S(i) <= RCOND*S(1) are treated as zero. 00089 * If RCOND < 0, machine precision is used instead. 00090 * 00091 * RANK (output) INTEGER 00092 * The effective rank of A, i.e., the number of singular values 00093 * which are greater than RCOND*S(1). 00094 * 00095 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00096 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00097 * 00098 * LWORK (input) INTEGER 00099 * The dimension of the array WORK. LWORK must be at least 1. 00100 * The exact minimum amount of workspace needed depends on M, 00101 * N and NRHS. As long as LWORK is at least 00102 * 2*N + N*NRHS 00103 * if M is greater than or equal to N or 00104 * 2*M + M*NRHS 00105 * if M is less than N, the code will execute correctly. 00106 * For good performance, LWORK should generally be larger. 00107 * 00108 * If LWORK = -1, then a workspace query is assumed; the routine 00109 * only calculates the optimal size of the array WORK and the 00110 * minimum sizes of the arrays RWORK and IWORK, and returns 00111 * these values as the first entries of the WORK, RWORK and 00112 * IWORK arrays, and no error message related to LWORK is issued 00113 * by XERBLA. 00114 * 00115 * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 00116 * LRWORK >= 00117 * 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00118 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00119 * if M is greater than or equal to N or 00120 * 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + 00121 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00122 * if M is less than N, the code will execute correctly. 00123 * SMLSIZ is returned by ILAENV and is equal to the maximum 00124 * size of the subproblems at the bottom of the computation 00125 * tree (usually about 25), and 00126 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 00127 * On exit, if INFO = 0, RWORK(1) returns the minimum LRWORK. 00128 * 00129 * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) 00130 * LIWORK >= max(1, 3*MINMN*NLVL + 11*MINMN), 00131 * where MINMN = MIN( M,N ). 00132 * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK. 00133 * 00134 * INFO (output) INTEGER 00135 * = 0: successful exit 00136 * < 0: if INFO = -i, the i-th argument had an illegal value. 00137 * > 0: the algorithm for computing the SVD failed to converge; 00138 * if INFO = i, i off-diagonal elements of an intermediate 00139 * bidiagonal form did not converge to zero. 00140 * 00141 * Further Details 00142 * =============== 00143 * 00144 * Based on contributions by 00145 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00146 * California at Berkeley, USA 00147 * Osni Marques, LBNL/NERSC, USA 00148 * 00149 * ===================================================================== 00150 * 00151 * .. Parameters .. 00152 DOUBLE PRECISION ZERO, ONE, TWO 00153 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00154 COMPLEX*16 CZERO 00155 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) ) 00156 * .. 00157 * .. Local Scalars .. 00158 LOGICAL LQUERY 00159 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ, 00160 $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN, 00161 $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ 00162 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF, 00166 $ ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR, 00167 $ ZUNMLQ, ZUNMQR 00168 * .. 00169 * .. External Functions .. 00170 INTEGER ILAENV 00171 DOUBLE PRECISION DLAMCH, ZLANGE 00172 EXTERNAL ILAENV, DLAMCH, ZLANGE 00173 * .. 00174 * .. Intrinsic Functions .. 00175 INTRINSIC INT, LOG, MAX, MIN, DBLE 00176 * .. 00177 * .. Executable Statements .. 00178 * 00179 * Test the input arguments. 00180 * 00181 INFO = 0 00182 MINMN = MIN( M, N ) 00183 MAXMN = MAX( M, N ) 00184 LQUERY = ( LWORK.EQ.-1 ) 00185 IF( M.LT.0 ) THEN 00186 INFO = -1 00187 ELSE IF( N.LT.0 ) THEN 00188 INFO = -2 00189 ELSE IF( NRHS.LT.0 ) THEN 00190 INFO = -3 00191 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00192 INFO = -5 00193 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00194 INFO = -7 00195 END IF 00196 * 00197 * Compute workspace. 00198 * (Note: Comments in the code beginning "Workspace:" describe the 00199 * minimal amount of workspace needed at that point in the code, 00200 * as well as the preferred amount for good performance. 00201 * NB refers to the optimal block size for the immediately 00202 * following subroutine, as returned by ILAENV.) 00203 * 00204 IF( INFO.EQ.0 ) THEN 00205 MINWRK = 1 00206 MAXWRK = 1 00207 LIWORK = 1 00208 LRWORK = 1 00209 IF( MINMN.GT.0 ) THEN 00210 SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 ) 00211 MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 ) 00212 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) / 00213 $ LOG( TWO ) ) + 1, 0 ) 00214 LIWORK = 3*MINMN*NLVL + 11*MINMN 00215 MM = M 00216 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00217 * 00218 * Path 1a - overdetermined, with many more rows than 00219 * columns. 00220 * 00221 MM = N 00222 MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N, 00223 $ -1, -1 ) ) 00224 MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M, 00225 $ NRHS, N, -1 ) ) 00226 END IF 00227 IF( M.GE.N ) THEN 00228 * 00229 * Path 1 - overdetermined or exactly determined. 00230 * 00231 LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00232 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00233 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, 00234 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) 00235 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', 00236 $ 'QLC', MM, NRHS, N, -1 ) ) 00237 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00238 $ 'ZUNMBR', 'PLN', N, NRHS, N, -1 ) ) 00239 MAXWRK = MAX( MAXWRK, 2*N + N*NRHS ) 00240 MINWRK = MAX( 2*N + MM, 2*N + N*NRHS ) 00241 END IF 00242 IF( N.GT.M ) THEN 00243 LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + 00244 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00245 IF( N.GE.MNTHR ) THEN 00246 * 00247 * Path 2a - underdetermined, with many more columns 00248 * than rows. 00249 * 00250 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, 00251 $ -1 ) 00252 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, 00253 $ 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00254 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, 00255 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) ) 00256 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1, 00257 $ 'ZUNMLQ', 'LC', N, NRHS, M, -1 ) ) 00258 IF( NRHS.GT.1 ) THEN 00259 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00260 ELSE 00261 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00262 END IF 00263 MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS ) 00264 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00265 ! calculation should use queries for all routines eventually. 00266 MAXWRK = MAX( MAXWRK, 00267 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00268 ELSE 00269 * 00270 * Path 2 - underdetermined. 00271 * 00272 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M, 00273 $ N, -1, -1 ) 00274 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR', 00275 $ 'QLC', M, NRHS, M, -1 ) ) 00276 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR', 00277 $ 'PLN', N, NRHS, M, -1 ) ) 00278 MAXWRK = MAX( MAXWRK, 2*M + M*NRHS ) 00279 END IF 00280 MINWRK = MAX( 2*M + N, 2*M + M*NRHS ) 00281 END IF 00282 END IF 00283 MINWRK = MIN( MINWRK, MAXWRK ) 00284 WORK( 1 ) = MAXWRK 00285 IWORK( 1 ) = LIWORK 00286 RWORK( 1 ) = LRWORK 00287 * 00288 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00289 INFO = -12 00290 END IF 00291 END IF 00292 * 00293 IF( INFO.NE.0 ) THEN 00294 CALL XERBLA( 'ZGELSD', -INFO ) 00295 RETURN 00296 ELSE IF( LQUERY ) THEN 00297 RETURN 00298 END IF 00299 * 00300 * Quick return if possible. 00301 * 00302 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00303 RANK = 0 00304 RETURN 00305 END IF 00306 * 00307 * Get machine parameters. 00308 * 00309 EPS = DLAMCH( 'P' ) 00310 SFMIN = DLAMCH( 'S' ) 00311 SMLNUM = SFMIN / EPS 00312 BIGNUM = ONE / SMLNUM 00313 CALL DLABAD( SMLNUM, BIGNUM ) 00314 * 00315 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00316 * 00317 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 00318 IASCL = 0 00319 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00320 * 00321 * Scale matrix norm up to SMLNUM 00322 * 00323 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00324 IASCL = 1 00325 ELSE IF( ANRM.GT.BIGNUM ) THEN 00326 * 00327 * Scale matrix norm down to BIGNUM. 00328 * 00329 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00330 IASCL = 2 00331 ELSE IF( ANRM.EQ.ZERO ) THEN 00332 * 00333 * Matrix all zero. Return zero solution. 00334 * 00335 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00336 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00337 RANK = 0 00338 GO TO 10 00339 END IF 00340 * 00341 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00342 * 00343 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00344 IBSCL = 0 00345 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00346 * 00347 * Scale matrix norm up to SMLNUM. 00348 * 00349 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00350 IBSCL = 1 00351 ELSE IF( BNRM.GT.BIGNUM ) THEN 00352 * 00353 * Scale matrix norm down to BIGNUM. 00354 * 00355 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00356 IBSCL = 2 00357 END IF 00358 * 00359 * If M < N make sure B(M+1:N,:) = 0 00360 * 00361 IF( M.LT.N ) 00362 $ CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00363 * 00364 * Overdetermined case. 00365 * 00366 IF( M.GE.N ) THEN 00367 * 00368 * Path 1 - overdetermined or exactly determined. 00369 * 00370 MM = M 00371 IF( M.GE.MNTHR ) THEN 00372 * 00373 * Path 1a - overdetermined, with many more rows than columns 00374 * 00375 MM = N 00376 ITAU = 1 00377 NWORK = ITAU + N 00378 * 00379 * Compute A=Q*R. 00380 * (RWorkspace: need N) 00381 * (CWorkspace: need N, prefer N*NB) 00382 * 00383 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00384 $ LWORK-NWORK+1, INFO ) 00385 * 00386 * Multiply B by transpose(Q). 00387 * (RWorkspace: need N) 00388 * (CWorkspace: need NRHS, prefer NRHS*NB) 00389 * 00390 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00391 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00392 * 00393 * Zero out below R. 00394 * 00395 IF( N.GT.1 ) THEN 00396 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00397 $ LDA ) 00398 END IF 00399 END IF 00400 * 00401 ITAUQ = 1 00402 ITAUP = ITAUQ + N 00403 NWORK = ITAUP + N 00404 IE = 1 00405 NRWORK = IE + N 00406 * 00407 * Bidiagonalize R in A. 00408 * (RWorkspace: need N) 00409 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB) 00410 * 00411 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00412 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00413 $ INFO ) 00414 * 00415 * Multiply B by transpose of left bidiagonalizing vectors of R. 00416 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB) 00417 * 00418 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00419 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00420 * 00421 * Solve the bidiagonal least squares problem. 00422 * 00423 CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB, 00424 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00425 $ IWORK, INFO ) 00426 IF( INFO.NE.0 ) THEN 00427 GO TO 10 00428 END IF 00429 * 00430 * Multiply B by right bidiagonalizing vectors of R. 00431 * 00432 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00433 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00434 * 00435 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00436 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00437 * 00438 * Path 2a - underdetermined, with many more columns than rows 00439 * and sufficient workspace for an efficient algorithm. 00440 * 00441 LDWORK = M 00442 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00443 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00444 ITAU = 1 00445 NWORK = M + 1 00446 * 00447 * Compute A=L*Q. 00448 * (CWorkspace: need 2*M, prefer M+M*NB) 00449 * 00450 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00451 $ LWORK-NWORK+1, INFO ) 00452 IL = NWORK 00453 * 00454 * Copy L to WORK(IL), zeroing out above its diagonal. 00455 * 00456 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00457 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ), 00458 $ LDWORK ) 00459 ITAUQ = IL + LDWORK*M 00460 ITAUP = ITAUQ + M 00461 NWORK = ITAUP + M 00462 IE = 1 00463 NRWORK = IE + M 00464 * 00465 * Bidiagonalize L in WORK(IL). 00466 * (RWorkspace: need M) 00467 * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB) 00468 * 00469 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ), 00470 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00471 $ LWORK-NWORK+1, INFO ) 00472 * 00473 * Multiply B by transpose of left bidiagonalizing vectors of L. 00474 * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00475 * 00476 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK, 00477 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00478 $ LWORK-NWORK+1, INFO ) 00479 * 00480 * Solve the bidiagonal least squares problem. 00481 * 00482 CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00483 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00484 $ IWORK, INFO ) 00485 IF( INFO.NE.0 ) THEN 00486 GO TO 10 00487 END IF 00488 * 00489 * Multiply B by right bidiagonalizing vectors of L. 00490 * 00491 CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00492 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00493 $ LWORK-NWORK+1, INFO ) 00494 * 00495 * Zero out below first M rows of B. 00496 * 00497 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00498 NWORK = ITAU + M 00499 * 00500 * Multiply transpose(Q) by B. 00501 * (CWorkspace: need NRHS, prefer NRHS*NB) 00502 * 00503 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00504 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00505 * 00506 ELSE 00507 * 00508 * Path 2 - remaining underdetermined cases. 00509 * 00510 ITAUQ = 1 00511 ITAUP = ITAUQ + M 00512 NWORK = ITAUP + M 00513 IE = 1 00514 NRWORK = IE + M 00515 * 00516 * Bidiagonalize A. 00517 * (RWorkspace: need M) 00518 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 00519 * 00520 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00521 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00522 $ INFO ) 00523 * 00524 * Multiply B by transpose of left bidiagonalizing vectors. 00525 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB) 00526 * 00527 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00528 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00529 * 00530 * Solve the bidiagonal least squares problem. 00531 * 00532 CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00533 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00534 $ IWORK, INFO ) 00535 IF( INFO.NE.0 ) THEN 00536 GO TO 10 00537 END IF 00538 * 00539 * Multiply B by right bidiagonalizing vectors of A. 00540 * 00541 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00542 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00543 * 00544 END IF 00545 * 00546 * Undo scaling. 00547 * 00548 IF( IASCL.EQ.1 ) THEN 00549 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00550 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00551 $ INFO ) 00552 ELSE IF( IASCL.EQ.2 ) THEN 00553 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00554 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00555 $ INFO ) 00556 END IF 00557 IF( IBSCL.EQ.1 ) THEN 00558 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00559 ELSE IF( IBSCL.EQ.2 ) THEN 00560 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00561 END IF 00562 * 00563 10 CONTINUE 00564 WORK( 1 ) = MAXWRK 00565 IWORK( 1 ) = LIWORK 00566 RWORK( 1 ) = LRWORK 00567 RETURN 00568 * 00569 * End of ZGELSD 00570 * 00571 END
1.7.3 
