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