![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZGELSD 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 ZGELSD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelsd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelsd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelsd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGELSD( 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 * DOUBLE PRECISION RCOND 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION RWORK( * ), S( * ) 00031 * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZGELSD 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] A 00094 *> \verbatim 00095 *> A is COMPLEX*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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 complex16GEsolve 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 ZGELSD( 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 DOUBLE PRECISION RCOND 00236 * .. 00237 * .. Array Arguments .. 00238 INTEGER IWORK( * ) 00239 DOUBLE PRECISION RWORK( * ), S( * ) 00240 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00241 * .. 00242 * 00243 * ===================================================================== 00244 * 00245 * .. Parameters .. 00246 DOUBLE PRECISION ZERO, ONE, TWO 00247 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00248 COMPLEX*16 CZERO 00249 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+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 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM 00257 * .. 00258 * .. External Subroutines .. 00259 EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZGEBRD, ZGELQF, 00260 $ ZGEQRF, ZLACPY, ZLALSD, ZLASCL, ZLASET, ZUNMBR, 00261 $ ZUNMLQ, ZUNMQR 00262 * .. 00263 * .. External Functions .. 00264 INTEGER ILAENV 00265 DOUBLE PRECISION DLAMCH, ZLANGE 00266 EXTERNAL ILAENV, DLAMCH, ZLANGE 00267 * .. 00268 * .. Intrinsic Functions .. 00269 INTRINSIC INT, LOG, MAX, MIN, DBLE 00270 * .. 00271 * .. Executable Statements .. 00272 * 00273 * Test the input arguments. 00274 * 00275 INFO = 0 00276 MINMN = MIN( M, N ) 00277 MAXMN = MAX( M, N ) 00278 LQUERY = ( LWORK.EQ.-1 ) 00279 IF( M.LT.0 ) THEN 00280 INFO = -1 00281 ELSE IF( N.LT.0 ) THEN 00282 INFO = -2 00283 ELSE IF( NRHS.LT.0 ) THEN 00284 INFO = -3 00285 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00286 INFO = -5 00287 ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN 00288 INFO = -7 00289 END IF 00290 * 00291 * Compute workspace. 00292 * (Note: Comments in the code beginning "Workspace:" describe the 00293 * minimal amount of workspace needed at that point in the code, 00294 * as well as the preferred amount for good performance. 00295 * NB refers to the optimal block size for the immediately 00296 * following subroutine, as returned by ILAENV.) 00297 * 00298 IF( INFO.EQ.0 ) THEN 00299 MINWRK = 1 00300 MAXWRK = 1 00301 LIWORK = 1 00302 LRWORK = 1 00303 IF( MINMN.GT.0 ) THEN 00304 SMLSIZ = ILAENV( 9, 'ZGELSD', ' ', 0, 0, 0, 0 ) 00305 MNTHR = ILAENV( 6, 'ZGELSD', ' ', M, N, NRHS, -1 ) 00306 NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ + 1 ) ) / 00307 $ LOG( TWO ) ) + 1, 0 ) 00308 LIWORK = 3*MINMN*NLVL + 11*MINMN 00309 MM = M 00310 IF( M.GE.N .AND. M.GE.MNTHR ) THEN 00311 * 00312 * Path 1a - overdetermined, with many more rows than 00313 * columns. 00314 * 00315 MM = N 00316 MAXWRK = MAX( MAXWRK, N*ILAENV( 1, 'ZGEQRF', ' ', M, N, 00317 $ -1, -1 ) ) 00318 MAXWRK = MAX( MAXWRK, NRHS*ILAENV( 1, 'ZUNMQR', 'LC', M, 00319 $ NRHS, N, -1 ) ) 00320 END IF 00321 IF( M.GE.N ) THEN 00322 * 00323 * Path 1 - overdetermined or exactly determined. 00324 * 00325 LRWORK = 10*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + 00326 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00327 MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1, 00328 $ 'ZGEBRD', ' ', MM, N, -1, -1 ) ) 00329 MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR', 00330 $ 'QLC', MM, NRHS, N, -1 ) ) 00331 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 00332 $ 'ZUNMBR', 'PLN', N, NRHS, N, -1 ) ) 00333 MAXWRK = MAX( MAXWRK, 2*N + N*NRHS ) 00334 MINWRK = MAX( 2*N + MM, 2*N + N*NRHS ) 00335 END IF 00336 IF( N.GT.M ) THEN 00337 LRWORK = 10*M + 2*M*SMLSIZ + 8*M*NLVL + 3*SMLSIZ*NRHS + 00338 $ MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ) 00339 IF( N.GE.MNTHR ) THEN 00340 * 00341 * Path 2a - underdetermined, with many more columns 00342 * than rows. 00343 * 00344 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, 00345 $ -1 ) 00346 MAXWRK = MAX( MAXWRK, M*M + 4*M + 2*M*ILAENV( 1, 00347 $ 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00348 MAXWRK = MAX( MAXWRK, M*M + 4*M + NRHS*ILAENV( 1, 00349 $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) ) 00350 MAXWRK = MAX( MAXWRK, M*M + 4*M + ( M - 1 )*ILAENV( 1, 00351 $ 'ZUNMLQ', 'LC', N, NRHS, M, -1 ) ) 00352 IF( NRHS.GT.1 ) THEN 00353 MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS ) 00354 ELSE 00355 MAXWRK = MAX( MAXWRK, M*M + 2*M ) 00356 END IF 00357 MAXWRK = MAX( MAXWRK, M*M + 4*M + M*NRHS ) 00358 ! XXX: Ensure the Path 2a case below is triggered. The workspace 00359 ! calculation should use queries for all routines eventually. 00360 MAXWRK = MAX( MAXWRK, 00361 $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) ) 00362 ELSE 00363 * 00364 * Path 2 - underdetermined. 00365 * 00366 MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M, 00367 $ N, -1, -1 ) 00368 MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR', 00369 $ 'QLC', M, NRHS, M, -1 ) ) 00370 MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNMBR', 00371 $ 'PLN', N, NRHS, M, -1 ) ) 00372 MAXWRK = MAX( MAXWRK, 2*M + M*NRHS ) 00373 END IF 00374 MINWRK = MAX( 2*M + N, 2*M + M*NRHS ) 00375 END IF 00376 END IF 00377 MINWRK = MIN( MINWRK, MAXWRK ) 00378 WORK( 1 ) = MAXWRK 00379 IWORK( 1 ) = LIWORK 00380 RWORK( 1 ) = LRWORK 00381 * 00382 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00383 INFO = -12 00384 END IF 00385 END IF 00386 * 00387 IF( INFO.NE.0 ) THEN 00388 CALL XERBLA( 'ZGELSD', -INFO ) 00389 RETURN 00390 ELSE IF( LQUERY ) THEN 00391 RETURN 00392 END IF 00393 * 00394 * Quick return if possible. 00395 * 00396 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00397 RANK = 0 00398 RETURN 00399 END IF 00400 * 00401 * Get machine parameters. 00402 * 00403 EPS = DLAMCH( 'P' ) 00404 SFMIN = DLAMCH( 'S' ) 00405 SMLNUM = SFMIN / EPS 00406 BIGNUM = ONE / SMLNUM 00407 CALL DLABAD( SMLNUM, BIGNUM ) 00408 * 00409 * Scale A if max entry outside range [SMLNUM,BIGNUM]. 00410 * 00411 ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 00412 IASCL = 0 00413 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00414 * 00415 * Scale matrix norm up to SMLNUM 00416 * 00417 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 00418 IASCL = 1 00419 ELSE IF( ANRM.GT.BIGNUM ) THEN 00420 * 00421 * Scale matrix norm down to BIGNUM. 00422 * 00423 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 00424 IASCL = 2 00425 ELSE IF( ANRM.EQ.ZERO ) THEN 00426 * 00427 * Matrix all zero. Return zero solution. 00428 * 00429 CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 00430 CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 ) 00431 RANK = 0 00432 GO TO 10 00433 END IF 00434 * 00435 * Scale B if max entry outside range [SMLNUM,BIGNUM]. 00436 * 00437 BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 00438 IBSCL = 0 00439 IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 00440 * 00441 * Scale matrix norm up to SMLNUM. 00442 * 00443 CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 00444 IBSCL = 1 00445 ELSE IF( BNRM.GT.BIGNUM ) THEN 00446 * 00447 * Scale matrix norm down to BIGNUM. 00448 * 00449 CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 00450 IBSCL = 2 00451 END IF 00452 * 00453 * If M < N make sure B(M+1:N,:) = 0 00454 * 00455 IF( M.LT.N ) 00456 $ CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00457 * 00458 * Overdetermined case. 00459 * 00460 IF( M.GE.N ) THEN 00461 * 00462 * Path 1 - overdetermined or exactly determined. 00463 * 00464 MM = M 00465 IF( M.GE.MNTHR ) THEN 00466 * 00467 * Path 1a - overdetermined, with many more rows than columns 00468 * 00469 MM = N 00470 ITAU = 1 00471 NWORK = ITAU + N 00472 * 00473 * Compute A=Q*R. 00474 * (RWorkspace: need N) 00475 * (CWorkspace: need N, prefer N*NB) 00476 * 00477 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00478 $ LWORK-NWORK+1, INFO ) 00479 * 00480 * Multiply B by transpose(Q). 00481 * (RWorkspace: need N) 00482 * (CWorkspace: need NRHS, prefer NRHS*NB) 00483 * 00484 CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B, 00485 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00486 * 00487 * Zero out below R. 00488 * 00489 IF( N.GT.1 ) THEN 00490 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00491 $ LDA ) 00492 END IF 00493 END IF 00494 * 00495 ITAUQ = 1 00496 ITAUP = ITAUQ + N 00497 NWORK = ITAUP + N 00498 IE = 1 00499 NRWORK = IE + N 00500 * 00501 * Bidiagonalize R in A. 00502 * (RWorkspace: need N) 00503 * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB) 00504 * 00505 CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00506 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00507 $ INFO ) 00508 * 00509 * Multiply B by transpose of left bidiagonalizing vectors of R. 00510 * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB) 00511 * 00512 CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ), 00513 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00514 * 00515 * Solve the bidiagonal least squares problem. 00516 * 00517 CALL ZLALSD( 'U', SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB, 00518 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00519 $ IWORK, INFO ) 00520 IF( INFO.NE.0 ) THEN 00521 GO TO 10 00522 END IF 00523 * 00524 * Multiply B by right bidiagonalizing vectors of R. 00525 * 00526 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ), 00527 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00528 * 00529 ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+ 00530 $ MAX( M, 2*M-4, NRHS, N-3*M ) ) THEN 00531 * 00532 * Path 2a - underdetermined, with many more columns than rows 00533 * and sufficient workspace for an efficient algorithm. 00534 * 00535 LDWORK = M 00536 IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ), 00537 $ M*LDA+M+M*NRHS ) )LDWORK = LDA 00538 ITAU = 1 00539 NWORK = M + 1 00540 * 00541 * Compute A=L*Q. 00542 * (CWorkspace: need 2*M, prefer M+M*NB) 00543 * 00544 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00545 $ LWORK-NWORK+1, INFO ) 00546 IL = NWORK 00547 * 00548 * Copy L to WORK(IL), zeroing out above its diagonal. 00549 * 00550 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK ) 00551 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ), 00552 $ LDWORK ) 00553 ITAUQ = IL + LDWORK*M 00554 ITAUP = ITAUQ + M 00555 NWORK = ITAUP + M 00556 IE = 1 00557 NRWORK = IE + M 00558 * 00559 * Bidiagonalize L in WORK(IL). 00560 * (RWorkspace: need M) 00561 * (CWorkspace: need M*M+4*M, prefer M*M+4*M+2*M*NB) 00562 * 00563 CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ), 00564 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00565 $ LWORK-NWORK+1, INFO ) 00566 * 00567 * Multiply B by transpose of left bidiagonalizing vectors of L. 00568 * (CWorkspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) 00569 * 00570 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK, 00571 $ WORK( ITAUQ ), B, LDB, WORK( NWORK ), 00572 $ LWORK-NWORK+1, INFO ) 00573 * 00574 * Solve the bidiagonal least squares problem. 00575 * 00576 CALL ZLALSD( 'U', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00577 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00578 $ IWORK, INFO ) 00579 IF( INFO.NE.0 ) THEN 00580 GO TO 10 00581 END IF 00582 * 00583 * Multiply B by right bidiagonalizing vectors of L. 00584 * 00585 CALL ZUNMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK, 00586 $ WORK( ITAUP ), B, LDB, WORK( NWORK ), 00587 $ LWORK-NWORK+1, INFO ) 00588 * 00589 * Zero out below first M rows of B. 00590 * 00591 CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB ) 00592 NWORK = ITAU + M 00593 * 00594 * Multiply transpose(Q) by B. 00595 * (CWorkspace: need NRHS, prefer NRHS*NB) 00596 * 00597 CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B, 00598 $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00599 * 00600 ELSE 00601 * 00602 * Path 2 - remaining underdetermined cases. 00603 * 00604 ITAUQ = 1 00605 ITAUP = ITAUQ + M 00606 NWORK = ITAUP + M 00607 IE = 1 00608 NRWORK = IE + M 00609 * 00610 * Bidiagonalize A. 00611 * (RWorkspace: need M) 00612 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 00613 * 00614 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00615 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00616 $ INFO ) 00617 * 00618 * Multiply B by transpose of left bidiagonalizing vectors. 00619 * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB) 00620 * 00621 CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ), 00622 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00623 * 00624 * Solve the bidiagonal least squares problem. 00625 * 00626 CALL ZLALSD( 'L', SMLSIZ, M, NRHS, S, RWORK( IE ), B, LDB, 00627 $ RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ), 00628 $ IWORK, INFO ) 00629 IF( INFO.NE.0 ) THEN 00630 GO TO 10 00631 END IF 00632 * 00633 * Multiply B by right bidiagonalizing vectors of A. 00634 * 00635 CALL ZUNMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ), 00636 $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO ) 00637 * 00638 END IF 00639 * 00640 * Undo scaling. 00641 * 00642 IF( IASCL.EQ.1 ) THEN 00643 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 00644 CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 00645 $ INFO ) 00646 ELSE IF( IASCL.EQ.2 ) THEN 00647 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 00648 CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 00649 $ INFO ) 00650 END IF 00651 IF( IBSCL.EQ.1 ) THEN 00652 CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 00653 ELSE IF( IBSCL.EQ.2 ) THEN 00654 CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 00655 END IF 00656 * 00657 10 CONTINUE 00658 WORK( 1 ) = MAXWRK 00659 IWORK( 1 ) = LIWORK 00660 RWORK( 1 ) = LRWORK 00661 RETURN 00662 * 00663 * End of ZGELSD 00664 * 00665 END