|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, 00002 $ WORK, LWORK, RWORK, 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 REAL RWORK( * ), S( * ) 00015 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGELSS computes the minimum norm solution to a complex linear 00022 * least squares problem: 00023 * 00024 * Minimize 2-norm(| b - A*x |). 00025 * 00026 * using the singular value decomposition (SVD) of A. A is an M-by-N 00027 * matrix which may be rank-deficient. 00028 * 00029 * Several right hand side vectors b and solution vectors x can be 00030 * handled in a single call; they are stored as the columns of the 00031 * M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix 00032 * X. 00033 * 00034 * The effective rank of A is determined by treating as zero those 00035 * singular values which are less than RCOND times the largest singular 00036 * value. 00037 * 00038 * Arguments 00039 * ========= 00040 * 00041 * M (input) INTEGER 00042 * The number of rows of the matrix A. M >= 0. 00043 * 00044 * N (input) INTEGER 00045 * The number of columns of the matrix A. N >= 0. 00046 * 00047 * NRHS (input) INTEGER 00048 * The number of right hand sides, i.e., the number of columns 00049 * of the matrices B and X. NRHS >= 0. 00050 * 00051 * A (input/output) COMPLEX array, dimension (LDA,N) 00052 * On entry, the M-by-N matrix A. 00053 * On exit, the first min(m,n) rows of A are overwritten with 00054 * its right singular vectors, stored rowwise. 00055 * 00056 * LDA (input) INTEGER 00057 * The leading dimension of the array A. LDA >= max(1,M). 00058 * 00059 * B (input/output) COMPLEX array, dimension (LDB,NRHS) 00060 * On entry, the M-by-NRHS right hand side matrix B. 00061 * On exit, B is overwritten by the N-by-NRHS solution matrix X. 00062 * If m >= n and RANK = n, the residual sum-of-squares for 00063 * the solution in the i-th column is given by the sum of 00064 * squares of the modulus of elements n+1:m in that column. 00065 * 00066 * LDB (input) INTEGER 00067 * The leading dimension of the array B. LDB >= max(1,M,N). 00068 * 00069 * S (output) REAL array, dimension (min(M,N)) 00070 * The singular values of A in decreasing order. 00071 * The condition number of A in the 2-norm = S(1)/S(min(m,n)). 00072 * 00073 * RCOND (input) REAL 00074 * RCOND is used to determine the effective rank of A. 00075 * Singular values S(i) <= RCOND*S(1) are treated as zero. 00076 * If RCOND < 0, machine precision is used instead. 00077 * 00078 * RANK (output) INTEGER 00079 * The effective rank of A, i.e., the number of singular values 00080 * which are greater than RCOND*S(1). 00081 * 00082 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00083 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 * 00085 * LWORK (input) INTEGER 00086 * The dimension of the array WORK. LWORK >= 1, and also: 00087 * LWORK >= 2*min(M,N) + max(M,N,NRHS) 00088 * For good performance, LWORK should generally be larger. 00089 * 00090 * If LWORK = -1, then a workspace query is assumed; the routine 00091 * only calculates the optimal size of the WORK array, returns 00092 * this value as the first entry of the WORK array, and no error 00093 * message related to LWORK is issued by XERBLA. 00094 * 00095 * RWORK (workspace) REAL array, dimension (5*min(M,N)) 00096 * 00097 * INFO (output) INTEGER 00098 * = 0: successful exit 00099 * < 0: if INFO = -i, the i-th argument had an illegal value. 00100 * > 0: the algorithm for computing the SVD failed to converge; 00101 * if INFO = i, i off-diagonal elements of an intermediate 00102 * bidiagonal form did not converge to zero. 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 REAL ZERO, ONE 00108 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00109 COMPLEX CZERO, CONE 00110 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00111 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00112 * .. 00113 * .. Local Scalars .. 00114 LOGICAL LQUERY 00115 INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK, 00116 $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN, 00117 $ MAXWRK, MINMN, MINWRK, MM, MNTHR 00118 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR 00119 * .. 00120 * .. Local Arrays .. 00121 COMPLEX VDUM( 1 ) 00122 * .. 00123 * .. External Subroutines .. 00124 EXTERNAL CBDSQR, CCOPY, CGEBRD, CGELQF, CGEMM, CGEMV, 00125 $ CGEQRF, CLACPY, CLASCL, CLASET, CSRSCL, CUNGBR, 00126 $ CUNMBR, CUNMLQ, CUNMQR, SLABAD, SLASCL, SLASET, 00127 $ XERBLA 00128 * .. 00129 * .. External Functions .. 00130 INTEGER ILAENV 00131 REAL CLANGE, SLAMCH 00132 EXTERNAL ILAENV, CLANGE, SLAMCH 00133 * .. 00134 * .. Intrinsic Functions .. 00135 INTRINSIC MAX, MIN 00136 * .. 00137 * .. Executable Statements .. 00138 * 00139 * Test the input arguments 00140 * 00141 INFO = 0 00142 MINMN = MIN( M, N ) 00143 MAXMN = MAX( M, N ) 00144 LQUERY = ( LWORK.EQ.-1 ) 00145 IF( M.LT.0 ) THEN 00146 INFO = -1 00147 ELSE IF( N.LT.0 ) THEN 00148 INFO = -2 00149 ELSE IF( NRHS.LT.0 ) THEN 00150 INFO = -3 00151 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00152 INFO = -5 00153 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00154 INFO = -7 00155 END IF 00156 * 00157 * Compute workspace 00158 * (Note: Comments in the code beginning "Workspace:" describe the 00159 * minimal amount of workspace needed at that point in the code, 00160 * as well as the preferred amount for good performance. 00161 * CWorkspace refers to complex workspace, and RWorkspace refers 00162 * to real workspace. NB refers to the optimal block size for the 00163 * immediately following subroutine, as returned by ILAENV.) 00164 * 00165 IF( INFO.EQ.0 ) THEN 00166 MINWRK = 1 00167 MAXWRK = 1 00168 IF( MINMN.GT.0 ) THEN 00169 MM = M 00170 MNTHR = ILAENV( 6, 'CGELSS', ' ', M, N, NRHS, -1 ) 00171 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00172 * 00173 * Path 1a - overdetermined, with many more rows than 00174 * columns 00175 * 00176 MM = N 00177 MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'CGEQRF', ' ', M, 00178 $ N, -1, -1 ) ) 00179 MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'CUNMQR', 'LC', 00180 $ M, NRHS, N, -1 ) ) 00181 END IF 00182 IF( M.GE.N ) THEN 00183 * 00184 * Path 1 - overdetermined or exactly determined 00185 * 00186 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, 00187 $ 'CGEBRD', ' ', MM, N, -1, -1 ) ) 00188 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'CUNMBR', 00189 $ 'QLC', MM, NRHS, N, -1 ) ) 00190 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00191 $ 'CUNGBR', 'P', N, N, N, -1 ) ) 00192 MAXWRK = MAX( MAXWRK, N*NRHS ) 00193 MINWRK = 2*N + MAX( NRHS, M ) 00194 END IF 00195 IF( N.GT.M ) THEN 00196 MINWRK = 2*M + MAX( NRHS, N ) 00197 IF( N.GE.MNTHR ) THEN 00198 * 00199 * Path 2a - underdetermined, with many more columns 00200 * than rows 00201 * 00202 MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, 00203 $ -1 ) 00204 MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1, 00205 $ 'CGEBRD', ' ', M, M, -1, -1 ) ) 00206 MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1, 00207 $ 'CUNMBR', 'QLC', M, NRHS, M, -1 ) ) 00208 MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1, 00209 $ 'CUNGBR', 'P', M, M, M, -1 ) ) 00210 IF( NRHS.GT.1 ) THEN 00211 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00212 ELSE 00213 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00214 END IF 00215 MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'CUNMLQ', 00216 $ 'LC', N, NRHS, M, -1 ) ) 00217 ELSE 00218 * 00219 * Path 2 - underdetermined 00220 * 00221 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'CGEBRD', ' ', M, 00222 $ N, -1, -1 ) 00223 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'CUNMBR', 00224 $ 'QLC', M, NRHS, M, -1 ) ) 00225 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'CUNGBR', 00226 $ 'P', M, N, M, -1 ) ) 00227 MAXWRK = MAX( MAXWRK, N*NRHS ) 00228 END IF 00229 END IF 00230 MAXWRK = MAX( MINWRK, MAXWRK ) 00231 END IF 00232 WORK( 1 ) = MAXWRK 00233 * 00234 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) 00235 $ INFO = -12 00236 END IF 00237 * 00238 IF( INFO.NE.0 ) THEN 00239 CALL XERBLA( 'CGELSS', -INFO ) 00240 RETURN 00241 ELSE IF( LQUERY ) THEN 00242 RETURN 00243 END IF 00244 * 00245 * Quick return if possible 00246 * 00247 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00248 RANK = 0 00249 RETURN 00250 END IF 00251 * 00252 * Get machine parameters 00253 * 00254 EPS = SLAMCH( 'P' ) 00255 SFMIN = SLAMCH( 'S' ) 00256 SMLNUM = SFMIN / EPS 00257 BIGNUM = ONE / SMLNUM 00258 CALL SLABAD( SMLNUM, BIGNUM ) 00259 * 00260 * Scale A if max element outside range [SMLNUM,BIGNUM] 00261 * 00262 ANRM = CLANGE( 'M', M, N, A, LDA, RWORK ) 00263 IASCL = 0 00264 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00265 * 00266 * Scale matrix norm up to SMLNUM 00267 * 00268 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00269 IASCL = 1 00270 ELSE IF( ANRM.GT.BIGNUM ) THEN 00271 * 00272 * Scale matrix norm down to BIGNUM 00273 * 00274 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00275 IASCL = 2 00276 ELSE IF( ANRM.EQ.ZERO ) THEN 00277 * 00278 * Matrix all zero. Return zero solution. 00279 * 00280 CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00281 CALL SLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN ) 00282 RANK = 0 00283 GO TO 70 00284 END IF 00285 * 00286 * Scale B if max element outside range [SMLNUM,BIGNUM] 00287 * 00288 BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00289 IBSCL = 0 00290 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00291 * 00292 * Scale matrix norm up to SMLNUM 00293 * 00294 CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00295 IBSCL = 1 00296 ELSE IF( BNRM.GT.BIGNUM ) THEN 00297 * 00298 * Scale matrix norm down to BIGNUM 00299 * 00300 CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00301 IBSCL = 2 00302 END IF 00303 * 00304 * Overdetermined case 00305 * 00306 IF( M.GE.N ) THEN 00307 * 00308 * Path 1 - overdetermined or exactly determined 00309 * 00310 MM = M 00311 IF( M.GE.MNTHR ) THEN 00312 * 00313 * Path 1a - overdetermined, with many more rows than columns 00314 * 00315 MM = N 00316 ITAU = 1 00317 IWORK = ITAU + N 00318 * 00319 * Compute A=Q*R 00320 * (CWorkspace: need 2*N, prefer N+N*NB) 00321 * (RWorkspace: none) 00322 * 00323 CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00324 $ LWORK-IWORK+1, INFO ) 00325 * 00326 * Multiply B by transpose(Q) 00327 * (CWorkspace: need N+NRHS, prefer N+NRHS*NB) 00328 * (RWorkspace: none) 00329 * 00330 CALL CUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00331 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00332 * 00333 * Zero out below R 00334 * 00335 IF( N.GT.1 ) 00336 $ CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00337 $ LDA ) 00338 END IF 00339 * 00340 IE = 1 00341 ITAUQ = 1 00342 ITAUP = ITAUQ + N 00343 IWORK = ITAUP + N 00344 * 00345 * Bidiagonalize R in A 00346 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB) 00347 * (RWorkspace: need N) 00348 * 00349 CALL CGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00350 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00351 $ INFO ) 00352 * 00353 * Multiply B by transpose of left bidiagonalizing vectors of R 00354 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB) 00355 * (RWorkspace: none) 00356 * 00357 CALL CUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00358 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00359 * 00360 * Generate right bidiagonalizing vectors of R in A 00361 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 00362 * (RWorkspace: none) 00363 * 00364 CALL CUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00365 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00366 IRWORK = IE + N 00367 * 00368 * Perform bidiagonal QR iteration 00369 * multiply B by transpose of left singular vectors 00370 * compute right singular vectors in A 00371 * (CWorkspace: none) 00372 * (RWorkspace: need BDSPAC) 00373 * 00374 CALL CBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, 00375 $ 1, B, LDB, RWORK( IRWORK ), INFO ) 00376 IF( INFO.NE.0 ) 00377 $ GO TO 70 00378 * 00379 * Multiply B by reciprocals of singular values 00380 * 00381 THR = MAX( RCOND*S( 1 ), SFMIN ) 00382 IF( RCOND.LT.ZERO ) 00383 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00384 RANK = 0 00385 DO 10 I = 1, N 00386 IF( S( I ).GT.THR ) THEN 00387 CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00388 RANK = RANK + 1 00389 ELSE 00390 CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 00391 END IF 00392 10 CONTINUE 00393 * 00394 * Multiply B by right singular vectors 00395 * (CWorkspace: need N, prefer N*NRHS) 00396 * (RWorkspace: none) 00397 * 00398 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00399 CALL CGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB, 00400 $ CZERO, WORK, LDB ) 00401 CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 00402 ELSE IF( NRHS.GT.1 ) THEN 00403 CHUNK = LWORK / N 00404 DO 20 I = 1, NRHS, CHUNK 00405 BL = MIN( NRHS-I+1, CHUNK ) 00406 CALL CGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ), 00407 $ LDB, CZERO, WORK, N ) 00408 CALL CLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB ) 00409 20 CONTINUE 00410 ELSE 00411 CALL CGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 ) 00412 CALL CCOPY( N, WORK, 1, B, 1 ) 00413 END IF 00414 * 00415 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) ) 00416 $ THEN 00417 * 00418 * Underdetermined case, M much less than N 00419 * 00420 * Path 2a - underdetermined, with many more columns than rows 00421 * and sufficient workspace for an efficient algorithm 00422 * 00423 LDWORK = M 00424 IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) ) 00425 $ LDWORK = LDA 00426 ITAU = 1 00427 IWORK = M + 1 00428 * 00429 * Compute A=L*Q 00430 * (CWorkspace: need 2*M, prefer M+M*NB) 00431 * (RWorkspace: none) 00432 * 00433 CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00434 $ LWORK-IWORK+1, INFO ) 00435 IL = IWORK 00436 * 00437 * Copy L to WORK(IL), zeroing out above it 00438 * 00439 CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00440 CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ), 00441 $ LDWORK ) 00442 IE = 1 00443 ITAUQ = IL + LDWORK*M 00444 ITAUP = ITAUQ + M 00445 IWORK = ITAUP + M 00446 * 00447 * Bidiagonalize L in WORK(IL) 00448 * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB) 00449 * (RWorkspace: need M) 00450 * 00451 CALL CGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ), 00452 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ), 00453 $ LWORK-IWORK+1, INFO ) 00454 * 00455 * Multiply B by transpose of left bidiagonalizing vectors of L 00456 * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB) 00457 * (RWorkspace: none) 00458 * 00459 CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK, 00460 $ WORK( ITAUQ ), B, LDB, WORK( IWORK ), 00461 $ LWORK-IWORK+1, INFO ) 00462 * 00463 * Generate right bidiagonalizing vectors of R in WORK(IL) 00464 * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB) 00465 * (RWorkspace: none) 00466 * 00467 CALL CUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ), 00468 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00469 IRWORK = IE + M 00470 * 00471 * Perform bidiagonal QR iteration, computing right singular 00472 * vectors of L in WORK(IL) and multiplying B by transpose of 00473 * left singular vectors 00474 * (CWorkspace: need M*M) 00475 * (RWorkspace: need BDSPAC) 00476 * 00477 CALL CBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ), 00478 $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO ) 00479 IF( INFO.NE.0 ) 00480 $ GO TO 70 00481 * 00482 * Multiply B by reciprocals of singular values 00483 * 00484 THR = MAX( RCOND*S( 1 ), SFMIN ) 00485 IF( RCOND.LT.ZERO ) 00486 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00487 RANK = 0 00488 DO 30 I = 1, M 00489 IF( S( I ).GT.THR ) THEN 00490 CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00491 RANK = RANK + 1 00492 ELSE 00493 CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 00494 END IF 00495 30 CONTINUE 00496 IWORK = IL + M*LDWORK 00497 * 00498 * Multiply B by right singular vectors of L in WORK(IL) 00499 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS) 00500 * (RWorkspace: none) 00501 * 00502 IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN 00503 CALL CGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK, 00504 $ B, LDB, CZERO, WORK( IWORK ), LDB ) 00505 CALL CLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB ) 00506 ELSE IF( NRHS.GT.1 ) THEN 00507 CHUNK = ( LWORK-IWORK+1 ) / M 00508 DO 40 I = 1, NRHS, CHUNK 00509 BL = MIN( NRHS-I+1, CHUNK ) 00510 CALL CGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK, 00511 $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M ) 00512 CALL CLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ), 00513 $ LDB ) 00514 40 CONTINUE 00515 ELSE 00516 CALL CGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ), 00517 $ 1, CZERO, WORK( IWORK ), 1 ) 00518 CALL CCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 ) 00519 END IF 00520 * 00521 * Zero out below first M rows of B 00522 * 00523 CALL CLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00524 IWORK = ITAU + M 00525 * 00526 * Multiply transpose(Q) by B 00527 * (CWorkspace: need M+NRHS, prefer M+NHRS*NB) 00528 * (RWorkspace: none) 00529 * 00530 CALL CUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00531 $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00532 * 00533 ELSE 00534 * 00535 * Path 2 - remaining underdetermined cases 00536 * 00537 IE = 1 00538 ITAUQ = 1 00539 ITAUP = ITAUQ + M 00540 IWORK = ITAUP + M 00541 * 00542 * Bidiagonalize A 00543 * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB) 00544 * (RWorkspace: need N) 00545 * 00546 CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00547 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00548 $ INFO ) 00549 * 00550 * Multiply B by transpose of left bidiagonalizing vectors 00551 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB) 00552 * (RWorkspace: none) 00553 * 00554 CALL CUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00555 $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO ) 00556 * 00557 * Generate right bidiagonalizing vectors in A 00558 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 00559 * (RWorkspace: none) 00560 * 00561 CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 00562 $ WORK( IWORK ), LWORK-IWORK+1, INFO ) 00563 IRWORK = IE + M 00564 * 00565 * Perform bidiagonal QR iteration, 00566 * computing right singular vectors of A in A and 00567 * multiplying B by transpose of left singular vectors 00568 * (CWorkspace: none) 00569 * (RWorkspace: need BDSPAC) 00570 * 00571 CALL CBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM, 00572 $ 1, B, LDB, RWORK( IRWORK ), INFO ) 00573 IF( INFO.NE.0 ) 00574 $ GO TO 70 00575 * 00576 * Multiply B by reciprocals of singular values 00577 * 00578 THR = MAX( RCOND*S( 1 ), SFMIN ) 00579 IF( RCOND.LT.ZERO ) 00580 $ THR = MAX( EPS*S( 1 ), SFMIN ) 00581 RANK = 0 00582 DO 50 I = 1, M 00583 IF( S( I ).GT.THR ) THEN 00584 CALL CSRSCL( NRHS, S( I ), B( I, 1 ), LDB ) 00585 RANK = RANK + 1 00586 ELSE 00587 CALL CLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 00588 END IF 00589 50 CONTINUE 00590 * 00591 * Multiply B by right singular vectors of A 00592 * (CWorkspace: need N, prefer N*NRHS) 00593 * (RWorkspace: none) 00594 * 00595 IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN 00596 CALL CGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB, 00597 $ CZERO, WORK, LDB ) 00598 CALL CLACPY( 'G', N, NRHS, WORK, LDB, B, LDB ) 00599 ELSE IF( NRHS.GT.1 ) THEN 00600 CHUNK = LWORK / N 00601 DO 60 I = 1, NRHS, CHUNK 00602 BL = MIN( NRHS-I+1, CHUNK ) 00603 CALL CGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ), 00604 $ LDB, CZERO, WORK, N ) 00605 CALL CLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB ) 00606 60 CONTINUE 00607 ELSE 00608 CALL CGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 ) 00609 CALL CCOPY( N, WORK, 1, B, 1 ) 00610 END IF 00611 END IF 00612 * 00613 * Undo scaling 00614 * 00615 IF( IASCL.EQ.1 ) THEN 00616 CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00617 CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00618 $ INFO ) 00619 ELSE IF( IASCL.EQ.2 ) THEN 00620 CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00621 CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00622 $ INFO ) 00623 END IF 00624 IF( IBSCL.EQ.1 ) THEN 00625 CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00626 ELSE IF( IBSCL.EQ.2 ) THEN 00627 CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00628 END IF 00629 70 CONTINUE 00630 WORK( 1 ) = MAXWRK 00631 RETURN 00632 * 00633 * End of CGELSS 00634 * 00635 END
1.7.3 
