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