![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZGESVD computes the singular value decomposition (SVD) 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 ZGESVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, 00022 * WORK, LWORK, RWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBU, JOBVT 00026 * INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION RWORK( * ), S( * ) 00030 * COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00031 * $ WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZGESVD computes the singular value decomposition (SVD) of a complex 00041 *> M-by-N matrix A, optionally computing the left and/or right singular 00042 *> vectors. The SVD is written 00043 *> 00044 *> A = U * SIGMA * conjugate-transpose(V) 00045 *> 00046 *> where SIGMA is an M-by-N matrix which is zero except for its 00047 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 00048 *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA 00049 *> are the singular values of A; they are real and non-negative, and 00050 *> are returned in descending order. The first min(m,n) columns of 00051 *> U and V are the left and right singular vectors of A. 00052 *> 00053 *> Note that the routine returns V**H, not V. 00054 *> \endverbatim 00055 * 00056 * Arguments: 00057 * ========== 00058 * 00059 *> \param[in] JOBU 00060 *> \verbatim 00061 *> JOBU is CHARACTER*1 00062 *> Specifies options for computing all or part of the matrix U: 00063 *> = 'A': all M columns of U are returned in array U: 00064 *> = 'S': the first min(m,n) columns of U (the left singular 00065 *> vectors) are returned in the array U; 00066 *> = 'O': the first min(m,n) columns of U (the left singular 00067 *> vectors) are overwritten on the array A; 00068 *> = 'N': no columns of U (no left singular vectors) are 00069 *> computed. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] JOBVT 00073 *> \verbatim 00074 *> JOBVT is CHARACTER*1 00075 *> Specifies options for computing all or part of the matrix 00076 *> V**H: 00077 *> = 'A': all N rows of V**H are returned in the array VT; 00078 *> = 'S': the first min(m,n) rows of V**H (the right singular 00079 *> vectors) are returned in the array VT; 00080 *> = 'O': the first min(m,n) rows of V**H (the right singular 00081 *> vectors) are overwritten on the array A; 00082 *> = 'N': no rows of V**H (no right singular vectors) are 00083 *> computed. 00084 *> 00085 *> JOBVT and JOBU cannot both be 'O'. 00086 *> \endverbatim 00087 *> 00088 *> \param[in] M 00089 *> \verbatim 00090 *> M is INTEGER 00091 *> The number of rows of the input matrix A. M >= 0. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] N 00095 *> \verbatim 00096 *> N is INTEGER 00097 *> The number of columns of the input matrix A. N >= 0. 00098 *> \endverbatim 00099 *> 00100 *> \param[in,out] A 00101 *> \verbatim 00102 *> A is COMPLEX*16 array, dimension (LDA,N) 00103 *> On entry, the M-by-N matrix A. 00104 *> On exit, 00105 *> if JOBU = 'O', A is overwritten with the first min(m,n) 00106 *> columns of U (the left singular vectors, 00107 *> stored columnwise); 00108 *> if JOBVT = 'O', A is overwritten with the first min(m,n) 00109 *> rows of V**H (the right singular vectors, 00110 *> stored rowwise); 00111 *> if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A 00112 *> are destroyed. 00113 *> \endverbatim 00114 *> 00115 *> \param[in] LDA 00116 *> \verbatim 00117 *> LDA is INTEGER 00118 *> The leading dimension of the array A. LDA >= max(1,M). 00119 *> \endverbatim 00120 *> 00121 *> \param[out] S 00122 *> \verbatim 00123 *> S is DOUBLE PRECISION array, dimension (min(M,N)) 00124 *> The singular values of A, sorted so that S(i) >= S(i+1). 00125 *> \endverbatim 00126 *> 00127 *> \param[out] U 00128 *> \verbatim 00129 *> U is COMPLEX*16 array, dimension (LDU,UCOL) 00130 *> (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. 00131 *> If JOBU = 'A', U contains the M-by-M unitary matrix U; 00132 *> if JOBU = 'S', U contains the first min(m,n) columns of U 00133 *> (the left singular vectors, stored columnwise); 00134 *> if JOBU = 'N' or 'O', U is not referenced. 00135 *> \endverbatim 00136 *> 00137 *> \param[in] LDU 00138 *> \verbatim 00139 *> LDU is INTEGER 00140 *> The leading dimension of the array U. LDU >= 1; if 00141 *> JOBU = 'S' or 'A', LDU >= M. 00142 *> \endverbatim 00143 *> 00144 *> \param[out] VT 00145 *> \verbatim 00146 *> VT is COMPLEX*16 array, dimension (LDVT,N) 00147 *> If JOBVT = 'A', VT contains the N-by-N unitary matrix 00148 *> V**H; 00149 *> if JOBVT = 'S', VT contains the first min(m,n) rows of 00150 *> V**H (the right singular vectors, stored rowwise); 00151 *> if JOBVT = 'N' or 'O', VT is not referenced. 00152 *> \endverbatim 00153 *> 00154 *> \param[in] LDVT 00155 *> \verbatim 00156 *> LDVT is INTEGER 00157 *> The leading dimension of the array VT. LDVT >= 1; if 00158 *> JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). 00159 *> \endverbatim 00160 *> 00161 *> \param[out] WORK 00162 *> \verbatim 00163 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 00164 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] LWORK 00168 *> \verbatim 00169 *> LWORK is INTEGER 00170 *> The dimension of the array WORK. 00171 *> LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). 00172 *> For good performance, LWORK should generally be larger. 00173 *> 00174 *> If LWORK = -1, then a workspace query is assumed; the routine 00175 *> only calculates the optimal size of the WORK array, returns 00176 *> this value as the first entry of the WORK array, and no error 00177 *> message related to LWORK is issued by XERBLA. 00178 *> \endverbatim 00179 *> 00180 *> \param[out] RWORK 00181 *> \verbatim 00182 *> RWORK is DOUBLE PRECISION array, dimension (5*min(M,N)) 00183 *> On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the 00184 *> unconverged superdiagonal elements of an upper bidiagonal 00185 *> matrix B whose diagonal is in S (not necessarily sorted). 00186 *> B satisfies A = U * B * VT, so it has the same singular 00187 *> values as A, and singular vectors related by U and VT. 00188 *> \endverbatim 00189 *> 00190 *> \param[out] INFO 00191 *> \verbatim 00192 *> INFO is INTEGER 00193 *> = 0: successful exit. 00194 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00195 *> > 0: if ZBDSQR did not converge, INFO specifies how many 00196 *> superdiagonals of an intermediate bidiagonal form B 00197 *> did not converge to zero. See the description of RWORK 00198 *> above for details. 00199 *> \endverbatim 00200 * 00201 * Authors: 00202 * ======== 00203 * 00204 *> \author Univ. of Tennessee 00205 *> \author Univ. of California Berkeley 00206 *> \author Univ. of Colorado Denver 00207 *> \author NAG Ltd. 00208 * 00209 *> \date April 2012 00210 * 00211 *> \ingroup complex16GEsing 00212 * 00213 * ===================================================================== 00214 SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, 00215 $ VT, LDVT, WORK, LWORK, RWORK, INFO ) 00216 * 00217 * -- LAPACK driver routine (version 3.4.1) -- 00218 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00220 * April 2012 00221 * 00222 * .. Scalar Arguments .. 00223 CHARACTER JOBU, JOBVT 00224 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00225 * .. 00226 * .. Array Arguments .. 00227 DOUBLE PRECISION RWORK( * ), S( * ) 00228 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00229 $ WORK( * ) 00230 * .. 00231 * 00232 * ===================================================================== 00233 * 00234 * .. Parameters .. 00235 COMPLEX*16 CZERO, CONE 00236 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), 00237 $ CONE = ( 1.0D0, 0.0D0 ) ) 00238 DOUBLE PRECISION ZERO, ONE 00239 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00240 * .. 00241 * .. Local Scalars .. 00242 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS, 00243 $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS 00244 INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL, 00245 $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU, 00246 $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU, 00247 $ NRVT, WRKBL 00248 INTEGER LWORK_ZGEQRF, LWORK_ZUNGQR_N, LWORK_ZUNGQR_M, 00249 $ LWORK_ZGEBRD, LWORK_ZUNGBR_P, LWORK_ZUNGBR_Q, 00250 $ LWORK_ZGELQF, LWORK_ZUNGLQ_N, LWORK_ZUNGLQ_M 00251 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM 00252 * .. 00253 * .. Local Arrays .. 00254 DOUBLE PRECISION DUM( 1 ) 00255 COMPLEX*16 CDUM( 1 ) 00256 * .. 00257 * .. External Subroutines .. 00258 EXTERNAL DLASCL, XERBLA, ZBDSQR, ZGEBRD, ZGELQF, ZGEMM, 00259 $ ZGEQRF, ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNGLQ, 00260 $ ZUNGQR, ZUNMBR 00261 * .. 00262 * .. External Functions .. 00263 LOGICAL LSAME 00264 INTEGER ILAENV 00265 DOUBLE PRECISION DLAMCH, ZLANGE 00266 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 00267 * .. 00268 * .. Intrinsic Functions .. 00269 INTRINSIC MAX, MIN, SQRT 00270 * .. 00271 * .. Executable Statements .. 00272 * 00273 * Test the input arguments 00274 * 00275 INFO = 0 00276 MINMN = MIN( M, N ) 00277 WNTUA = LSAME( JOBU, 'A' ) 00278 WNTUS = LSAME( JOBU, 'S' ) 00279 WNTUAS = WNTUA .OR. WNTUS 00280 WNTUO = LSAME( JOBU, 'O' ) 00281 WNTUN = LSAME( JOBU, 'N' ) 00282 WNTVA = LSAME( JOBVT, 'A' ) 00283 WNTVS = LSAME( JOBVT, 'S' ) 00284 WNTVAS = WNTVA .OR. WNTVS 00285 WNTVO = LSAME( JOBVT, 'O' ) 00286 WNTVN = LSAME( JOBVT, 'N' ) 00287 LQUERY = ( LWORK.EQ.-1 ) 00288 * 00289 IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN 00290 INFO = -1 00291 ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR. 00292 $ ( WNTVO .AND. WNTUO ) ) THEN 00293 INFO = -2 00294 ELSE IF( M.LT.0 ) THEN 00295 INFO = -3 00296 ELSE IF( N.LT.0 ) THEN 00297 INFO = -4 00298 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00299 INFO = -6 00300 ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN 00301 INFO = -9 00302 ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR. 00303 $ ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN 00304 INFO = -11 00305 END IF 00306 * 00307 * Compute workspace 00308 * (Note: Comments in the code beginning "Workspace:" describe the 00309 * minimal amount of workspace needed at that point in the code, 00310 * as well as the preferred amount for good performance. 00311 * CWorkspace refers to complex workspace, and RWorkspace to 00312 * real workspace. NB refers to the optimal block size for the 00313 * immediately following subroutine, as returned by ILAENV.) 00314 * 00315 IF( INFO.EQ.0 ) THEN 00316 MINWRK = 1 00317 MAXWRK = 1 00318 IF( M.GE.N .AND. MINMN.GT.0 ) THEN 00319 * 00320 * Space needed for ZBDSQR is BDSPAC = 5*N 00321 * 00322 MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) 00323 * Compute space needed for ZGEQRF 00324 CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00325 LWORK_ZGEQRF=DUM(1) 00326 * Compute space needed for ZUNGQR 00327 CALL ZUNGQR( M, N, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00328 LWORK_ZUNGQR_N=DUM(1) 00329 CALL ZUNGQR( M, M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00330 LWORK_ZUNGQR_M=DUM(1) 00331 * Compute space needed for ZGEBRD 00332 CALL ZGEBRD( N, N, A, LDA, S, DUM(1), DUM(1), 00333 $ DUM(1), DUM(1), -1, IERR ) 00334 LWORK_ZGEBRD=DUM(1) 00335 * Compute space needed for ZUNGBR 00336 CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1), 00337 $ DUM(1), -1, IERR ) 00338 LWORK_ZUNGBR_P=DUM(1) 00339 CALL ZUNGBR( 'Q', N, N, N, A, LDA, DUM(1), 00340 $ DUM(1), -1, IERR ) 00341 LWORK_ZUNGBR_Q=DUM(1) 00342 * 00343 IF( M.GE.MNTHR ) THEN 00344 IF( WNTUN ) THEN 00345 * 00346 * Path 1 (M much larger than N, JOBU='N') 00347 * 00348 MAXWRK = N + LWORK_ZGEQRF 00349 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZGEBRD ) 00350 IF( WNTVO .OR. WNTVAS ) 00351 $ MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) 00352 MINWRK = 3*N 00353 ELSE IF( WNTUO .AND. WNTVN ) THEN 00354 * 00355 * Path 2 (M much larger than N, JOBU='O', JOBVT='N') 00356 * 00357 WRKBL = N + LWORK_ZGEQRF 00358 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 00359 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00360 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00361 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 00362 MINWRK = 2*N + M 00363 ELSE IF( WNTUO .AND. WNTVAS ) THEN 00364 * 00365 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 00366 * 'A') 00367 * 00368 WRKBL = N + LWORK_ZGEQRF 00369 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 00370 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00371 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00372 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 00373 MAXWRK = MAX( N*N+WRKBL, N*N+M*N ) 00374 MINWRK = 2*N + M 00375 ELSE IF( WNTUS .AND. WNTVN ) THEN 00376 * 00377 * Path 4 (M much larger than N, JOBU='S', JOBVT='N') 00378 * 00379 WRKBL = N + LWORK_ZGEQRF 00380 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 00381 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00382 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00383 MAXWRK = N*N + WRKBL 00384 MINWRK = 2*N + M 00385 ELSE IF( WNTUS .AND. WNTVO ) THEN 00386 * 00387 * Path 5 (M much larger than N, JOBU='S', JOBVT='O') 00388 * 00389 WRKBL = N + LWORK_ZGEQRF 00390 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 00391 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00392 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00393 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 00394 MAXWRK = 2*N*N + WRKBL 00395 MINWRK = 2*N + M 00396 ELSE IF( WNTUS .AND. WNTVAS ) THEN 00397 * 00398 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or 00399 * 'A') 00400 * 00401 WRKBL = N + LWORK_ZGEQRF 00402 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_N ) 00403 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00404 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00405 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 00406 MAXWRK = N*N + WRKBL 00407 MINWRK = 2*N + M 00408 ELSE IF( WNTUA .AND. WNTVN ) THEN 00409 * 00410 * Path 7 (M much larger than N, JOBU='A', JOBVT='N') 00411 * 00412 WRKBL = N + LWORK_ZGEQRF 00413 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 00414 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00415 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00416 MAXWRK = N*N + WRKBL 00417 MINWRK = 2*N + M 00418 ELSE IF( WNTUA .AND. WNTVO ) THEN 00419 * 00420 * Path 8 (M much larger than N, JOBU='A', JOBVT='O') 00421 * 00422 WRKBL = N + LWORK_ZGEQRF 00423 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 00424 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00425 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00426 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 00427 MAXWRK = 2*N*N + WRKBL 00428 MINWRK = 2*N + M 00429 ELSE IF( WNTUA .AND. WNTVAS ) THEN 00430 * 00431 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or 00432 * 'A') 00433 * 00434 WRKBL = N + LWORK_ZGEQRF 00435 WRKBL = MAX( WRKBL, N+LWORK_ZUNGQR_M ) 00436 WRKBL = MAX( WRKBL, 2*N+LWORK_ZGEBRD ) 00437 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_Q ) 00438 WRKBL = MAX( WRKBL, 2*N+LWORK_ZUNGBR_P ) 00439 MAXWRK = N*N + WRKBL 00440 MINWRK = 2*N + M 00441 END IF 00442 ELSE 00443 * 00444 * Path 10 (M at least N, but not much larger) 00445 * 00446 CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00447 $ DUM(1), DUM(1), -1, IERR ) 00448 LWORK_ZGEBRD=DUM(1) 00449 MAXWRK = 2*N + LWORK_ZGEBRD 00450 IF( WNTUS .OR. WNTUO ) THEN 00451 CALL ZUNGBR( 'Q', M, N, N, A, LDA, DUM(1), 00452 $ DUM(1), -1, IERR ) 00453 LWORK_ZUNGBR_Q=DUM(1) 00454 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) 00455 END IF 00456 IF( WNTUA ) THEN 00457 CALL ZUNGBR( 'Q', M, M, N, A, LDA, DUM(1), 00458 $ DUM(1), -1, IERR ) 00459 LWORK_ZUNGBR_Q=DUM(1) 00460 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_Q ) 00461 END IF 00462 IF( .NOT.WNTVN ) THEN 00463 MAXWRK = MAX( MAXWRK, 2*N+LWORK_ZUNGBR_P ) 00464 MINWRK = 2*N + M 00465 END IF 00466 END IF 00467 ELSE IF( MINMN.GT.0 ) THEN 00468 * 00469 * Space needed for ZBDSQR is BDSPAC = 5*M 00470 * 00471 MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 ) 00472 * Compute space needed for ZGELQF 00473 CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1), -1, IERR ) 00474 LWORK_ZGELQF=DUM(1) 00475 * Compute space needed for ZUNGLQ 00476 CALL ZUNGLQ( N, N, M, DUM(1), N, DUM(1), DUM(1), -1, IERR ) 00477 LWORK_ZUNGLQ_N=DUM(1) 00478 CALL ZUNGLQ( M, N, M, A, LDA, DUM(1), DUM(1), -1, IERR ) 00479 LWORK_ZUNGLQ_M=DUM(1) 00480 * Compute space needed for ZGEBRD 00481 CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1), 00482 $ DUM(1), DUM(1), -1, IERR ) 00483 LWORK_ZGEBRD=DUM(1) 00484 * Compute space needed for ZUNGBR P 00485 CALL ZUNGBR( 'P', M, M, M, A, N, DUM(1), 00486 $ DUM(1), -1, IERR ) 00487 LWORK_ZUNGBR_P=DUM(1) 00488 * Compute space needed for ZUNGBR Q 00489 CALL ZUNGBR( 'Q', M, M, M, A, N, DUM(1), 00490 $ DUM(1), -1, IERR ) 00491 LWORK_ZUNGBR_Q=DUM(1) 00492 IF( N.GE.MNTHR ) THEN 00493 IF( WNTVN ) THEN 00494 * 00495 * Path 1t(N much larger than M, JOBVT='N') 00496 * 00497 MAXWRK = M + LWORK_ZGELQF 00498 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZGEBRD ) 00499 IF( WNTUO .OR. WNTUAS ) 00500 $ MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) 00501 MINWRK = 3*M 00502 ELSE IF( WNTVO .AND. WNTUN ) THEN 00503 * 00504 * Path 2t(N much larger than M, JOBU='N', JOBVT='O') 00505 * 00506 WRKBL = M + LWORK_ZGELQF 00507 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 00508 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00509 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00510 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 00511 MINWRK = 2*M + N 00512 ELSE IF( WNTVO .AND. WNTUAS ) THEN 00513 * 00514 * Path 3t(N much larger than M, JOBU='S' or 'A', 00515 * JOBVT='O') 00516 * 00517 WRKBL = M + LWORK_ZGELQF 00518 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 00519 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00520 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00521 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 00522 MAXWRK = MAX( M*M+WRKBL, M*M+M*N ) 00523 MINWRK = 2*M + N 00524 ELSE IF( WNTVS .AND. WNTUN ) THEN 00525 * 00526 * Path 4t(N much larger than M, JOBU='N', JOBVT='S') 00527 * 00528 WRKBL = M + LWORK_ZGELQF 00529 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 00530 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00531 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00532 MAXWRK = M*M + WRKBL 00533 MINWRK = 2*M + N 00534 ELSE IF( WNTVS .AND. WNTUO ) THEN 00535 * 00536 * Path 5t(N much larger than M, JOBU='O', JOBVT='S') 00537 * 00538 WRKBL = M + LWORK_ZGELQF 00539 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 00540 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00541 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00542 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 00543 MAXWRK = 2*M*M + WRKBL 00544 MINWRK = 2*M + N 00545 ELSE IF( WNTVS .AND. WNTUAS ) THEN 00546 * 00547 * Path 6t(N much larger than M, JOBU='S' or 'A', 00548 * JOBVT='S') 00549 * 00550 WRKBL = M + LWORK_ZGELQF 00551 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_M ) 00552 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00553 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00554 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 00555 MAXWRK = M*M + WRKBL 00556 MINWRK = 2*M + N 00557 ELSE IF( WNTVA .AND. WNTUN ) THEN 00558 * 00559 * Path 7t(N much larger than M, JOBU='N', JOBVT='A') 00560 * 00561 WRKBL = M + LWORK_ZGELQF 00562 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 00563 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00564 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00565 MAXWRK = M*M + WRKBL 00566 MINWRK = 2*M + N 00567 ELSE IF( WNTVA .AND. WNTUO ) THEN 00568 * 00569 * Path 8t(N much larger than M, JOBU='O', JOBVT='A') 00570 * 00571 WRKBL = M + LWORK_ZGELQF 00572 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 00573 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00574 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00575 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 00576 MAXWRK = 2*M*M + WRKBL 00577 MINWRK = 2*M + N 00578 ELSE IF( WNTVA .AND. WNTUAS ) THEN 00579 * 00580 * Path 9t(N much larger than M, JOBU='S' or 'A', 00581 * JOBVT='A') 00582 * 00583 WRKBL = M + LWORK_ZGELQF 00584 WRKBL = MAX( WRKBL, M+LWORK_ZUNGLQ_N ) 00585 WRKBL = MAX( WRKBL, 2*M+LWORK_ZGEBRD ) 00586 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_P ) 00587 WRKBL = MAX( WRKBL, 2*M+LWORK_ZUNGBR_Q ) 00588 MAXWRK = M*M + WRKBL 00589 MINWRK = 2*M + N 00590 END IF 00591 ELSE 00592 * 00593 * Path 10t(N greater than M, but not much larger) 00594 * 00595 CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1), 00596 $ DUM(1), DUM(1), -1, IERR ) 00597 LWORK_ZGEBRD=DUM(1) 00598 MAXWRK = 2*M + LWORK_ZGEBRD 00599 IF( WNTVS .OR. WNTVO ) THEN 00600 * Compute space needed for ZUNGBR P 00601 CALL ZUNGBR( 'P', M, N, M, A, N, DUM(1), 00602 $ DUM(1), -1, IERR ) 00603 LWORK_ZUNGBR_P=DUM(1) 00604 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) 00605 END IF 00606 IF( WNTVA ) THEN 00607 CALL ZUNGBR( 'P', N, N, M, A, N, DUM(1), 00608 $ DUM(1), -1, IERR ) 00609 LWORK_ZUNGBR_P=DUM(1) 00610 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_P ) 00611 END IF 00612 IF( .NOT.WNTUN ) THEN 00613 MAXWRK = MAX( MAXWRK, 2*M+LWORK_ZUNGBR_Q ) 00614 MINWRK = 2*M + N 00615 END IF 00616 END IF 00617 END IF 00618 MAXWRK = MAX( MAXWRK, MINWRK ) 00619 WORK( 1 ) = MAXWRK 00620 * 00621 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00622 INFO = -13 00623 END IF 00624 END IF 00625 * 00626 IF( INFO.NE.0 ) THEN 00627 CALL XERBLA( 'ZGESVD', -INFO ) 00628 RETURN 00629 ELSE IF( LQUERY ) THEN 00630 RETURN 00631 END IF 00632 * 00633 * Quick return if possible 00634 * 00635 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00636 RETURN 00637 END IF 00638 * 00639 * Get machine constants 00640 * 00641 EPS = DLAMCH( 'P' ) 00642 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS 00643 BIGNUM = ONE / SMLNUM 00644 * 00645 * Scale A if max element outside range [SMLNUM,BIGNUM] 00646 * 00647 ANRM = ZLANGE( 'M', M, N, A, LDA, DUM ) 00648 ISCL = 0 00649 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00650 ISCL = 1 00651 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 00652 ELSE IF( ANRM.GT.BIGNUM ) THEN 00653 ISCL = 1 00654 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 00655 END IF 00656 * 00657 IF( M.GE.N ) THEN 00658 * 00659 * A has at least as many rows as columns. If A has sufficiently 00660 * more rows than columns, first reduce using the QR 00661 * decomposition (if sufficient workspace available) 00662 * 00663 IF( M.GE.MNTHR ) THEN 00664 * 00665 IF( WNTUN ) THEN 00666 * 00667 * Path 1 (M much larger than N, JOBU='N') 00668 * No left singular vectors to be computed 00669 * 00670 ITAU = 1 00671 IWORK = ITAU + N 00672 * 00673 * Compute A=Q*R 00674 * (CWorkspace: need 2*N, prefer N+N*NB) 00675 * (RWorkspace: need 0) 00676 * 00677 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 00678 $ LWORK-IWORK+1, IERR ) 00679 * 00680 * Zero out below R 00681 * 00682 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00683 $ LDA ) 00684 IE = 1 00685 ITAUQ = 1 00686 ITAUP = ITAUQ + N 00687 IWORK = ITAUP + N 00688 * 00689 * Bidiagonalize R in A 00690 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00691 * (RWorkspace: need N) 00692 * 00693 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00694 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 00695 $ IERR ) 00696 NCVT = 0 00697 IF( WNTVO .OR. WNTVAS ) THEN 00698 * 00699 * If right singular vectors desired, generate P'. 00700 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 00701 * (RWorkspace: 0) 00702 * 00703 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 00704 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00705 NCVT = N 00706 END IF 00707 IRWORK = IE + N 00708 * 00709 * Perform bidiagonal QR iteration, computing right 00710 * singular vectors of A in A if desired 00711 * (CWorkspace: 0) 00712 * (RWorkspace: need BDSPAC) 00713 * 00714 CALL ZBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA, 00715 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 00716 * 00717 * If right singular vectors desired in VT, copy them there 00718 * 00719 IF( WNTVAS ) 00720 $ CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT ) 00721 * 00722 ELSE IF( WNTUO .AND. WNTVN ) THEN 00723 * 00724 * Path 2 (M much larger than N, JOBU='O', JOBVT='N') 00725 * N left singular vectors to be overwritten on A and 00726 * no right singular vectors to be computed 00727 * 00728 IF( LWORK.GE.N*N+3*N ) THEN 00729 * 00730 * Sufficient workspace for a fast algorithm 00731 * 00732 IR = 1 00733 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 00734 * 00735 * WORK(IU) is LDA by N, WORK(IR) is LDA by N 00736 * 00737 LDWRKU = LDA 00738 LDWRKR = LDA 00739 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 00740 * 00741 * WORK(IU) is LDA by N, WORK(IR) is N by N 00742 * 00743 LDWRKU = LDA 00744 LDWRKR = N 00745 ELSE 00746 * 00747 * WORK(IU) is LDWRKU by N, WORK(IR) is N by N 00748 * 00749 LDWRKU = ( LWORK-N*N ) / N 00750 LDWRKR = N 00751 END IF 00752 ITAU = IR + LDWRKR*N 00753 IWORK = ITAU + N 00754 * 00755 * Compute A=Q*R 00756 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00757 * (RWorkspace: 0) 00758 * 00759 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 00760 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00761 * 00762 * Copy R to WORK(IR) and zero out below it 00763 * 00764 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00765 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 00766 $ WORK( IR+1 ), LDWRKR ) 00767 * 00768 * Generate Q in A 00769 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00770 * (RWorkspace: 0) 00771 * 00772 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00773 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00774 IE = 1 00775 ITAUQ = ITAU 00776 ITAUP = ITAUQ + N 00777 IWORK = ITAUP + N 00778 * 00779 * Bidiagonalize R in WORK(IR) 00780 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 00781 * (RWorkspace: need N) 00782 * 00783 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 00784 $ WORK( ITAUQ ), WORK( ITAUP ), 00785 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00786 * 00787 * Generate left vectors bidiagonalizing R 00788 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00789 * (RWorkspace: need 0) 00790 * 00791 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 00792 $ WORK( ITAUQ ), WORK( IWORK ), 00793 $ LWORK-IWORK+1, IERR ) 00794 IRWORK = IE + N 00795 * 00796 * Perform bidiagonal QR iteration, computing left 00797 * singular vectors of R in WORK(IR) 00798 * (CWorkspace: need N*N) 00799 * (RWorkspace: need BDSPAC) 00800 * 00801 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1, 00802 $ WORK( IR ), LDWRKR, CDUM, 1, 00803 $ RWORK( IRWORK ), INFO ) 00804 IU = ITAUQ 00805 * 00806 * Multiply Q in A by left singular vectors of R in 00807 * WORK(IR), storing result in WORK(IU) and copying to A 00808 * (CWorkspace: need N*N+N, prefer N*N+M*N) 00809 * (RWorkspace: 0) 00810 * 00811 DO 10 I = 1, M, LDWRKU 00812 CHUNK = MIN( M-I+1, LDWRKU ) 00813 CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00814 $ LDA, WORK( IR ), LDWRKR, CZERO, 00815 $ WORK( IU ), LDWRKU ) 00816 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 00817 $ A( I, 1 ), LDA ) 00818 10 CONTINUE 00819 * 00820 ELSE 00821 * 00822 * Insufficient workspace for a fast algorithm 00823 * 00824 IE = 1 00825 ITAUQ = 1 00826 ITAUP = ITAUQ + N 00827 IWORK = ITAUP + N 00828 * 00829 * Bidiagonalize A 00830 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 00831 * (RWorkspace: N) 00832 * 00833 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), 00834 $ WORK( ITAUQ ), WORK( ITAUP ), 00835 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00836 * 00837 * Generate left vectors bidiagonalizing A 00838 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 00839 * (RWorkspace: 0) 00840 * 00841 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 00842 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00843 IRWORK = IE + N 00844 * 00845 * Perform bidiagonal QR iteration, computing left 00846 * singular vectors of A in A 00847 * (CWorkspace: need 0) 00848 * (RWorkspace: need BDSPAC) 00849 * 00850 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1, 00851 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 00852 * 00853 END IF 00854 * 00855 ELSE IF( WNTUO .AND. WNTVAS ) THEN 00856 * 00857 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A') 00858 * N left singular vectors to be overwritten on A and 00859 * N right singular vectors to be computed in VT 00860 * 00861 IF( LWORK.GE.N*N+3*N ) THEN 00862 * 00863 * Sufficient workspace for a fast algorithm 00864 * 00865 IR = 1 00866 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN 00867 * 00868 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 00869 * 00870 LDWRKU = LDA 00871 LDWRKR = LDA 00872 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN 00873 * 00874 * WORK(IU) is LDA by N and WORK(IR) is N by N 00875 * 00876 LDWRKU = LDA 00877 LDWRKR = N 00878 ELSE 00879 * 00880 * WORK(IU) is LDWRKU by N and WORK(IR) is N by N 00881 * 00882 LDWRKU = ( LWORK-N*N ) / N 00883 LDWRKR = N 00884 END IF 00885 ITAU = IR + LDWRKR*N 00886 IWORK = ITAU + N 00887 * 00888 * Compute A=Q*R 00889 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00890 * (RWorkspace: 0) 00891 * 00892 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 00893 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00894 * 00895 * Copy R to VT, zeroing out below it 00896 * 00897 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00898 IF( N.GT.1 ) 00899 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 00900 $ VT( 2, 1 ), LDVT ) 00901 * 00902 * Generate Q in A 00903 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00904 * (RWorkspace: 0) 00905 * 00906 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00907 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00908 IE = 1 00909 ITAUQ = ITAU 00910 ITAUP = ITAUQ + N 00911 IWORK = ITAUP + N 00912 * 00913 * Bidiagonalize R in VT, copying result to WORK(IR) 00914 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 00915 * (RWorkspace: need N) 00916 * 00917 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 00918 $ WORK( ITAUQ ), WORK( ITAUP ), 00919 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00920 CALL ZLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR ) 00921 * 00922 * Generate left vectors bidiagonalizing R in WORK(IR) 00923 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00924 * (RWorkspace: 0) 00925 * 00926 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 00927 $ WORK( ITAUQ ), WORK( IWORK ), 00928 $ LWORK-IWORK+1, IERR ) 00929 * 00930 * Generate right vectors bidiagonalizing R in VT 00931 * (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB) 00932 * (RWorkspace: 0) 00933 * 00934 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00935 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00936 IRWORK = IE + N 00937 * 00938 * Perform bidiagonal QR iteration, computing left 00939 * singular vectors of R in WORK(IR) and computing right 00940 * singular vectors of R in VT 00941 * (CWorkspace: need N*N) 00942 * (RWorkspace: need BDSPAC) 00943 * 00944 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 00945 $ LDVT, WORK( IR ), LDWRKR, CDUM, 1, 00946 $ RWORK( IRWORK ), INFO ) 00947 IU = ITAUQ 00948 * 00949 * Multiply Q in A by left singular vectors of R in 00950 * WORK(IR), storing result in WORK(IU) and copying to A 00951 * (CWorkspace: need N*N+N, prefer N*N+M*N) 00952 * (RWorkspace: 0) 00953 * 00954 DO 20 I = 1, M, LDWRKU 00955 CHUNK = MIN( M-I+1, LDWRKU ) 00956 CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00957 $ LDA, WORK( IR ), LDWRKR, CZERO, 00958 $ WORK( IU ), LDWRKU ) 00959 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 00960 $ A( I, 1 ), LDA ) 00961 20 CONTINUE 00962 * 00963 ELSE 00964 * 00965 * Insufficient workspace for a fast algorithm 00966 * 00967 ITAU = 1 00968 IWORK = ITAU + N 00969 * 00970 * Compute A=Q*R 00971 * (CWorkspace: need 2*N, prefer N+N*NB) 00972 * (RWorkspace: 0) 00973 * 00974 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 00975 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00976 * 00977 * Copy R to VT, zeroing out below it 00978 * 00979 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00980 IF( N.GT.1 ) 00981 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 00982 $ VT( 2, 1 ), LDVT ) 00983 * 00984 * Generate Q in A 00985 * (CWorkspace: need 2*N, prefer N+N*NB) 00986 * (RWorkspace: 0) 00987 * 00988 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00989 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 00990 IE = 1 00991 ITAUQ = ITAU 00992 ITAUP = ITAUQ + N 00993 IWORK = ITAUP + N 00994 * 00995 * Bidiagonalize R in VT 00996 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00997 * (RWorkspace: N) 00998 * 00999 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 01000 $ WORK( ITAUQ ), WORK( ITAUP ), 01001 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01002 * 01003 * Multiply Q in A by left vectors bidiagonalizing R 01004 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01005 * (RWorkspace: 0) 01006 * 01007 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 01008 $ WORK( ITAUQ ), A, LDA, WORK( IWORK ), 01009 $ LWORK-IWORK+1, IERR ) 01010 * 01011 * Generate right vectors bidiagonalizing R in VT 01012 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01013 * (RWorkspace: 0) 01014 * 01015 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01016 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01017 IRWORK = IE + N 01018 * 01019 * Perform bidiagonal QR iteration, computing left 01020 * singular vectors of A in A and computing right 01021 * singular vectors of A in VT 01022 * (CWorkspace: 0) 01023 * (RWorkspace: need BDSPAC) 01024 * 01025 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 01026 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 01027 $ INFO ) 01028 * 01029 END IF 01030 * 01031 ELSE IF( WNTUS ) THEN 01032 * 01033 IF( WNTVN ) THEN 01034 * 01035 * Path 4 (M much larger than N, JOBU='S', JOBVT='N') 01036 * N left singular vectors to be computed in U and 01037 * no right singular vectors to be computed 01038 * 01039 IF( LWORK.GE.N*N+3*N ) THEN 01040 * 01041 * Sufficient workspace for a fast algorithm 01042 * 01043 IR = 1 01044 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01045 * 01046 * WORK(IR) is LDA by N 01047 * 01048 LDWRKR = LDA 01049 ELSE 01050 * 01051 * WORK(IR) is N by N 01052 * 01053 LDWRKR = N 01054 END IF 01055 ITAU = IR + LDWRKR*N 01056 IWORK = ITAU + N 01057 * 01058 * Compute A=Q*R 01059 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01060 * (RWorkspace: 0) 01061 * 01062 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01063 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01064 * 01065 * Copy R to WORK(IR), zeroing out below it 01066 * 01067 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), 01068 $ LDWRKR ) 01069 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01070 $ WORK( IR+1 ), LDWRKR ) 01071 * 01072 * Generate Q in A 01073 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01074 * (RWorkspace: 0) 01075 * 01076 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01077 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01078 IE = 1 01079 ITAUQ = ITAU 01080 ITAUP = ITAUQ + N 01081 IWORK = ITAUP + N 01082 * 01083 * Bidiagonalize R in WORK(IR) 01084 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01085 * (RWorkspace: need N) 01086 * 01087 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, 01088 $ RWORK( IE ), WORK( ITAUQ ), 01089 $ WORK( ITAUP ), WORK( IWORK ), 01090 $ LWORK-IWORK+1, IERR ) 01091 * 01092 * Generate left vectors bidiagonalizing R in WORK(IR) 01093 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01094 * (RWorkspace: 0) 01095 * 01096 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 01097 $ WORK( ITAUQ ), WORK( IWORK ), 01098 $ LWORK-IWORK+1, IERR ) 01099 IRWORK = IE + N 01100 * 01101 * Perform bidiagonal QR iteration, computing left 01102 * singular vectors of R in WORK(IR) 01103 * (CWorkspace: need N*N) 01104 * (RWorkspace: need BDSPAC) 01105 * 01106 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 01107 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 01108 $ RWORK( IRWORK ), INFO ) 01109 * 01110 * Multiply Q in A by left singular vectors of R in 01111 * WORK(IR), storing result in U 01112 * (CWorkspace: need N*N) 01113 * (RWorkspace: 0) 01114 * 01115 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01116 $ WORK( IR ), LDWRKR, CZERO, U, LDU ) 01117 * 01118 ELSE 01119 * 01120 * Insufficient workspace for a fast algorithm 01121 * 01122 ITAU = 1 01123 IWORK = ITAU + N 01124 * 01125 * Compute A=Q*R, copying result to U 01126 * (CWorkspace: need 2*N, prefer N+N*NB) 01127 * (RWorkspace: 0) 01128 * 01129 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01130 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01131 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01132 * 01133 * Generate Q in U 01134 * (CWorkspace: need 2*N, prefer N+N*NB) 01135 * (RWorkspace: 0) 01136 * 01137 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01138 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01139 IE = 1 01140 ITAUQ = ITAU 01141 ITAUP = ITAUQ + N 01142 IWORK = ITAUP + N 01143 * 01144 * Zero out below R in A 01145 * 01146 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01147 $ A( 2, 1 ), LDA ) 01148 * 01149 * Bidiagonalize R in A 01150 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01151 * (RWorkspace: need N) 01152 * 01153 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 01154 $ WORK( ITAUQ ), WORK( ITAUP ), 01155 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01156 * 01157 * Multiply Q in U by left vectors bidiagonalizing R 01158 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01159 * (RWorkspace: 0) 01160 * 01161 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01162 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01163 $ LWORK-IWORK+1, IERR ) 01164 IRWORK = IE + N 01165 * 01166 * Perform bidiagonal QR iteration, computing left 01167 * singular vectors of A in U 01168 * (CWorkspace: 0) 01169 * (RWorkspace: need BDSPAC) 01170 * 01171 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 01172 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 01173 $ INFO ) 01174 * 01175 END IF 01176 * 01177 ELSE IF( WNTVO ) THEN 01178 * 01179 * Path 5 (M much larger than N, JOBU='S', JOBVT='O') 01180 * N left singular vectors to be computed in U and 01181 * N right singular vectors to be overwritten on A 01182 * 01183 IF( LWORK.GE.2*N*N+3*N ) THEN 01184 * 01185 * Sufficient workspace for a fast algorithm 01186 * 01187 IU = 1 01188 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 01189 * 01190 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 01191 * 01192 LDWRKU = LDA 01193 IR = IU + LDWRKU*N 01194 LDWRKR = LDA 01195 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 01196 * 01197 * WORK(IU) is LDA by N and WORK(IR) is N by N 01198 * 01199 LDWRKU = LDA 01200 IR = IU + LDWRKU*N 01201 LDWRKR = N 01202 ELSE 01203 * 01204 * WORK(IU) is N by N and WORK(IR) is N by N 01205 * 01206 LDWRKU = N 01207 IR = IU + LDWRKU*N 01208 LDWRKR = N 01209 END IF 01210 ITAU = IR + LDWRKR*N 01211 IWORK = ITAU + N 01212 * 01213 * Compute A=Q*R 01214 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01215 * (RWorkspace: 0) 01216 * 01217 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01218 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01219 * 01220 * Copy R to WORK(IU), zeroing out below it 01221 * 01222 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 01223 $ LDWRKU ) 01224 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01225 $ WORK( IU+1 ), LDWRKU ) 01226 * 01227 * Generate Q in A 01228 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01229 * (RWorkspace: 0) 01230 * 01231 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01232 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01233 IE = 1 01234 ITAUQ = ITAU 01235 ITAUP = ITAUQ + N 01236 IWORK = ITAUP + N 01237 * 01238 * Bidiagonalize R in WORK(IU), copying result to 01239 * WORK(IR) 01240 * (CWorkspace: need 2*N*N+3*N, 01241 * prefer 2*N*N+2*N+2*N*NB) 01242 * (RWorkspace: need N) 01243 * 01244 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 01245 $ RWORK( IE ), WORK( ITAUQ ), 01246 $ WORK( ITAUP ), WORK( IWORK ), 01247 $ LWORK-IWORK+1, IERR ) 01248 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, 01249 $ WORK( IR ), LDWRKR ) 01250 * 01251 * Generate left bidiagonalizing vectors in WORK(IU) 01252 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 01253 * (RWorkspace: 0) 01254 * 01255 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01256 $ WORK( ITAUQ ), WORK( IWORK ), 01257 $ LWORK-IWORK+1, IERR ) 01258 * 01259 * Generate right bidiagonalizing vectors in WORK(IR) 01260 * (CWorkspace: need 2*N*N+3*N-1, 01261 * prefer 2*N*N+2*N+(N-1)*NB) 01262 * (RWorkspace: 0) 01263 * 01264 CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 01265 $ WORK( ITAUP ), WORK( IWORK ), 01266 $ LWORK-IWORK+1, IERR ) 01267 IRWORK = IE + N 01268 * 01269 * Perform bidiagonal QR iteration, computing left 01270 * singular vectors of R in WORK(IU) and computing 01271 * right singular vectors of R in WORK(IR) 01272 * (CWorkspace: need 2*N*N) 01273 * (RWorkspace: need BDSPAC) 01274 * 01275 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 01276 $ WORK( IR ), LDWRKR, WORK( IU ), 01277 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 01278 $ INFO ) 01279 * 01280 * Multiply Q in A by left singular vectors of R in 01281 * WORK(IU), storing result in U 01282 * (CWorkspace: need N*N) 01283 * (RWorkspace: 0) 01284 * 01285 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01286 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 01287 * 01288 * Copy right singular vectors of R to A 01289 * (CWorkspace: need N*N) 01290 * (RWorkspace: 0) 01291 * 01292 CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 01293 $ LDA ) 01294 * 01295 ELSE 01296 * 01297 * Insufficient workspace for a fast algorithm 01298 * 01299 ITAU = 1 01300 IWORK = ITAU + N 01301 * 01302 * Compute A=Q*R, copying result to U 01303 * (CWorkspace: need 2*N, prefer N+N*NB) 01304 * (RWorkspace: 0) 01305 * 01306 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01307 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01308 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01309 * 01310 * Generate Q in U 01311 * (CWorkspace: need 2*N, prefer N+N*NB) 01312 * (RWorkspace: 0) 01313 * 01314 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01315 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01316 IE = 1 01317 ITAUQ = ITAU 01318 ITAUP = ITAUQ + N 01319 IWORK = ITAUP + N 01320 * 01321 * Zero out below R in A 01322 * 01323 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01324 $ A( 2, 1 ), LDA ) 01325 * 01326 * Bidiagonalize R in A 01327 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01328 * (RWorkspace: need N) 01329 * 01330 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 01331 $ WORK( ITAUQ ), WORK( ITAUP ), 01332 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01333 * 01334 * Multiply Q in U by left vectors bidiagonalizing R 01335 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01336 * (RWorkspace: 0) 01337 * 01338 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01339 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01340 $ LWORK-IWORK+1, IERR ) 01341 * 01342 * Generate right vectors bidiagonalizing R in A 01343 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01344 * (RWorkspace: 0) 01345 * 01346 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 01347 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01348 IRWORK = IE + N 01349 * 01350 * Perform bidiagonal QR iteration, computing left 01351 * singular vectors of A in U and computing right 01352 * singular vectors of A in A 01353 * (CWorkspace: 0) 01354 * (RWorkspace: need BDSPAC) 01355 * 01356 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 01357 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 01358 $ INFO ) 01359 * 01360 END IF 01361 * 01362 ELSE IF( WNTVAS ) THEN 01363 * 01364 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' 01365 * or 'A') 01366 * N left singular vectors to be computed in U and 01367 * N right singular vectors to be computed in VT 01368 * 01369 IF( LWORK.GE.N*N+3*N ) THEN 01370 * 01371 * Sufficient workspace for a fast algorithm 01372 * 01373 IU = 1 01374 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01375 * 01376 * WORK(IU) is LDA by N 01377 * 01378 LDWRKU = LDA 01379 ELSE 01380 * 01381 * WORK(IU) is N by N 01382 * 01383 LDWRKU = N 01384 END IF 01385 ITAU = IU + LDWRKU*N 01386 IWORK = ITAU + N 01387 * 01388 * Compute A=Q*R 01389 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01390 * (RWorkspace: 0) 01391 * 01392 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01393 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01394 * 01395 * Copy R to WORK(IU), zeroing out below it 01396 * 01397 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 01398 $ LDWRKU ) 01399 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01400 $ WORK( IU+1 ), LDWRKU ) 01401 * 01402 * Generate Q in A 01403 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01404 * (RWorkspace: 0) 01405 * 01406 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 01407 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01408 IE = 1 01409 ITAUQ = ITAU 01410 ITAUP = ITAUQ + N 01411 IWORK = ITAUP + N 01412 * 01413 * Bidiagonalize R in WORK(IU), copying result to VT 01414 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01415 * (RWorkspace: need N) 01416 * 01417 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 01418 $ RWORK( IE ), WORK( ITAUQ ), 01419 $ WORK( ITAUP ), WORK( IWORK ), 01420 $ LWORK-IWORK+1, IERR ) 01421 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 01422 $ LDVT ) 01423 * 01424 * Generate left bidiagonalizing vectors in WORK(IU) 01425 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01426 * (RWorkspace: 0) 01427 * 01428 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01429 $ WORK( ITAUQ ), WORK( IWORK ), 01430 $ LWORK-IWORK+1, IERR ) 01431 * 01432 * Generate right bidiagonalizing vectors in VT 01433 * (CWorkspace: need N*N+3*N-1, 01434 * prefer N*N+2*N+(N-1)*NB) 01435 * (RWorkspace: 0) 01436 * 01437 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01438 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01439 IRWORK = IE + N 01440 * 01441 * Perform bidiagonal QR iteration, computing left 01442 * singular vectors of R in WORK(IU) and computing 01443 * right singular vectors of R in VT 01444 * (CWorkspace: need N*N) 01445 * (RWorkspace: need BDSPAC) 01446 * 01447 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 01448 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 01449 $ RWORK( IRWORK ), INFO ) 01450 * 01451 * Multiply Q in A by left singular vectors of R in 01452 * WORK(IU), storing result in U 01453 * (CWorkspace: need N*N) 01454 * (RWorkspace: 0) 01455 * 01456 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, 01457 $ WORK( IU ), LDWRKU, CZERO, U, LDU ) 01458 * 01459 ELSE 01460 * 01461 * Insufficient workspace for a fast algorithm 01462 * 01463 ITAU = 1 01464 IWORK = ITAU + N 01465 * 01466 * Compute A=Q*R, copying result to U 01467 * (CWorkspace: need 2*N, prefer N+N*NB) 01468 * (RWorkspace: 0) 01469 * 01470 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01471 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01472 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01473 * 01474 * Generate Q in U 01475 * (CWorkspace: need 2*N, prefer N+N*NB) 01476 * (RWorkspace: 0) 01477 * 01478 CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ), 01479 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01480 * 01481 * Copy R to VT, zeroing out below it 01482 * 01483 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 01484 IF( N.GT.1 ) 01485 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01486 $ VT( 2, 1 ), LDVT ) 01487 IE = 1 01488 ITAUQ = ITAU 01489 ITAUP = ITAUQ + N 01490 IWORK = ITAUP + N 01491 * 01492 * Bidiagonalize R in VT 01493 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01494 * (RWorkspace: need N) 01495 * 01496 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 01497 $ WORK( ITAUQ ), WORK( ITAUP ), 01498 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01499 * 01500 * Multiply Q in U by left bidiagonalizing vectors 01501 * in VT 01502 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01503 * (RWorkspace: 0) 01504 * 01505 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 01506 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01507 $ LWORK-IWORK+1, IERR ) 01508 * 01509 * Generate right bidiagonalizing vectors in VT 01510 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01511 * (RWorkspace: 0) 01512 * 01513 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01514 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01515 IRWORK = IE + N 01516 * 01517 * Perform bidiagonal QR iteration, computing left 01518 * singular vectors of A in U and computing right 01519 * singular vectors of A in VT 01520 * (CWorkspace: 0) 01521 * (RWorkspace: need BDSPAC) 01522 * 01523 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 01524 $ LDVT, U, LDU, CDUM, 1, 01525 $ RWORK( IRWORK ), INFO ) 01526 * 01527 END IF 01528 * 01529 END IF 01530 * 01531 ELSE IF( WNTUA ) THEN 01532 * 01533 IF( WNTVN ) THEN 01534 * 01535 * Path 7 (M much larger than N, JOBU='A', JOBVT='N') 01536 * M left singular vectors to be computed in U and 01537 * no right singular vectors to be computed 01538 * 01539 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 01540 * 01541 * Sufficient workspace for a fast algorithm 01542 * 01543 IR = 1 01544 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01545 * 01546 * WORK(IR) is LDA by N 01547 * 01548 LDWRKR = LDA 01549 ELSE 01550 * 01551 * WORK(IR) is N by N 01552 * 01553 LDWRKR = N 01554 END IF 01555 ITAU = IR + LDWRKR*N 01556 IWORK = ITAU + N 01557 * 01558 * Compute A=Q*R, copying result to U 01559 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01560 * (RWorkspace: 0) 01561 * 01562 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01563 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01564 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01565 * 01566 * Copy R to WORK(IR), zeroing out below it 01567 * 01568 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), 01569 $ LDWRKR ) 01570 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01571 $ WORK( IR+1 ), LDWRKR ) 01572 * 01573 * Generate Q in U 01574 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 01575 * (RWorkspace: 0) 01576 * 01577 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01578 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01579 IE = 1 01580 ITAUQ = ITAU 01581 ITAUP = ITAUQ + N 01582 IWORK = ITAUP + N 01583 * 01584 * Bidiagonalize R in WORK(IR) 01585 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01586 * (RWorkspace: need N) 01587 * 01588 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, 01589 $ RWORK( IE ), WORK( ITAUQ ), 01590 $ WORK( ITAUP ), WORK( IWORK ), 01591 $ LWORK-IWORK+1, IERR ) 01592 * 01593 * Generate left bidiagonalizing vectors in WORK(IR) 01594 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01595 * (RWorkspace: 0) 01596 * 01597 CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR, 01598 $ WORK( ITAUQ ), WORK( IWORK ), 01599 $ LWORK-IWORK+1, IERR ) 01600 IRWORK = IE + N 01601 * 01602 * Perform bidiagonal QR iteration, computing left 01603 * singular vectors of R in WORK(IR) 01604 * (CWorkspace: need N*N) 01605 * (RWorkspace: need BDSPAC) 01606 * 01607 CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 01608 $ 1, WORK( IR ), LDWRKR, CDUM, 1, 01609 $ RWORK( IRWORK ), INFO ) 01610 * 01611 * Multiply Q in U by left singular vectors of R in 01612 * WORK(IR), storing result in A 01613 * (CWorkspace: need N*N) 01614 * (RWorkspace: 0) 01615 * 01616 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01617 $ WORK( IR ), LDWRKR, CZERO, A, LDA ) 01618 * 01619 * Copy left singular vectors of A from A to U 01620 * 01621 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 01622 * 01623 ELSE 01624 * 01625 * Insufficient workspace for a fast algorithm 01626 * 01627 ITAU = 1 01628 IWORK = ITAU + N 01629 * 01630 * Compute A=Q*R, copying result to U 01631 * (CWorkspace: need 2*N, prefer N+N*NB) 01632 * (RWorkspace: 0) 01633 * 01634 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01635 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01636 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01637 * 01638 * Generate Q in U 01639 * (CWorkspace: need N+M, prefer N+M*NB) 01640 * (RWorkspace: 0) 01641 * 01642 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01643 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01644 IE = 1 01645 ITAUQ = ITAU 01646 ITAUP = ITAUQ + N 01647 IWORK = ITAUP + N 01648 * 01649 * Zero out below R in A 01650 * 01651 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01652 $ A( 2, 1 ), LDA ) 01653 * 01654 * Bidiagonalize R in A 01655 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01656 * (RWorkspace: need N) 01657 * 01658 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 01659 $ WORK( ITAUQ ), WORK( ITAUP ), 01660 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01661 * 01662 * Multiply Q in U by left bidiagonalizing vectors 01663 * in A 01664 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01665 * (RWorkspace: 0) 01666 * 01667 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01668 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01669 $ LWORK-IWORK+1, IERR ) 01670 IRWORK = IE + N 01671 * 01672 * Perform bidiagonal QR iteration, computing left 01673 * singular vectors of A in U 01674 * (CWorkspace: 0) 01675 * (RWorkspace: need BDSPAC) 01676 * 01677 CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 01678 $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ), 01679 $ INFO ) 01680 * 01681 END IF 01682 * 01683 ELSE IF( WNTVO ) THEN 01684 * 01685 * Path 8 (M much larger than N, JOBU='A', JOBVT='O') 01686 * M left singular vectors to be computed in U and 01687 * N right singular vectors to be overwritten on A 01688 * 01689 IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN 01690 * 01691 * Sufficient workspace for a fast algorithm 01692 * 01693 IU = 1 01694 IF( LWORK.GE.WRKBL+2*LDA*N ) THEN 01695 * 01696 * WORK(IU) is LDA by N and WORK(IR) is LDA by N 01697 * 01698 LDWRKU = LDA 01699 IR = IU + LDWRKU*N 01700 LDWRKR = LDA 01701 ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN 01702 * 01703 * WORK(IU) is LDA by N and WORK(IR) is N by N 01704 * 01705 LDWRKU = LDA 01706 IR = IU + LDWRKU*N 01707 LDWRKR = N 01708 ELSE 01709 * 01710 * WORK(IU) is N by N and WORK(IR) is N by N 01711 * 01712 LDWRKU = N 01713 IR = IU + LDWRKU*N 01714 LDWRKR = N 01715 END IF 01716 ITAU = IR + LDWRKR*N 01717 IWORK = ITAU + N 01718 * 01719 * Compute A=Q*R, copying result to U 01720 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB) 01721 * (RWorkspace: 0) 01722 * 01723 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01724 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01725 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01726 * 01727 * Generate Q in U 01728 * (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB) 01729 * (RWorkspace: 0) 01730 * 01731 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01732 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01733 * 01734 * Copy R to WORK(IU), zeroing out below it 01735 * 01736 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 01737 $ LDWRKU ) 01738 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01739 $ WORK( IU+1 ), LDWRKU ) 01740 IE = 1 01741 ITAUQ = ITAU 01742 ITAUP = ITAUQ + N 01743 IWORK = ITAUP + N 01744 * 01745 * Bidiagonalize R in WORK(IU), copying result to 01746 * WORK(IR) 01747 * (CWorkspace: need 2*N*N+3*N, 01748 * prefer 2*N*N+2*N+2*N*NB) 01749 * (RWorkspace: need N) 01750 * 01751 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 01752 $ RWORK( IE ), WORK( ITAUQ ), 01753 $ WORK( ITAUP ), WORK( IWORK ), 01754 $ LWORK-IWORK+1, IERR ) 01755 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, 01756 $ WORK( IR ), LDWRKR ) 01757 * 01758 * Generate left bidiagonalizing vectors in WORK(IU) 01759 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB) 01760 * (RWorkspace: 0) 01761 * 01762 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01763 $ WORK( ITAUQ ), WORK( IWORK ), 01764 $ LWORK-IWORK+1, IERR ) 01765 * 01766 * Generate right bidiagonalizing vectors in WORK(IR) 01767 * (CWorkspace: need 2*N*N+3*N-1, 01768 * prefer 2*N*N+2*N+(N-1)*NB) 01769 * (RWorkspace: 0) 01770 * 01771 CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR, 01772 $ WORK( ITAUP ), WORK( IWORK ), 01773 $ LWORK-IWORK+1, IERR ) 01774 IRWORK = IE + N 01775 * 01776 * Perform bidiagonal QR iteration, computing left 01777 * singular vectors of R in WORK(IU) and computing 01778 * right singular vectors of R in WORK(IR) 01779 * (CWorkspace: need 2*N*N) 01780 * (RWorkspace: need BDSPAC) 01781 * 01782 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), 01783 $ WORK( IR ), LDWRKR, WORK( IU ), 01784 $ LDWRKU, CDUM, 1, RWORK( IRWORK ), 01785 $ INFO ) 01786 * 01787 * Multiply Q in U by left singular vectors of R in 01788 * WORK(IU), storing result in A 01789 * (CWorkspace: need N*N) 01790 * (RWorkspace: 0) 01791 * 01792 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01793 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 01794 * 01795 * Copy left singular vectors of A from A to U 01796 * 01797 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 01798 * 01799 * Copy right singular vectors of R from WORK(IR) to A 01800 * 01801 CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A, 01802 $ LDA ) 01803 * 01804 ELSE 01805 * 01806 * Insufficient workspace for a fast algorithm 01807 * 01808 ITAU = 1 01809 IWORK = ITAU + N 01810 * 01811 * Compute A=Q*R, copying result to U 01812 * (CWorkspace: need 2*N, prefer N+N*NB) 01813 * (RWorkspace: 0) 01814 * 01815 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01816 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01817 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01818 * 01819 * Generate Q in U 01820 * (CWorkspace: need N+M, prefer N+M*NB) 01821 * (RWorkspace: 0) 01822 * 01823 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01824 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01825 IE = 1 01826 ITAUQ = ITAU 01827 ITAUP = ITAUQ + N 01828 IWORK = ITAUP + N 01829 * 01830 * Zero out below R in A 01831 * 01832 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01833 $ A( 2, 1 ), LDA ) 01834 * 01835 * Bidiagonalize R in A 01836 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 01837 * (RWorkspace: need N) 01838 * 01839 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), 01840 $ WORK( ITAUQ ), WORK( ITAUP ), 01841 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01842 * 01843 * Multiply Q in U by left bidiagonalizing vectors 01844 * in A 01845 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01846 * (RWorkspace: 0) 01847 * 01848 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA, 01849 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 01850 $ LWORK-IWORK+1, IERR ) 01851 * 01852 * Generate right bidiagonalizing vectors in A 01853 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 01854 * (RWorkspace: 0) 01855 * 01856 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 01857 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01858 IRWORK = IE + N 01859 * 01860 * Perform bidiagonal QR iteration, computing left 01861 * singular vectors of A in U and computing right 01862 * singular vectors of A in A 01863 * (CWorkspace: 0) 01864 * (RWorkspace: need BDSPAC) 01865 * 01866 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A, 01867 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 01868 $ INFO ) 01869 * 01870 END IF 01871 * 01872 ELSE IF( WNTVAS ) THEN 01873 * 01874 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' 01875 * or 'A') 01876 * M left singular vectors to be computed in U and 01877 * N right singular vectors to be computed in VT 01878 * 01879 IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN 01880 * 01881 * Sufficient workspace for a fast algorithm 01882 * 01883 IU = 1 01884 IF( LWORK.GE.WRKBL+LDA*N ) THEN 01885 * 01886 * WORK(IU) is LDA by N 01887 * 01888 LDWRKU = LDA 01889 ELSE 01890 * 01891 * WORK(IU) is N by N 01892 * 01893 LDWRKU = N 01894 END IF 01895 ITAU = IU + LDWRKU*N 01896 IWORK = ITAU + N 01897 * 01898 * Compute A=Q*R, copying result to U 01899 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 01900 * (RWorkspace: 0) 01901 * 01902 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01903 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01904 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01905 * 01906 * Generate Q in U 01907 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB) 01908 * (RWorkspace: 0) 01909 * 01910 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01911 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01912 * 01913 * Copy R to WORK(IU), zeroing out below it 01914 * 01915 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ), 01916 $ LDWRKU ) 01917 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 01918 $ WORK( IU+1 ), LDWRKU ) 01919 IE = 1 01920 ITAUQ = ITAU 01921 ITAUP = ITAUQ + N 01922 IWORK = ITAUP + N 01923 * 01924 * Bidiagonalize R in WORK(IU), copying result to VT 01925 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 01926 * (RWorkspace: need N) 01927 * 01928 CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S, 01929 $ RWORK( IE ), WORK( ITAUQ ), 01930 $ WORK( ITAUP ), WORK( IWORK ), 01931 $ LWORK-IWORK+1, IERR ) 01932 CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT, 01933 $ LDVT ) 01934 * 01935 * Generate left bidiagonalizing vectors in WORK(IU) 01936 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 01937 * (RWorkspace: 0) 01938 * 01939 CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU, 01940 $ WORK( ITAUQ ), WORK( IWORK ), 01941 $ LWORK-IWORK+1, IERR ) 01942 * 01943 * Generate right bidiagonalizing vectors in VT 01944 * (CWorkspace: need N*N+3*N-1, 01945 * prefer N*N+2*N+(N-1)*NB) 01946 * (RWorkspace: need 0) 01947 * 01948 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 01949 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01950 IRWORK = IE + N 01951 * 01952 * Perform bidiagonal QR iteration, computing left 01953 * singular vectors of R in WORK(IU) and computing 01954 * right singular vectors of R in VT 01955 * (CWorkspace: need N*N) 01956 * (RWorkspace: need BDSPAC) 01957 * 01958 CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT, 01959 $ LDVT, WORK( IU ), LDWRKU, CDUM, 1, 01960 $ RWORK( IRWORK ), INFO ) 01961 * 01962 * Multiply Q in U by left singular vectors of R in 01963 * WORK(IU), storing result in A 01964 * (CWorkspace: need N*N) 01965 * (RWorkspace: 0) 01966 * 01967 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, 01968 $ WORK( IU ), LDWRKU, CZERO, A, LDA ) 01969 * 01970 * Copy left singular vectors of A from A to U 01971 * 01972 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 01973 * 01974 ELSE 01975 * 01976 * Insufficient workspace for a fast algorithm 01977 * 01978 ITAU = 1 01979 IWORK = ITAU + N 01980 * 01981 * Compute A=Q*R, copying result to U 01982 * (CWorkspace: need 2*N, prefer N+N*NB) 01983 * (RWorkspace: 0) 01984 * 01985 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), 01986 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01987 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 01988 * 01989 * Generate Q in U 01990 * (CWorkspace: need N+M, prefer N+M*NB) 01991 * (RWorkspace: 0) 01992 * 01993 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 01994 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 01995 * 01996 * Copy R from A to VT, zeroing out below it 01997 * 01998 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 01999 IF( N.GT.1 ) 02000 $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, 02001 $ VT( 2, 1 ), LDVT ) 02002 IE = 1 02003 ITAUQ = ITAU 02004 ITAUP = ITAUQ + N 02005 IWORK = ITAUP + N 02006 * 02007 * Bidiagonalize R in VT 02008 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 02009 * (RWorkspace: need N) 02010 * 02011 CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ), 02012 $ WORK( ITAUQ ), WORK( ITAUP ), 02013 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02014 * 02015 * Multiply Q in U by left bidiagonalizing vectors 02016 * in VT 02017 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 02018 * (RWorkspace: 0) 02019 * 02020 CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT, 02021 $ WORK( ITAUQ ), U, LDU, WORK( IWORK ), 02022 $ LWORK-IWORK+1, IERR ) 02023 * 02024 * Generate right bidiagonalizing vectors in VT 02025 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02026 * (RWorkspace: 0) 02027 * 02028 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 02029 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02030 IRWORK = IE + N 02031 * 02032 * Perform bidiagonal QR iteration, computing left 02033 * singular vectors of A in U and computing right 02034 * singular vectors of A in VT 02035 * (CWorkspace: 0) 02036 * (RWorkspace: need BDSPAC) 02037 * 02038 CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT, 02039 $ LDVT, U, LDU, CDUM, 1, 02040 $ RWORK( IRWORK ), INFO ) 02041 * 02042 END IF 02043 * 02044 END IF 02045 * 02046 END IF 02047 * 02048 ELSE 02049 * 02050 * M .LT. MNTHR 02051 * 02052 * Path 10 (M at least N, but not much larger) 02053 * Reduce to bidiagonal form without QR decomposition 02054 * 02055 IE = 1 02056 ITAUQ = 1 02057 ITAUP = ITAUQ + N 02058 IWORK = ITAUP + N 02059 * 02060 * Bidiagonalize A 02061 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 02062 * (RWorkspace: need N) 02063 * 02064 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 02065 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 02066 $ IERR ) 02067 IF( WNTUAS ) THEN 02068 * 02069 * If left singular vectors desired in U, copy result to U 02070 * and generate left bidiagonalizing vectors in U 02071 * (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB) 02072 * (RWorkspace: 0) 02073 * 02074 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 02075 IF( WNTUS ) 02076 $ NCU = N 02077 IF( WNTUA ) 02078 $ NCU = M 02079 CALL ZUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ), 02080 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02081 END IF 02082 IF( WNTVAS ) THEN 02083 * 02084 * If right singular vectors desired in VT, copy result to 02085 * VT and generate right bidiagonalizing vectors in VT 02086 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02087 * (RWorkspace: 0) 02088 * 02089 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 02090 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 02091 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02092 END IF 02093 IF( WNTUO ) THEN 02094 * 02095 * If left singular vectors desired in A, generate left 02096 * bidiagonalizing vectors in A 02097 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 02098 * (RWorkspace: 0) 02099 * 02100 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 02101 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02102 END IF 02103 IF( WNTVO ) THEN 02104 * 02105 * If right singular vectors desired in A, generate right 02106 * bidiagonalizing vectors in A 02107 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB) 02108 * (RWorkspace: 0) 02109 * 02110 CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ), 02111 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02112 END IF 02113 IRWORK = IE + N 02114 IF( WNTUAS .OR. WNTUO ) 02115 $ NRU = M 02116 IF( WNTUN ) 02117 $ NRU = 0 02118 IF( WNTVAS .OR. WNTVO ) 02119 $ NCVT = N 02120 IF( WNTVN ) 02121 $ NCVT = 0 02122 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 02123 * 02124 * Perform bidiagonal QR iteration, if desired, computing 02125 * left singular vectors in U and computing right singular 02126 * vectors in VT 02127 * (CWorkspace: 0) 02128 * (RWorkspace: need BDSPAC) 02129 * 02130 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 02131 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 02132 $ INFO ) 02133 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 02134 * 02135 * Perform bidiagonal QR iteration, if desired, computing 02136 * left singular vectors in U and computing right singular 02137 * vectors in A 02138 * (CWorkspace: 0) 02139 * (RWorkspace: need BDSPAC) 02140 * 02141 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A, 02142 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 02143 $ INFO ) 02144 ELSE 02145 * 02146 * Perform bidiagonal QR iteration, if desired, computing 02147 * left singular vectors in A and computing right singular 02148 * vectors in VT 02149 * (CWorkspace: 0) 02150 * (RWorkspace: need BDSPAC) 02151 * 02152 CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT, 02153 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 02154 $ INFO ) 02155 END IF 02156 * 02157 END IF 02158 * 02159 ELSE 02160 * 02161 * A has more columns than rows. If A has sufficiently more 02162 * columns than rows, first reduce using the LQ decomposition (if 02163 * sufficient workspace available) 02164 * 02165 IF( N.GE.MNTHR ) THEN 02166 * 02167 IF( WNTVN ) THEN 02168 * 02169 * Path 1t(N much larger than M, JOBVT='N') 02170 * No right singular vectors to be computed 02171 * 02172 ITAU = 1 02173 IWORK = ITAU + M 02174 * 02175 * Compute A=L*Q 02176 * (CWorkspace: need 2*M, prefer M+M*NB) 02177 * (RWorkspace: 0) 02178 * 02179 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ), 02180 $ LWORK-IWORK+1, IERR ) 02181 * 02182 * Zero out above L 02183 * 02184 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 02185 $ LDA ) 02186 IE = 1 02187 ITAUQ = 1 02188 ITAUP = ITAUQ + M 02189 IWORK = ITAUP + M 02190 * 02191 * Bidiagonalize L in A 02192 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02193 * (RWorkspace: need M) 02194 * 02195 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 02196 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 02197 $ IERR ) 02198 IF( WNTUO .OR. WNTUAS ) THEN 02199 * 02200 * If left singular vectors desired, generate Q 02201 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02202 * (RWorkspace: 0) 02203 * 02204 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 02205 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02206 END IF 02207 IRWORK = IE + M 02208 NRU = 0 02209 IF( WNTUO .OR. WNTUAS ) 02210 $ NRU = M 02211 * 02212 * Perform bidiagonal QR iteration, computing left singular 02213 * vectors of A in A if desired 02214 * (CWorkspace: 0) 02215 * (RWorkspace: need BDSPAC) 02216 * 02217 CALL ZBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1, 02218 $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO ) 02219 * 02220 * If left singular vectors desired in U, copy them there 02221 * 02222 IF( WNTUAS ) 02223 $ CALL ZLACPY( 'F', M, M, A, LDA, U, LDU ) 02224 * 02225 ELSE IF( WNTVO .AND. WNTUN ) THEN 02226 * 02227 * Path 2t(N much larger than M, JOBU='N', JOBVT='O') 02228 * M right singular vectors to be overwritten on A and 02229 * no left singular vectors to be computed 02230 * 02231 IF( LWORK.GE.M*M+3*M ) THEN 02232 * 02233 * Sufficient workspace for a fast algorithm 02234 * 02235 IR = 1 02236 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 02237 * 02238 * WORK(IU) is LDA by N and WORK(IR) is LDA by M 02239 * 02240 LDWRKU = LDA 02241 CHUNK = N 02242 LDWRKR = LDA 02243 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 02244 * 02245 * WORK(IU) is LDA by N and WORK(IR) is M by M 02246 * 02247 LDWRKU = LDA 02248 CHUNK = N 02249 LDWRKR = M 02250 ELSE 02251 * 02252 * WORK(IU) is M by CHUNK and WORK(IR) is M by M 02253 * 02254 LDWRKU = M 02255 CHUNK = ( LWORK-M*M ) / M 02256 LDWRKR = M 02257 END IF 02258 ITAU = IR + LDWRKR*M 02259 IWORK = ITAU + M 02260 * 02261 * Compute A=L*Q 02262 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02263 * (RWorkspace: 0) 02264 * 02265 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02266 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02267 * 02268 * Copy L to WORK(IR) and zero out above it 02269 * 02270 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR ) 02271 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02272 $ WORK( IR+LDWRKR ), LDWRKR ) 02273 * 02274 * Generate Q in A 02275 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02276 * (RWorkspace: 0) 02277 * 02278 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02279 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02280 IE = 1 02281 ITAUQ = ITAU 02282 ITAUP = ITAUQ + M 02283 IWORK = ITAUP + M 02284 * 02285 * Bidiagonalize L in WORK(IR) 02286 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02287 * (RWorkspace: need M) 02288 * 02289 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ), 02290 $ WORK( ITAUQ ), WORK( ITAUP ), 02291 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02292 * 02293 * Generate right vectors bidiagonalizing L 02294 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 02295 * (RWorkspace: 0) 02296 * 02297 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02298 $ WORK( ITAUP ), WORK( IWORK ), 02299 $ LWORK-IWORK+1, IERR ) 02300 IRWORK = IE + M 02301 * 02302 * Perform bidiagonal QR iteration, computing right 02303 * singular vectors of L in WORK(IR) 02304 * (CWorkspace: need M*M) 02305 * (RWorkspace: need BDSPAC) 02306 * 02307 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 02308 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 02309 $ RWORK( IRWORK ), INFO ) 02310 IU = ITAUQ 02311 * 02312 * Multiply right singular vectors of L in WORK(IR) by Q 02313 * in A, storing result in WORK(IU) and copying to A 02314 * (CWorkspace: need M*M+M, prefer M*M+M*N) 02315 * (RWorkspace: 0) 02316 * 02317 DO 30 I = 1, N, CHUNK 02318 BLK = MIN( N-I+1, CHUNK ) 02319 CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 02320 $ LDWRKR, A( 1, I ), LDA, CZERO, 02321 $ WORK( IU ), LDWRKU ) 02322 CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 02323 $ A( 1, I ), LDA ) 02324 30 CONTINUE 02325 * 02326 ELSE 02327 * 02328 * Insufficient workspace for a fast algorithm 02329 * 02330 IE = 1 02331 ITAUQ = 1 02332 ITAUP = ITAUQ + M 02333 IWORK = ITAUP + M 02334 * 02335 * Bidiagonalize A 02336 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 02337 * (RWorkspace: need M) 02338 * 02339 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), 02340 $ WORK( ITAUQ ), WORK( ITAUP ), 02341 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02342 * 02343 * Generate right vectors bidiagonalizing A 02344 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02345 * (RWorkspace: 0) 02346 * 02347 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 02348 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02349 IRWORK = IE + M 02350 * 02351 * Perform bidiagonal QR iteration, computing right 02352 * singular vectors of A in A 02353 * (CWorkspace: 0) 02354 * (RWorkspace: need BDSPAC) 02355 * 02356 CALL ZBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA, 02357 $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO ) 02358 * 02359 END IF 02360 * 02361 ELSE IF( WNTVO .AND. WNTUAS ) THEN 02362 * 02363 * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O') 02364 * M right singular vectors to be overwritten on A and 02365 * M left singular vectors to be computed in U 02366 * 02367 IF( LWORK.GE.M*M+3*M ) THEN 02368 * 02369 * Sufficient workspace for a fast algorithm 02370 * 02371 IR = 1 02372 IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN 02373 * 02374 * WORK(IU) is LDA by N and WORK(IR) is LDA by M 02375 * 02376 LDWRKU = LDA 02377 CHUNK = N 02378 LDWRKR = LDA 02379 ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN 02380 * 02381 * WORK(IU) is LDA by N and WORK(IR) is M by M 02382 * 02383 LDWRKU = LDA 02384 CHUNK = N 02385 LDWRKR = M 02386 ELSE 02387 * 02388 * WORK(IU) is M by CHUNK and WORK(IR) is M by M 02389 * 02390 LDWRKU = M 02391 CHUNK = ( LWORK-M*M ) / M 02392 LDWRKR = M 02393 END IF 02394 ITAU = IR + LDWRKR*M 02395 IWORK = ITAU + M 02396 * 02397 * Compute A=L*Q 02398 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02399 * (RWorkspace: 0) 02400 * 02401 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02402 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02403 * 02404 * Copy L to U, zeroing about above it 02405 * 02406 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 02407 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 02408 $ LDU ) 02409 * 02410 * Generate Q in A 02411 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02412 * (RWorkspace: 0) 02413 * 02414 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02415 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02416 IE = 1 02417 ITAUQ = ITAU 02418 ITAUP = ITAUQ + M 02419 IWORK = ITAUP + M 02420 * 02421 * Bidiagonalize L in U, copying result to WORK(IR) 02422 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02423 * (RWorkspace: need M) 02424 * 02425 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 02426 $ WORK( ITAUQ ), WORK( ITAUP ), 02427 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02428 CALL ZLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR ) 02429 * 02430 * Generate right vectors bidiagonalizing L in WORK(IR) 02431 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB) 02432 * (RWorkspace: 0) 02433 * 02434 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02435 $ WORK( ITAUP ), WORK( IWORK ), 02436 $ LWORK-IWORK+1, IERR ) 02437 * 02438 * Generate left vectors bidiagonalizing L in U 02439 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 02440 * (RWorkspace: 0) 02441 * 02442 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02443 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02444 IRWORK = IE + M 02445 * 02446 * Perform bidiagonal QR iteration, computing left 02447 * singular vectors of L in U, and computing right 02448 * singular vectors of L in WORK(IR) 02449 * (CWorkspace: need M*M) 02450 * (RWorkspace: need BDSPAC) 02451 * 02452 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02453 $ WORK( IR ), LDWRKR, U, LDU, CDUM, 1, 02454 $ RWORK( IRWORK ), INFO ) 02455 IU = ITAUQ 02456 * 02457 * Multiply right singular vectors of L in WORK(IR) by Q 02458 * in A, storing result in WORK(IU) and copying to A 02459 * (CWorkspace: need M*M+M, prefer M*M+M*N)) 02460 * (RWorkspace: 0) 02461 * 02462 DO 40 I = 1, N, CHUNK 02463 BLK = MIN( N-I+1, CHUNK ) 02464 CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ), 02465 $ LDWRKR, A( 1, I ), LDA, CZERO, 02466 $ WORK( IU ), LDWRKU ) 02467 CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU, 02468 $ A( 1, I ), LDA ) 02469 40 CONTINUE 02470 * 02471 ELSE 02472 * 02473 * Insufficient workspace for a fast algorithm 02474 * 02475 ITAU = 1 02476 IWORK = ITAU + M 02477 * 02478 * Compute A=L*Q 02479 * (CWorkspace: need 2*M, prefer M+M*NB) 02480 * (RWorkspace: 0) 02481 * 02482 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02483 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02484 * 02485 * Copy L to U, zeroing out above it 02486 * 02487 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 02488 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ), 02489 $ LDU ) 02490 * 02491 * Generate Q in A 02492 * (CWorkspace: need 2*M, prefer M+M*NB) 02493 * (RWorkspace: 0) 02494 * 02495 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02496 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02497 IE = 1 02498 ITAUQ = ITAU 02499 ITAUP = ITAUQ + M 02500 IWORK = ITAUP + M 02501 * 02502 * Bidiagonalize L in U 02503 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02504 * (RWorkspace: need M) 02505 * 02506 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 02507 $ WORK( ITAUQ ), WORK( ITAUP ), 02508 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02509 * 02510 * Multiply right vectors bidiagonalizing L by Q in A 02511 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02512 * (RWorkspace: 0) 02513 * 02514 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 02515 $ WORK( ITAUP ), A, LDA, WORK( IWORK ), 02516 $ LWORK-IWORK+1, IERR ) 02517 * 02518 * Generate left vectors bidiagonalizing L in U 02519 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02520 * (RWorkspace: 0) 02521 * 02522 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02523 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02524 IRWORK = IE + M 02525 * 02526 * Perform bidiagonal QR iteration, computing left 02527 * singular vectors of A in U and computing right 02528 * singular vectors of A in A 02529 * (CWorkspace: 0) 02530 * (RWorkspace: need BDSPAC) 02531 * 02532 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA, 02533 $ U, LDU, CDUM, 1, RWORK( IRWORK ), INFO ) 02534 * 02535 END IF 02536 * 02537 ELSE IF( WNTVS ) THEN 02538 * 02539 IF( WNTUN ) THEN 02540 * 02541 * Path 4t(N much larger than M, JOBU='N', JOBVT='S') 02542 * M right singular vectors to be computed in VT and 02543 * no left singular vectors to be computed 02544 * 02545 IF( LWORK.GE.M*M+3*M ) THEN 02546 * 02547 * Sufficient workspace for a fast algorithm 02548 * 02549 IR = 1 02550 IF( LWORK.GE.WRKBL+LDA*M ) THEN 02551 * 02552 * WORK(IR) is LDA by M 02553 * 02554 LDWRKR = LDA 02555 ELSE 02556 * 02557 * WORK(IR) is M by M 02558 * 02559 LDWRKR = M 02560 END IF 02561 ITAU = IR + LDWRKR*M 02562 IWORK = ITAU + M 02563 * 02564 * Compute A=L*Q 02565 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02566 * (RWorkspace: 0) 02567 * 02568 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02569 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02570 * 02571 * Copy L to WORK(IR), zeroing out above it 02572 * 02573 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), 02574 $ LDWRKR ) 02575 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02576 $ WORK( IR+LDWRKR ), LDWRKR ) 02577 * 02578 * Generate Q in A 02579 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02580 * (RWorkspace: 0) 02581 * 02582 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02583 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02584 IE = 1 02585 ITAUQ = ITAU 02586 ITAUP = ITAUQ + M 02587 IWORK = ITAUP + M 02588 * 02589 * Bidiagonalize L in WORK(IR) 02590 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02591 * (RWorkspace: need M) 02592 * 02593 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, 02594 $ RWORK( IE ), WORK( ITAUQ ), 02595 $ WORK( ITAUP ), WORK( IWORK ), 02596 $ LWORK-IWORK+1, IERR ) 02597 * 02598 * Generate right vectors bidiagonalizing L in 02599 * WORK(IR) 02600 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 02601 * (RWorkspace: 0) 02602 * 02603 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 02604 $ WORK( ITAUP ), WORK( IWORK ), 02605 $ LWORK-IWORK+1, IERR ) 02606 IRWORK = IE + M 02607 * 02608 * Perform bidiagonal QR iteration, computing right 02609 * singular vectors of L in WORK(IR) 02610 * (CWorkspace: need M*M) 02611 * (RWorkspace: need BDSPAC) 02612 * 02613 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 02614 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 02615 $ RWORK( IRWORK ), INFO ) 02616 * 02617 * Multiply right singular vectors of L in WORK(IR) by 02618 * Q in A, storing result in VT 02619 * (CWorkspace: need M*M) 02620 * (RWorkspace: 0) 02621 * 02622 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 02623 $ LDWRKR, A, LDA, CZERO, VT, LDVT ) 02624 * 02625 ELSE 02626 * 02627 * Insufficient workspace for a fast algorithm 02628 * 02629 ITAU = 1 02630 IWORK = ITAU + M 02631 * 02632 * Compute A=L*Q 02633 * (CWorkspace: need 2*M, prefer M+M*NB) 02634 * (RWorkspace: 0) 02635 * 02636 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02637 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02638 * 02639 * Copy result to VT 02640 * 02641 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02642 * 02643 * Generate Q in VT 02644 * (CWorkspace: need 2*M, prefer M+M*NB) 02645 * (RWorkspace: 0) 02646 * 02647 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02648 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02649 IE = 1 02650 ITAUQ = ITAU 02651 ITAUP = ITAUQ + M 02652 IWORK = ITAUP + M 02653 * 02654 * Zero out above L in A 02655 * 02656 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02657 $ A( 1, 2 ), LDA ) 02658 * 02659 * Bidiagonalize L in A 02660 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02661 * (RWorkspace: need M) 02662 * 02663 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 02664 $ WORK( ITAUQ ), WORK( ITAUP ), 02665 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02666 * 02667 * Multiply right vectors bidiagonalizing L by Q in VT 02668 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02669 * (RWorkspace: 0) 02670 * 02671 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 02672 $ WORK( ITAUP ), VT, LDVT, 02673 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02674 IRWORK = IE + M 02675 * 02676 * Perform bidiagonal QR iteration, computing right 02677 * singular vectors of A in VT 02678 * (CWorkspace: 0) 02679 * (RWorkspace: need BDSPAC) 02680 * 02681 CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 02682 $ LDVT, CDUM, 1, CDUM, 1, 02683 $ RWORK( IRWORK ), INFO ) 02684 * 02685 END IF 02686 * 02687 ELSE IF( WNTUO ) THEN 02688 * 02689 * Path 5t(N much larger than M, JOBU='O', JOBVT='S') 02690 * M right singular vectors to be computed in VT and 02691 * M left singular vectors to be overwritten on A 02692 * 02693 IF( LWORK.GE.2*M*M+3*M ) THEN 02694 * 02695 * Sufficient workspace for a fast algorithm 02696 * 02697 IU = 1 02698 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 02699 * 02700 * WORK(IU) is LDA by M and WORK(IR) is LDA by M 02701 * 02702 LDWRKU = LDA 02703 IR = IU + LDWRKU*M 02704 LDWRKR = LDA 02705 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 02706 * 02707 * WORK(IU) is LDA by M and WORK(IR) is M by M 02708 * 02709 LDWRKU = LDA 02710 IR = IU + LDWRKU*M 02711 LDWRKR = M 02712 ELSE 02713 * 02714 * WORK(IU) is M by M and WORK(IR) is M by M 02715 * 02716 LDWRKU = M 02717 IR = IU + LDWRKU*M 02718 LDWRKR = M 02719 END IF 02720 ITAU = IR + LDWRKR*M 02721 IWORK = ITAU + M 02722 * 02723 * Compute A=L*Q 02724 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 02725 * (RWorkspace: 0) 02726 * 02727 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02728 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02729 * 02730 * Copy L to WORK(IU), zeroing out below it 02731 * 02732 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 02733 $ LDWRKU ) 02734 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02735 $ WORK( IU+LDWRKU ), LDWRKU ) 02736 * 02737 * Generate Q in A 02738 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 02739 * (RWorkspace: 0) 02740 * 02741 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02742 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02743 IE = 1 02744 ITAUQ = ITAU 02745 ITAUP = ITAUQ + M 02746 IWORK = ITAUP + M 02747 * 02748 * Bidiagonalize L in WORK(IU), copying result to 02749 * WORK(IR) 02750 * (CWorkspace: need 2*M*M+3*M, 02751 * prefer 2*M*M+2*M+2*M*NB) 02752 * (RWorkspace: need M) 02753 * 02754 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 02755 $ RWORK( IE ), WORK( ITAUQ ), 02756 $ WORK( ITAUP ), WORK( IWORK ), 02757 $ LWORK-IWORK+1, IERR ) 02758 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, 02759 $ WORK( IR ), LDWRKR ) 02760 * 02761 * Generate right bidiagonalizing vectors in WORK(IU) 02762 * (CWorkspace: need 2*M*M+3*M-1, 02763 * prefer 2*M*M+2*M+(M-1)*NB) 02764 * (RWorkspace: 0) 02765 * 02766 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 02767 $ WORK( ITAUP ), WORK( IWORK ), 02768 $ LWORK-IWORK+1, IERR ) 02769 * 02770 * Generate left bidiagonalizing vectors in WORK(IR) 02771 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 02772 * (RWorkspace: 0) 02773 * 02774 CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 02775 $ WORK( ITAUQ ), WORK( IWORK ), 02776 $ LWORK-IWORK+1, IERR ) 02777 IRWORK = IE + M 02778 * 02779 * Perform bidiagonal QR iteration, computing left 02780 * singular vectors of L in WORK(IR) and computing 02781 * right singular vectors of L in WORK(IU) 02782 * (CWorkspace: need 2*M*M) 02783 * (RWorkspace: need BDSPAC) 02784 * 02785 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02786 $ WORK( IU ), LDWRKU, WORK( IR ), 02787 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 02788 $ INFO ) 02789 * 02790 * Multiply right singular vectors of L in WORK(IU) by 02791 * Q in A, storing result in VT 02792 * (CWorkspace: need M*M) 02793 * (RWorkspace: 0) 02794 * 02795 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 02796 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 02797 * 02798 * Copy left singular vectors of L to A 02799 * (CWorkspace: need M*M) 02800 * (RWorkspace: 0) 02801 * 02802 CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 02803 $ LDA ) 02804 * 02805 ELSE 02806 * 02807 * Insufficient workspace for a fast algorithm 02808 * 02809 ITAU = 1 02810 IWORK = ITAU + M 02811 * 02812 * Compute A=L*Q, copying result to VT 02813 * (CWorkspace: need 2*M, prefer M+M*NB) 02814 * (RWorkspace: 0) 02815 * 02816 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02817 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02818 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02819 * 02820 * Generate Q in VT 02821 * (CWorkspace: need 2*M, prefer M+M*NB) 02822 * (RWorkspace: 0) 02823 * 02824 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02825 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02826 IE = 1 02827 ITAUQ = ITAU 02828 ITAUP = ITAUQ + M 02829 IWORK = ITAUP + M 02830 * 02831 * Zero out above L in A 02832 * 02833 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02834 $ A( 1, 2 ), LDA ) 02835 * 02836 * Bidiagonalize L in A 02837 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 02838 * (RWorkspace: need M) 02839 * 02840 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 02841 $ WORK( ITAUQ ), WORK( ITAUP ), 02842 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02843 * 02844 * Multiply right vectors bidiagonalizing L by Q in VT 02845 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 02846 * (RWorkspace: 0) 02847 * 02848 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 02849 $ WORK( ITAUP ), VT, LDVT, 02850 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02851 * 02852 * Generate left bidiagonalizing vectors of L in A 02853 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 02854 * (RWorkspace: 0) 02855 * 02856 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 02857 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02858 IRWORK = IE + M 02859 * 02860 * Perform bidiagonal QR iteration, computing left 02861 * singular vectors of A in A and computing right 02862 * singular vectors of A in VT 02863 * (CWorkspace: 0) 02864 * (RWorkspace: need BDSPAC) 02865 * 02866 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 02867 $ LDVT, A, LDA, CDUM, 1, 02868 $ RWORK( IRWORK ), INFO ) 02869 * 02870 END IF 02871 * 02872 ELSE IF( WNTUAS ) THEN 02873 * 02874 * Path 6t(N much larger than M, JOBU='S' or 'A', 02875 * JOBVT='S') 02876 * M right singular vectors to be computed in VT and 02877 * M left singular vectors to be computed in U 02878 * 02879 IF( LWORK.GE.M*M+3*M ) THEN 02880 * 02881 * Sufficient workspace for a fast algorithm 02882 * 02883 IU = 1 02884 IF( LWORK.GE.WRKBL+LDA*M ) THEN 02885 * 02886 * WORK(IU) is LDA by N 02887 * 02888 LDWRKU = LDA 02889 ELSE 02890 * 02891 * WORK(IU) is LDA by M 02892 * 02893 LDWRKU = M 02894 END IF 02895 ITAU = IU + LDWRKU*M 02896 IWORK = ITAU + M 02897 * 02898 * Compute A=L*Q 02899 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02900 * (RWorkspace: 0) 02901 * 02902 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02903 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02904 * 02905 * Copy L to WORK(IU), zeroing out above it 02906 * 02907 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 02908 $ LDWRKU ) 02909 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02910 $ WORK( IU+LDWRKU ), LDWRKU ) 02911 * 02912 * Generate Q in A 02913 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 02914 * (RWorkspace: 0) 02915 * 02916 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 02917 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02918 IE = 1 02919 ITAUQ = ITAU 02920 ITAUP = ITAUQ + M 02921 IWORK = ITAUP + M 02922 * 02923 * Bidiagonalize L in WORK(IU), copying result to U 02924 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 02925 * (RWorkspace: need M) 02926 * 02927 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 02928 $ RWORK( IE ), WORK( ITAUQ ), 02929 $ WORK( ITAUP ), WORK( IWORK ), 02930 $ LWORK-IWORK+1, IERR ) 02931 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 02932 $ LDU ) 02933 * 02934 * Generate right bidiagonalizing vectors in WORK(IU) 02935 * (CWorkspace: need M*M+3*M-1, 02936 * prefer M*M+2*M+(M-1)*NB) 02937 * (RWorkspace: 0) 02938 * 02939 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 02940 $ WORK( ITAUP ), WORK( IWORK ), 02941 $ LWORK-IWORK+1, IERR ) 02942 * 02943 * Generate left bidiagonalizing vectors in U 02944 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 02945 * (RWorkspace: 0) 02946 * 02947 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 02948 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02949 IRWORK = IE + M 02950 * 02951 * Perform bidiagonal QR iteration, computing left 02952 * singular vectors of L in U and computing right 02953 * singular vectors of L in WORK(IU) 02954 * (CWorkspace: need M*M) 02955 * (RWorkspace: need BDSPAC) 02956 * 02957 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 02958 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 02959 $ RWORK( IRWORK ), INFO ) 02960 * 02961 * Multiply right singular vectors of L in WORK(IU) by 02962 * Q in A, storing result in VT 02963 * (CWorkspace: need M*M) 02964 * (RWorkspace: 0) 02965 * 02966 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 02967 $ LDWRKU, A, LDA, CZERO, VT, LDVT ) 02968 * 02969 ELSE 02970 * 02971 * Insufficient workspace for a fast algorithm 02972 * 02973 ITAU = 1 02974 IWORK = ITAU + M 02975 * 02976 * Compute A=L*Q, copying result to VT 02977 * (CWorkspace: need 2*M, prefer M+M*NB) 02978 * (RWorkspace: 0) 02979 * 02980 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 02981 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02982 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 02983 * 02984 * Generate Q in VT 02985 * (CWorkspace: need 2*M, prefer M+M*NB) 02986 * (RWorkspace: 0) 02987 * 02988 CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ), 02989 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 02990 * 02991 * Copy L to U, zeroing out above it 02992 * 02993 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 02994 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 02995 $ U( 1, 2 ), LDU ) 02996 IE = 1 02997 ITAUQ = ITAU 02998 ITAUP = ITAUQ + M 02999 IWORK = ITAUP + M 03000 * 03001 * Bidiagonalize L in U 03002 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03003 * (RWorkspace: need M) 03004 * 03005 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 03006 $ WORK( ITAUQ ), WORK( ITAUP ), 03007 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03008 * 03009 * Multiply right bidiagonalizing vectors in U by Q 03010 * in VT 03011 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03012 * (RWorkspace: 0) 03013 * 03014 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 03015 $ WORK( ITAUP ), VT, LDVT, 03016 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03017 * 03018 * Generate left bidiagonalizing vectors in U 03019 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03020 * (RWorkspace: 0) 03021 * 03022 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03023 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03024 IRWORK = IE + M 03025 * 03026 * Perform bidiagonal QR iteration, computing left 03027 * singular vectors of A in U and computing right 03028 * singular vectors of A in VT 03029 * (CWorkspace: 0) 03030 * (RWorkspace: need BDSPAC) 03031 * 03032 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03033 $ LDVT, U, LDU, CDUM, 1, 03034 $ RWORK( IRWORK ), INFO ) 03035 * 03036 END IF 03037 * 03038 END IF 03039 * 03040 ELSE IF( WNTVA ) THEN 03041 * 03042 IF( WNTUN ) THEN 03043 * 03044 * Path 7t(N much larger than M, JOBU='N', JOBVT='A') 03045 * N right singular vectors to be computed in VT and 03046 * no left singular vectors to be computed 03047 * 03048 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 03049 * 03050 * Sufficient workspace for a fast algorithm 03051 * 03052 IR = 1 03053 IF( LWORK.GE.WRKBL+LDA*M ) THEN 03054 * 03055 * WORK(IR) is LDA by M 03056 * 03057 LDWRKR = LDA 03058 ELSE 03059 * 03060 * WORK(IR) is M by M 03061 * 03062 LDWRKR = M 03063 END IF 03064 ITAU = IR + LDWRKR*M 03065 IWORK = ITAU + M 03066 * 03067 * Compute A=L*Q, copying result to VT 03068 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 03069 * (RWorkspace: 0) 03070 * 03071 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03072 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03073 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03074 * 03075 * Copy L to WORK(IR), zeroing out above it 03076 * 03077 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), 03078 $ LDWRKR ) 03079 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03080 $ WORK( IR+LDWRKR ), LDWRKR ) 03081 * 03082 * Generate Q in VT 03083 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 03084 * (RWorkspace: 0) 03085 * 03086 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03087 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03088 IE = 1 03089 ITAUQ = ITAU 03090 ITAUP = ITAUQ + M 03091 IWORK = ITAUP + M 03092 * 03093 * Bidiagonalize L in WORK(IR) 03094 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 03095 * (RWorkspace: need M) 03096 * 03097 CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, 03098 $ RWORK( IE ), WORK( ITAUQ ), 03099 $ WORK( ITAUP ), WORK( IWORK ), 03100 $ LWORK-IWORK+1, IERR ) 03101 * 03102 * Generate right bidiagonalizing vectors in WORK(IR) 03103 * (CWorkspace: need M*M+3*M-1, 03104 * prefer M*M+2*M+(M-1)*NB) 03105 * (RWorkspace: 0) 03106 * 03107 CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR, 03108 $ WORK( ITAUP ), WORK( IWORK ), 03109 $ LWORK-IWORK+1, IERR ) 03110 IRWORK = IE + M 03111 * 03112 * Perform bidiagonal QR iteration, computing right 03113 * singular vectors of L in WORK(IR) 03114 * (CWorkspace: need M*M) 03115 * (RWorkspace: need BDSPAC) 03116 * 03117 CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ), 03118 $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1, 03119 $ RWORK( IRWORK ), INFO ) 03120 * 03121 * Multiply right singular vectors of L in WORK(IR) by 03122 * Q in VT, storing result in A 03123 * (CWorkspace: need M*M) 03124 * (RWorkspace: 0) 03125 * 03126 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ), 03127 $ LDWRKR, VT, LDVT, CZERO, A, LDA ) 03128 * 03129 * Copy right singular vectors of A from A to VT 03130 * 03131 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03132 * 03133 ELSE 03134 * 03135 * Insufficient workspace for a fast algorithm 03136 * 03137 ITAU = 1 03138 IWORK = ITAU + M 03139 * 03140 * Compute A=L*Q, copying result to VT 03141 * (CWorkspace: need 2*M, prefer M+M*NB) 03142 * (RWorkspace: 0) 03143 * 03144 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03145 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03146 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03147 * 03148 * Generate Q in VT 03149 * (CWorkspace: need M+N, prefer M+N*NB) 03150 * (RWorkspace: 0) 03151 * 03152 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03153 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03154 IE = 1 03155 ITAUQ = ITAU 03156 ITAUP = ITAUQ + M 03157 IWORK = ITAUP + M 03158 * 03159 * Zero out above L in A 03160 * 03161 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03162 $ A( 1, 2 ), LDA ) 03163 * 03164 * Bidiagonalize L in A 03165 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03166 * (RWorkspace: need M) 03167 * 03168 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 03169 $ WORK( ITAUQ ), WORK( ITAUP ), 03170 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03171 * 03172 * Multiply right bidiagonalizing vectors in A by Q 03173 * in VT 03174 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03175 * (RWorkspace: 0) 03176 * 03177 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 03178 $ WORK( ITAUP ), VT, LDVT, 03179 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03180 IRWORK = IE + M 03181 * 03182 * Perform bidiagonal QR iteration, computing right 03183 * singular vectors of A in VT 03184 * (CWorkspace: 0) 03185 * (RWorkspace: need BDSPAC) 03186 * 03187 CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT, 03188 $ LDVT, CDUM, 1, CDUM, 1, 03189 $ RWORK( IRWORK ), INFO ) 03190 * 03191 END IF 03192 * 03193 ELSE IF( WNTUO ) THEN 03194 * 03195 * Path 8t(N much larger than M, JOBU='O', JOBVT='A') 03196 * N right singular vectors to be computed in VT and 03197 * M left singular vectors to be overwritten on A 03198 * 03199 IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN 03200 * 03201 * Sufficient workspace for a fast algorithm 03202 * 03203 IU = 1 03204 IF( LWORK.GE.WRKBL+2*LDA*M ) THEN 03205 * 03206 * WORK(IU) is LDA by M and WORK(IR) is LDA by M 03207 * 03208 LDWRKU = LDA 03209 IR = IU + LDWRKU*M 03210 LDWRKR = LDA 03211 ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN 03212 * 03213 * WORK(IU) is LDA by M and WORK(IR) is M by M 03214 * 03215 LDWRKU = LDA 03216 IR = IU + LDWRKU*M 03217 LDWRKR = M 03218 ELSE 03219 * 03220 * WORK(IU) is M by M and WORK(IR) is M by M 03221 * 03222 LDWRKU = M 03223 IR = IU + LDWRKU*M 03224 LDWRKR = M 03225 END IF 03226 ITAU = IR + LDWRKR*M 03227 IWORK = ITAU + M 03228 * 03229 * Compute A=L*Q, copying result to VT 03230 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB) 03231 * (RWorkspace: 0) 03232 * 03233 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03234 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03235 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03236 * 03237 * Generate Q in VT 03238 * (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB) 03239 * (RWorkspace: 0) 03240 * 03241 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03242 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03243 * 03244 * Copy L to WORK(IU), zeroing out above it 03245 * 03246 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 03247 $ LDWRKU ) 03248 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03249 $ WORK( IU+LDWRKU ), LDWRKU ) 03250 IE = 1 03251 ITAUQ = ITAU 03252 ITAUP = ITAUQ + M 03253 IWORK = ITAUP + M 03254 * 03255 * Bidiagonalize L in WORK(IU), copying result to 03256 * WORK(IR) 03257 * (CWorkspace: need 2*M*M+3*M, 03258 * prefer 2*M*M+2*M+2*M*NB) 03259 * (RWorkspace: need M) 03260 * 03261 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 03262 $ RWORK( IE ), WORK( ITAUQ ), 03263 $ WORK( ITAUP ), WORK( IWORK ), 03264 $ LWORK-IWORK+1, IERR ) 03265 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, 03266 $ WORK( IR ), LDWRKR ) 03267 * 03268 * Generate right bidiagonalizing vectors in WORK(IU) 03269 * (CWorkspace: need 2*M*M+3*M-1, 03270 * prefer 2*M*M+2*M+(M-1)*NB) 03271 * (RWorkspace: 0) 03272 * 03273 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 03274 $ WORK( ITAUP ), WORK( IWORK ), 03275 $ LWORK-IWORK+1, IERR ) 03276 * 03277 * Generate left bidiagonalizing vectors in WORK(IR) 03278 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB) 03279 * (RWorkspace: 0) 03280 * 03281 CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR, 03282 $ WORK( ITAUQ ), WORK( IWORK ), 03283 $ LWORK-IWORK+1, IERR ) 03284 IRWORK = IE + M 03285 * 03286 * Perform bidiagonal QR iteration, computing left 03287 * singular vectors of L in WORK(IR) and computing 03288 * right singular vectors of L in WORK(IU) 03289 * (CWorkspace: need 2*M*M) 03290 * (RWorkspace: need BDSPAC) 03291 * 03292 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 03293 $ WORK( IU ), LDWRKU, WORK( IR ), 03294 $ LDWRKR, CDUM, 1, RWORK( IRWORK ), 03295 $ INFO ) 03296 * 03297 * Multiply right singular vectors of L in WORK(IU) by 03298 * Q in VT, storing result in A 03299 * (CWorkspace: need M*M) 03300 * (RWorkspace: 0) 03301 * 03302 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 03303 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 03304 * 03305 * Copy right singular vectors of A from A to VT 03306 * 03307 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03308 * 03309 * Copy left singular vectors of A from WORK(IR) to A 03310 * 03311 CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A, 03312 $ LDA ) 03313 * 03314 ELSE 03315 * 03316 * Insufficient workspace for a fast algorithm 03317 * 03318 ITAU = 1 03319 IWORK = ITAU + M 03320 * 03321 * Compute A=L*Q, copying result to VT 03322 * (CWorkspace: need 2*M, prefer M+M*NB) 03323 * (RWorkspace: 0) 03324 * 03325 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03326 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03327 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03328 * 03329 * Generate Q in VT 03330 * (CWorkspace: need M+N, prefer M+N*NB) 03331 * (RWorkspace: 0) 03332 * 03333 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03334 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03335 IE = 1 03336 ITAUQ = ITAU 03337 ITAUP = ITAUQ + M 03338 IWORK = ITAUP + M 03339 * 03340 * Zero out above L in A 03341 * 03342 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03343 $ A( 1, 2 ), LDA ) 03344 * 03345 * Bidiagonalize L in A 03346 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03347 * (RWorkspace: need M) 03348 * 03349 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), 03350 $ WORK( ITAUQ ), WORK( ITAUP ), 03351 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03352 * 03353 * Multiply right bidiagonalizing vectors in A by Q 03354 * in VT 03355 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03356 * (RWorkspace: 0) 03357 * 03358 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA, 03359 $ WORK( ITAUP ), VT, LDVT, 03360 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03361 * 03362 * Generate left bidiagonalizing vectors in A 03363 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03364 * (RWorkspace: 0) 03365 * 03366 CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ), 03367 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03368 IRWORK = IE + M 03369 * 03370 * Perform bidiagonal QR iteration, computing left 03371 * singular vectors of A in A and computing right 03372 * singular vectors of A in VT 03373 * (CWorkspace: 0) 03374 * (RWorkspace: need BDSPAC) 03375 * 03376 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03377 $ LDVT, A, LDA, CDUM, 1, 03378 $ RWORK( IRWORK ), INFO ) 03379 * 03380 END IF 03381 * 03382 ELSE IF( WNTUAS ) THEN 03383 * 03384 * Path 9t(N much larger than M, JOBU='S' or 'A', 03385 * JOBVT='A') 03386 * N right singular vectors to be computed in VT and 03387 * M left singular vectors to be computed in U 03388 * 03389 IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN 03390 * 03391 * Sufficient workspace for a fast algorithm 03392 * 03393 IU = 1 03394 IF( LWORK.GE.WRKBL+LDA*M ) THEN 03395 * 03396 * WORK(IU) is LDA by M 03397 * 03398 LDWRKU = LDA 03399 ELSE 03400 * 03401 * WORK(IU) is M by M 03402 * 03403 LDWRKU = M 03404 END IF 03405 ITAU = IU + LDWRKU*M 03406 IWORK = ITAU + M 03407 * 03408 * Compute A=L*Q, copying result to VT 03409 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 03410 * (RWorkspace: 0) 03411 * 03412 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03413 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03414 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03415 * 03416 * Generate Q in VT 03417 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB) 03418 * (RWorkspace: 0) 03419 * 03420 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03421 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03422 * 03423 * Copy L to WORK(IU), zeroing out above it 03424 * 03425 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ), 03426 $ LDWRKU ) 03427 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03428 $ WORK( IU+LDWRKU ), LDWRKU ) 03429 IE = 1 03430 ITAUQ = ITAU 03431 ITAUP = ITAUQ + M 03432 IWORK = ITAUP + M 03433 * 03434 * Bidiagonalize L in WORK(IU), copying result to U 03435 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 03436 * (RWorkspace: need M) 03437 * 03438 CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S, 03439 $ RWORK( IE ), WORK( ITAUQ ), 03440 $ WORK( ITAUP ), WORK( IWORK ), 03441 $ LWORK-IWORK+1, IERR ) 03442 CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U, 03443 $ LDU ) 03444 * 03445 * Generate right bidiagonalizing vectors in WORK(IU) 03446 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB) 03447 * (RWorkspace: 0) 03448 * 03449 CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU, 03450 $ WORK( ITAUP ), WORK( IWORK ), 03451 $ LWORK-IWORK+1, IERR ) 03452 * 03453 * Generate left bidiagonalizing vectors in U 03454 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 03455 * (RWorkspace: 0) 03456 * 03457 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03458 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03459 IRWORK = IE + M 03460 * 03461 * Perform bidiagonal QR iteration, computing left 03462 * singular vectors of L in U and computing right 03463 * singular vectors of L in WORK(IU) 03464 * (CWorkspace: need M*M) 03465 * (RWorkspace: need BDSPAC) 03466 * 03467 CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ), 03468 $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1, 03469 $ RWORK( IRWORK ), INFO ) 03470 * 03471 * Multiply right singular vectors of L in WORK(IU) by 03472 * Q in VT, storing result in A 03473 * (CWorkspace: need M*M) 03474 * (RWorkspace: 0) 03475 * 03476 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ), 03477 $ LDWRKU, VT, LDVT, CZERO, A, LDA ) 03478 * 03479 * Copy right singular vectors of A from A to VT 03480 * 03481 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 03482 * 03483 ELSE 03484 * 03485 * Insufficient workspace for a fast algorithm 03486 * 03487 ITAU = 1 03488 IWORK = ITAU + M 03489 * 03490 * Compute A=L*Q, copying result to VT 03491 * (CWorkspace: need 2*M, prefer M+M*NB) 03492 * (RWorkspace: 0) 03493 * 03494 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), 03495 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03496 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03497 * 03498 * Generate Q in VT 03499 * (CWorkspace: need M+N, prefer M+N*NB) 03500 * (RWorkspace: 0) 03501 * 03502 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 03503 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03504 * 03505 * Copy L to U, zeroing out above it 03506 * 03507 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 03508 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 03509 $ U( 1, 2 ), LDU ) 03510 IE = 1 03511 ITAUQ = ITAU 03512 ITAUP = ITAUQ + M 03513 IWORK = ITAUP + M 03514 * 03515 * Bidiagonalize L in U 03516 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 03517 * (RWorkspace: need M) 03518 * 03519 CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ), 03520 $ WORK( ITAUQ ), WORK( ITAUP ), 03521 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03522 * 03523 * Multiply right bidiagonalizing vectors in U by Q 03524 * in VT 03525 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 03526 * (RWorkspace: 0) 03527 * 03528 CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU, 03529 $ WORK( ITAUP ), VT, LDVT, 03530 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03531 * 03532 * Generate left bidiagonalizing vectors in U 03533 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03534 * (RWorkspace: 0) 03535 * 03536 CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ), 03537 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03538 IRWORK = IE + M 03539 * 03540 * Perform bidiagonal QR iteration, computing left 03541 * singular vectors of A in U and computing right 03542 * singular vectors of A in VT 03543 * (CWorkspace: 0) 03544 * (RWorkspace: need BDSPAC) 03545 * 03546 CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT, 03547 $ LDVT, U, LDU, CDUM, 1, 03548 $ RWORK( IRWORK ), INFO ) 03549 * 03550 END IF 03551 * 03552 END IF 03553 * 03554 END IF 03555 * 03556 ELSE 03557 * 03558 * N .LT. MNTHR 03559 * 03560 * Path 10t(N greater than M, but not much larger) 03561 * Reduce to bidiagonal form without LQ decomposition 03562 * 03563 IE = 1 03564 ITAUQ = 1 03565 ITAUP = ITAUQ + M 03566 IWORK = ITAUP + M 03567 * 03568 * Bidiagonalize A 03569 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 03570 * (RWorkspace: M) 03571 * 03572 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 03573 $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1, 03574 $ IERR ) 03575 IF( WNTUAS ) THEN 03576 * 03577 * If left singular vectors desired in U, copy result to U 03578 * and generate left bidiagonalizing vectors in U 03579 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 03580 * (RWorkspace: 0) 03581 * 03582 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 03583 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 03584 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03585 END IF 03586 IF( WNTVAS ) THEN 03587 * 03588 * If right singular vectors desired in VT, copy result to 03589 * VT and generate right bidiagonalizing vectors in VT 03590 * (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB) 03591 * (RWorkspace: 0) 03592 * 03593 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 03594 IF( WNTVA ) 03595 $ NRVT = N 03596 IF( WNTVS ) 03597 $ NRVT = M 03598 CALL ZUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ), 03599 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03600 END IF 03601 IF( WNTUO ) THEN 03602 * 03603 * If left singular vectors desired in A, generate left 03604 * bidiagonalizing vectors in A 03605 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB) 03606 * (RWorkspace: 0) 03607 * 03608 CALL ZUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ), 03609 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03610 END IF 03611 IF( WNTVO ) THEN 03612 * 03613 * If right singular vectors desired in A, generate right 03614 * bidiagonalizing vectors in A 03615 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 03616 * (RWorkspace: 0) 03617 * 03618 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 03619 $ WORK( IWORK ), LWORK-IWORK+1, IERR ) 03620 END IF 03621 IRWORK = IE + M 03622 IF( WNTUAS .OR. WNTUO ) 03623 $ NRU = M 03624 IF( WNTUN ) 03625 $ NRU = 0 03626 IF( WNTVAS .OR. WNTVO ) 03627 $ NCVT = N 03628 IF( WNTVN ) 03629 $ NCVT = 0 03630 IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN 03631 * 03632 * Perform bidiagonal QR iteration, if desired, computing 03633 * left singular vectors in U and computing right singular 03634 * vectors in VT 03635 * (CWorkspace: 0) 03636 * (RWorkspace: need BDSPAC) 03637 * 03638 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 03639 $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ), 03640 $ INFO ) 03641 ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN 03642 * 03643 * Perform bidiagonal QR iteration, if desired, computing 03644 * left singular vectors in U and computing right singular 03645 * vectors in A 03646 * (CWorkspace: 0) 03647 * (RWorkspace: need BDSPAC) 03648 * 03649 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A, 03650 $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ), 03651 $ INFO ) 03652 ELSE 03653 * 03654 * Perform bidiagonal QR iteration, if desired, computing 03655 * left singular vectors in A and computing right singular 03656 * vectors in VT 03657 * (CWorkspace: 0) 03658 * (RWorkspace: need BDSPAC) 03659 * 03660 CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT, 03661 $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ), 03662 $ INFO ) 03663 END IF 03664 * 03665 END IF 03666 * 03667 END IF 03668 * 03669 * Undo scaling if necessary 03670 * 03671 IF( ISCL.EQ.1 ) THEN 03672 IF( ANRM.GT.BIGNUM ) 03673 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 03674 $ IERR ) 03675 IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) 03676 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, 03677 $ RWORK( IE ), MINMN, IERR ) 03678 IF( ANRM.LT.SMLNUM ) 03679 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 03680 $ IERR ) 03681 IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) 03682 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, 03683 $ RWORK( IE ), MINMN, IERR ) 03684 END IF 03685 * 03686 * Return optimal workspace in WORK(1) 03687 * 03688 WORK( 1 ) = MAXWRK 03689 * 03690 RETURN 03691 * 03692 * End of ZGESVD 03693 * 03694 END