 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGELSD( 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 REAL RCOND 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL RWORK( * ), S( * ) 00016 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CGELSD 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/output) COMPLEX 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 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) REAL 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) REAL 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 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) REAL 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 REAL ZERO, ONE, TWO 00153 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00154 COMPLEX CZERO 00155 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+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 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL CGEBRD, CGELQF, CGEQRF, CLACPY, 00166 $ CLALSD, CLASCL, CLASET, CUNMBR, 00167 $ CUNMLQ, CUNMQR, SLABAD, SLASCL, 00168 $ SLASET, XERBLA 00169 * .. 00170 * .. External Functions .. 00171 INTEGER ILAENV 00172 REAL CLANGE, SLAMCH 00173 EXTERNAL CLANGE, SLAMCH, ILAENV 00174 * .. 00175 * .. Intrinsic Functions .. 00176 INTRINSIC INT, LOG, MAX, MIN, REAL 00177 * .. 00178 * .. Executable Statements .. 00179 * 00180 * Test the input arguments. 00181 * 00182 INFO = 0 00183 MINMN = MIN( M, N ) 00184 MAXMN = MAX( M, N ) 00185 LQUERY = ( LWORK.EQ.-1 ) 00186 IF( M.LT.0 ) THEN 00187 INFO = -1 00188 ELSE IF( N.LT.0 ) THEN 00189 INFO = -2 00190 ELSE IF( NRHS.LT.0 ) THEN 00191 INFO = -3 00192 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00193 INFO = -5 00194 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00195 INFO = -7 00196 END IF 00197 * 00198 * Compute workspace. 00199 * (Note: Comments in the code beginning "Workspace:" describe the 00200 * minimal amount of workspace needed at that point in the code, 00201 * as well as the preferred amount for good performance. 00202 * NB refers to the optimal block size for the immediately 00203 * following subroutine, as returned by ILAENV.) 00204 * 00205 IF( INFO.EQ.0 ) THEN 00206 MINWRK = 1 00207 MAXWRK = 1 00208 LIWORK = 1 00209 LRWORK = 1 00210 IF( MINMN.GT.0 ) THEN 00211 SMLSIZ = ILAENV( 9, 'CGELSD', ' ', 0, 0, 0, 0 ) 00212 MNTHR = ILAENV( 6, 'CGELSD', ' ', M, N, NRHS, -1 ) 00213 NLVL = MAX( INT( LOG( REAL( MINMN ) / REAL( SMLSIZ + 1 ) ) / 00214 $ LOG( TWO ) ) + 1, 0 ) 00215 LIWORK = 3*MINMN*NLVL + 11*MINMN 00216 MM = M 00217 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00218 * 00219 * Path 1a - overdetermined, with many more rows than 00220 * columns. 00221 * 00222 MM = N 00223 MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'CGEQRF', ' ', M, N, 00224 $ -1, -1 ) ) 00225 MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'CUNMQR', 'LC', M, 00226 $ NRHS, N, -1 ) ) 00227 END IF 00228 IF( M.GE.N ) THEN 00229 * 00230 * Path 1 - overdetermined or exactly determined. 00231 * 00232 LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00233 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00234 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, 00235 $ 'CGEBRD', ' ', MM, N, -1, -1 ) ) 00236 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR', 00237 $ 'QLC', MM, NRHS, N, -1 ) ) 00238 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00239 $ 'CUNMBR', 'PLN', N, NRHS, N, -1 ) ) 00240 MAXWRK = MAX( MAXWRK, 2*N + N*NRHS ) 00241 MINWRK = MAX( 2*N + MM, 2*N + N*NRHS ) 00242 END IF 00243 IF( N.GT.M ) THEN 00244 LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + 00245 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00246 IF( N.GE.MNTHR ) THEN 00247 * 00248 * Path 2a - underdetermined, with many more columns 00249 * than rows. 00250 * 00251 MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, 00252 $ -1 ) 00253 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, 00254 $ 'CGEBRD', ' ', M, M, -1, -1 ) ) 00255 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, 00256 $ 'CUNMBR', 'QLC', M, NRHS, M, -1 ) ) 00257 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1, 00258 $ 'CUNMLQ', 'LC', N, NRHS, M, -1 ) ) 00259 IF( NRHS.GT.1 ) THEN 00260 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00261 ELSE 00262 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00263 END IF 00264 MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS ) 00265 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00266 ! calculation should use queries for all routines eventually. 00267 MAXWRK = MAX( MAXWRK, 00268 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00269 ELSE 00270 * 00271 * Path 2 - underdetermined. 00272 * 00273 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M, 00274 $ N, -1, -1 ) 00275 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR', 00276 $ 'QLC', M, NRHS, M, -1 ) ) 00277 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNMBR', 00278 $ 'PLN', N, NRHS, M, -1 ) ) 00279 MAXWRK = MAX( MAXWRK, 2*M + M*NRHS ) 00280 END IF 00281 MINWRK = MAX( 2*M + N, 2*M + M*NRHS ) 00282 END IF 00283 END IF 00284 MINWRK = MIN( MINWRK, MAXWRK ) 00285 WORK( 1 ) = MAXWRK 00286 IWORK( 1 ) = LIWORK 00287 RWORK( 1 ) = LRWORK 00288 * 00289 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00290 INFO = -12 00291 END IF 00292 END IF 00293 * 00294 IF( INFO.NE.0 ) THEN 00295 CALL XERBLA( 'CGELSD', -INFO ) 00296 RETURN 00297 ELSE IF( LQUERY ) THEN 00298 RETURN 00299 END IF 00300 * 00301 * Quick return if possible. 00302 * 00303 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00304 RANK = 0 00305 RETURN 00306 END IF 00307 * 00308 * Get machine parameters. 00309 * 00310 EPS = SLAMCH( 'P' ) 00311 SFMIN = SLAMCH( 'S' ) 00312 SMLNUM = SFMIN / EPS 00313 BIGNUM = ONE / SMLNUM 00314 CALL SLABAD( SMLNUM, BIGNUM ) 00315 * 00316 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00317 * 00318 ANRM = CLANGE( 'M', M, N, A, LDA, RWORK ) 00319 IASCL = 0 00320 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00321 * 00322 * Scale matrix norm up to SMLNUM 00323 * 00324 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00325 IASCL = 1 00326 ELSE IF( ANRM.GT.BIGNUM ) THEN 00327 * 00328 * Scale matrix norm down to BIGNUM. 00329 * 00330 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00331 IASCL = 2 00332 ELSE IF( ANRM.EQ.ZERO ) THEN 00333 * 00334 * Matrix all zero. Return zero solution. 00335 * 00336 CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00337 CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00338 RANK = 0 00339 GO TO 10 00340 END IF 00341 * 00342 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00343 * 00344 BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00345 IBSCL = 0 00346 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00347 * 00348 * Scale matrix norm up to SMLNUM. 00349 * 00350 CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00351 IBSCL = 1 00352 ELSE IF( BNRM.GT.BIGNUM ) THEN 00353 * 00354 * Scale matrix norm down to BIGNUM. 00355 * 00356 CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00357 IBSCL = 2 00358 END IF 00359 * 00360 * If M < N make sure B(M+1:N,:) = 0 00361 * 00362 IF( M.LT.N ) 00363 $ CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00364 * 00365 * Overdetermined case. 00366 * 00367 IF( M.GE.N ) THEN 00368 * 00369 * Path 1 - overdetermined or exactly determined. 00370 * 00371 MM = M 00372 IF( M.GE.MNTHR ) THEN 00373 * 00374 * Path 1a - overdetermined, with many more rows than columns 00375 * 00376 MM = N 00377 ITAU = 1 00378 NWORK = ITAU + N 00379 * 00380 * Compute A=Q*R. 00381 * (RWorkspace: need N) 00382 * (CWorkspace: need N, prefer N*NB) 00383 * 00384 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00385 $ LWORK-NWORK+1, INFO ) 00386 * 00387 * Multiply B by transpose(Q). 00388 * (RWorkspace: need N) 00389 * (CWorkspace: need NRHS, prefer NRHS*NB) 00390 * 00391 CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00392 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00393 * 00394 * Zero out below R. 00395 * 00396 IF( N.GT.1 ) THEN 00397 CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00398 $ LDA ) 00399 END IF 00400 END IF 00401 * 00402 ITAUQ = 1 00403 ITAUP = ITAUQ + N 00404 NWORK = ITAUP + N 00405 IE = 1 00406 NRWORK = IE + N 00407 * 00408 * Bidiagonalize R in A. 00409 * (RWorkspace: need N) 00410 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB) 00411 * 00412 CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00413 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00414 $ INFO ) 00415 * 00416 * Multiply B by transpose of left bidiagonalizing vectors of R. 00417 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB) 00418 * 00419 CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00420 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00421 * 00422 * Solve the bidiagonal least squares problem. 00423 * 00424 CALL CLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB, 00425 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00426 $ IWORK, INFO ) 00427 IF( INFO.NE.0 ) THEN 00428 GO TO 10 00429 END IF 00430 * 00431 * Multiply B by right bidiagonalizing vectors of R. 00432 * 00433 CALL CUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00434 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00435 * 00436 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00437 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00438 * 00439 * Path 2a - underdetermined, with many more columns than rows 00440 * and sufficient workspace for an efficient algorithm. 00441 * 00442 LDWORK = M 00443 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00444 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00445 ITAU = 1 00446 NWORK = M + 1 00447 * 00448 * Compute A=L*Q. 00449 * (CWorkspace: need 2*M, prefer M+M*NB) 00450 * 00451 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00452 $ LWORK-NWORK+1, INFO ) 00453 IL = NWORK 00454 * 00455 * Copy L to WORK(IL), zeroing out above its diagonal. 00456 * 00457 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00458 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ), 00459 $ LDWORK ) 00460 ITAUQ = IL + LDWORK*M 00461 ITAUP = ITAUQ + M 00462 NWORK = ITAUP + M 00463 IE = 1 00464 NRWORK = IE + M 00465 * 00466 * Bidiagonalize L in WORK(IL). 00467 * (RWorkspace: need M) 00468 * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB) 00469 * 00470 CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ), 00471 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00472 $ LWORK-NWORK+1, INFO ) 00473 * 00474 * Multiply B by transpose of left bidiagonalizing vectors of L. 00475 * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00476 * 00477 CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK, 00478 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00479 $ LWORK-NWORK+1, INFO ) 00480 * 00481 * Solve the bidiagonal least squares problem. 00482 * 00483 CALL CLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00484 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00485 $ IWORK, INFO ) 00486 IF( INFO.NE.0 ) THEN 00487 GO TO 10 00488 END IF 00489 * 00490 * Multiply B by right bidiagonalizing vectors of L. 00491 * 00492 CALL CUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00493 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00494 $ LWORK-NWORK+1, INFO ) 00495 * 00496 * Zero out below first M rows of B. 00497 * 00498 CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00499 NWORK = ITAU + M 00500 * 00501 * Multiply transpose(Q) by B. 00502 * (CWorkspace: need NRHS, prefer NRHS*NB) 00503 * 00504 CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00505 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00506 * 00507 ELSE 00508 * 00509 * Path 2 - remaining underdetermined cases. 00510 * 00511 ITAUQ = 1 00512 ITAUP = ITAUQ + M 00513 NWORK = ITAUP + M 00514 IE = 1 00515 NRWORK = IE + M 00516 * 00517 * Bidiagonalize A. 00518 * (RWorkspace: need M) 00519 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 00520 * 00521 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00522 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00523 $ INFO ) 00524 * 00525 * Multiply B by transpose of left bidiagonalizing vectors. 00526 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB) 00527 * 00528 CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00529 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00530 * 00531 * Solve the bidiagonal least squares problem. 00532 * 00533 CALL CLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00534 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00535 $ IWORK, INFO ) 00536 IF( INFO.NE.0 ) THEN 00537 GO TO 10 00538 END IF 00539 * 00540 * Multiply B by right bidiagonalizing vectors of A. 00541 * 00542 CALL CUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00543 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00544 * 00545 END IF 00546 * 00547 * Undo scaling. 00548 * 00549 IF( IASCL.EQ.1 ) THEN 00550 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00551 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00552 $ INFO ) 00553 ELSE IF( IASCL.EQ.2 ) THEN 00554 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00555 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00556 $ INFO ) 00557 END IF 00558 IF( IBSCL.EQ.1 ) THEN 00559 CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00560 ELSE IF( IBSCL.EQ.2 ) THEN 00561 CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00562 END IF 00563 * 00564 10 CONTINUE 00565 WORK( 1 ) = MAXWRK 00566 IWORK( 1 ) = LIWORK 00567 RWORK( 1 ) = LRWORK 00568 RETURN 00569 * 00570 * End of CGELSD 00571 * 00572 END
1.7.3 
