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