 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, 00002 $ LWORK, RWORK, IWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2.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 * June 2010 00008 * 8-15-00: Improve consistency of WS calculations (eca) 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBZ 00012 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IWORK( * ) 00016 DOUBLE PRECISION RWORK( * ), S( * ) 00017 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ), 00018 $ WORK( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZGESDD computes the singular value decomposition (SVD) of a complex 00025 * M-by-N matrix A, optionally computing the left and/or right singular 00026 * vectors, by using divide-and-conquer method. The SVD is written 00027 * 00028 * A = U * SIGMA * conjugate-transpose(V) 00029 * 00030 * where SIGMA is an M-by-N matrix which is zero except for its 00031 * min(m,n) diagonal elements, U is an M-by-M unitary matrix, and 00032 * V is an N-by-N unitary matrix. The diagonal elements of SIGMA 00033 * are the singular values of A; they are real and non-negative, and 00034 * are returned in descending order. The first min(m,n) columns of 00035 * U and V are the left and right singular vectors of A. 00036 * 00037 * Note that the routine returns VT = V**H, not V. 00038 * 00039 * The divide and conquer algorithm makes very mild assumptions about 00040 * floating point arithmetic. It will work on machines with a guard 00041 * digit in add/subtract, or on those binary machines without guard 00042 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00043 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00044 * without guard digits, but we know of none. 00045 * 00046 * Arguments 00047 * ========= 00048 * 00049 * JOBZ (input) CHARACTER*1 00050 * Specifies options for computing all or part of the matrix U: 00051 * = 'A': all M columns of U and all N rows of V**H are 00052 * returned in the arrays U and VT; 00053 * = 'S': the first min(M,N) columns of U and the first 00054 * min(M,N) rows of V**H are returned in the arrays U 00055 * and VT; 00056 * = 'O': If M >= N, the first N columns of U are overwritten 00057 * in the array A and all rows of V**H are returned in 00058 * the array VT; 00059 * otherwise, all columns of U are returned in the 00060 * array U and the first M rows of V**H are overwritten 00061 * in the array A; 00062 * = 'N': no columns of U or rows of V**H are computed. 00063 * 00064 * M (input) INTEGER 00065 * The number of rows of the input matrix A. M >= 0. 00066 * 00067 * N (input) INTEGER 00068 * The number of columns of the input matrix A. N >= 0. 00069 * 00070 * A (input/output) COMPLEX*16 array, dimension (LDA,N) 00071 * On entry, the M-by-N matrix A. 00072 * On exit, 00073 * if JOBZ = 'O', A is overwritten with the first N columns 00074 * of U (the left singular vectors, stored 00075 * columnwise) if M >= N; 00076 * A is overwritten with the first M rows 00077 * of V**H (the right singular vectors, stored 00078 * rowwise) otherwise. 00079 * if JOBZ .ne. 'O', the contents of A are destroyed. 00080 * 00081 * LDA (input) INTEGER 00082 * The leading dimension of the array A. LDA >= max(1,M). 00083 * 00084 * S (output) DOUBLE PRECISION array, dimension (min(M,N)) 00085 * The singular values of A, sorted so that S(i) >= S(i+1). 00086 * 00087 * U (output) COMPLEX*16 array, dimension (LDU,UCOL) 00088 * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N; 00089 * UCOL = min(M,N) if JOBZ = 'S'. 00090 * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M 00091 * unitary matrix U; 00092 * if JOBZ = 'S', U contains the first min(M,N) columns of U 00093 * (the left singular vectors, stored columnwise); 00094 * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced. 00095 * 00096 * LDU (input) INTEGER 00097 * The leading dimension of the array U. LDU >= 1; if 00098 * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M. 00099 * 00100 * VT (output) COMPLEX*16 array, dimension (LDVT,N) 00101 * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the 00102 * N-by-N unitary matrix V**H; 00103 * if JOBZ = 'S', VT contains the first min(M,N) rows of 00104 * V**H (the right singular vectors, stored rowwise); 00105 * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced. 00106 * 00107 * LDVT (input) INTEGER 00108 * The leading dimension of the array VT. LDVT >= 1; if 00109 * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N; 00110 * if JOBZ = 'S', LDVT >= min(M,N). 00111 * 00112 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00113 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00114 * 00115 * LWORK (input) INTEGER 00116 * The dimension of the array WORK. LWORK >= 1. 00117 * if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N). 00118 * if JOBZ = 'O', 00119 * LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 00120 * if JOBZ = 'S' or 'A', 00121 * LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N). 00122 * For good performance, LWORK should generally be larger. 00123 * 00124 * If LWORK = -1, a workspace query is assumed. The optimal 00125 * size for the WORK array is calculated and stored in WORK(1), 00126 * and no other work except argument checking is performed. 00127 * 00128 * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 00129 * If JOBZ = 'N', LRWORK >= 5*min(M,N). 00130 * Otherwise, 00131 * LRWORK >= min(M,N)*max(5*min(M,N)+7,2*max(M,N)+2*min(M,N)+1) 00132 * 00133 * IWORK (workspace) INTEGER array, dimension (8*min(M,N)) 00134 * 00135 * INFO (output) INTEGER 00136 * = 0: successful exit. 00137 * < 0: if INFO = -i, the i-th argument had an illegal value. 00138 * > 0: The updating process of DBDSDC did not converge. 00139 * 00140 * Further Details 00141 * =============== 00142 * 00143 * Based on contributions by 00144 * Ming Gu and Huan Ren, Computer Science Division, University of 00145 * California at Berkeley, USA 00146 * 00147 * ===================================================================== 00148 * 00149 * .. Parameters .. 00150 INTEGER LQUERV 00151 PARAMETER ( LQUERV = -1 ) 00152 COMPLEX*16 CZERO, CONE 00153 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00154 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00155 DOUBLE PRECISION ZERO, ONE 00156 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00157 * .. 00158 * .. Local Scalars .. 00159 LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS 00160 INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT, 00161 $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT, 00162 $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK, 00163 $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL 00164 DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM 00165 * .. 00166 * .. Local Arrays .. 00167 INTEGER IDUM( 1 ) 00168 DOUBLE PRECISION DUM( 1 ) 00169 * .. 00170 * .. External Subroutines .. 00171 EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM, 00172 $ ZGEQRF, ZLACP2, ZLACPY, ZLACRM, ZLARCM, ZLASCL, 00173 $ ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR 00174 * .. 00175 * .. External Functions .. 00176 LOGICAL LSAME 00177 INTEGER ILAENV 00178 DOUBLE PRECISION DLAMCH, ZLANGE 00179 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE 00180 * .. 00181 * .. Intrinsic Functions .. 00182 INTRINSIC INT, MAX, MIN, SQRT 00183 * .. 00184 * .. Executable Statements .. 00185 * 00186 * Test the input arguments 00187 * 00188 INFO = 0 00189 MINMN = MIN( M, N ) 00190 MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 ) 00191 MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 ) 00192 WNTQA = LSAME( JOBZ, 'A' ) 00193 WNTQS = LSAME( JOBZ, 'S' ) 00194 WNTQAS = WNTQA .OR. WNTQS 00195 WNTQO = LSAME( JOBZ, 'O' ) 00196 WNTQN = LSAME( JOBZ, 'N' ) 00197 MINWRK = 1 00198 MAXWRK = 1 00199 * 00200 IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN 00201 INFO = -1 00202 ELSE IF( M.LT.0 ) THEN 00203 INFO = -2 00204 ELSE IF( N.LT.0 ) THEN 00205 INFO = -3 00206 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00207 INFO = -5 00208 ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR. 00209 $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN 00210 INFO = -8 00211 ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR. 00212 $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR. 00213 $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN 00214 INFO = -10 00215 END IF 00216 * 00217 * Compute workspace 00218 * (Note: Comments in the code beginning "Workspace:" describe the 00219 * minimal amount of workspace needed at that point in the code, 00220 * as well as the preferred amount for good performance. 00221 * CWorkspace refers to complex workspace, and RWorkspace to 00222 * real workspace. NB refers to the optimal block size for the 00223 * immediately following subroutine, as returned by ILAENV.) 00224 * 00225 IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN 00226 IF( M.GE.N ) THEN 00227 * 00228 * There is no complex work space needed for bidiagonal SVD 00229 * The real work space needed for bidiagonal SVD is BDSPAC 00230 * for computing singular values and singular vectors; BDSPAN 00231 * for computing singular values only. 00232 * BDSPAC = 5*N*N + 7*N 00233 * BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8)) 00234 * 00235 IF( M.GE.MNTHR1 ) THEN 00236 IF( WNTQN ) THEN 00237 * 00238 * Path 1 (M much larger than N, JOBZ='N') 00239 * 00240 MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, 00241 $ -1 ) 00242 MAXWRK = MAX( MAXWRK, 2*N+2*N* 00243 $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) 00244 MINWRK = 3*N 00245 ELSE IF( WNTQO ) THEN 00246 * 00247 * Path 2 (M much larger than N, JOBZ='O') 00248 * 00249 WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) 00250 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, 00251 $ N, N, -1 ) ) 00252 WRKBL = MAX( WRKBL, 2*N+2*N* 00253 $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) 00254 WRKBL = MAX( WRKBL, 2*N+N* 00255 $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) 00256 WRKBL = MAX( WRKBL, 2*N+N* 00257 $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) 00258 MAXWRK = M*N + N*N + WRKBL 00259 MINWRK = 2*N*N + 3*N 00260 ELSE IF( WNTQS ) THEN 00261 * 00262 * Path 3 (M much larger than N, JOBZ='S') 00263 * 00264 WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) 00265 WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M, 00266 $ N, N, -1 ) ) 00267 WRKBL = MAX( WRKBL, 2*N+2*N* 00268 $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) 00269 WRKBL = MAX( WRKBL, 2*N+N* 00270 $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) 00271 WRKBL = MAX( WRKBL, 2*N+N* 00272 $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) 00273 MAXWRK = N*N + WRKBL 00274 MINWRK = N*N + 3*N 00275 ELSE IF( WNTQA ) THEN 00276 * 00277 * Path 4 (M much larger than N, JOBZ='A') 00278 * 00279 WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) 00280 WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M, 00281 $ M, N, -1 ) ) 00282 WRKBL = MAX( WRKBL, 2*N+2*N* 00283 $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) ) 00284 WRKBL = MAX( WRKBL, 2*N+N* 00285 $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) ) 00286 WRKBL = MAX( WRKBL, 2*N+N* 00287 $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) 00288 MAXWRK = N*N + WRKBL 00289 MINWRK = N*N + 2*N + M 00290 END IF 00291 ELSE IF( M.GE.MNTHR2 ) THEN 00292 * 00293 * Path 5 (M much larger than N, but not as much as MNTHR1) 00294 * 00295 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, 00296 $ -1, -1 ) 00297 MINWRK = 2*N + M 00298 IF( WNTQO ) THEN 00299 MAXWRK = MAX( MAXWRK, 2*N+N* 00300 $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) 00301 MAXWRK = MAX( MAXWRK, 2*N+N* 00302 $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) ) 00303 MAXWRK = MAXWRK + M*N 00304 MINWRK = MINWRK + N*N 00305 ELSE IF( WNTQS ) THEN 00306 MAXWRK = MAX( MAXWRK, 2*N+N* 00307 $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) 00308 MAXWRK = MAX( MAXWRK, 2*N+N* 00309 $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) ) 00310 ELSE IF( WNTQA ) THEN 00311 MAXWRK = MAX( MAXWRK, 2*N+N* 00312 $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) ) 00313 MAXWRK = MAX( MAXWRK, 2*N+M* 00314 $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) 00315 END IF 00316 ELSE 00317 * 00318 * Path 6 (M at least N, but not much larger) 00319 * 00320 MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, 00321 $ -1, -1 ) 00322 MINWRK = 2*N + M 00323 IF( WNTQO ) THEN 00324 MAXWRK = MAX( MAXWRK, 2*N+N* 00325 $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) 00326 MAXWRK = MAX( MAXWRK, 2*N+N* 00327 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) ) 00328 MAXWRK = MAXWRK + M*N 00329 MINWRK = MINWRK + N*N 00330 ELSE IF( WNTQS ) THEN 00331 MAXWRK = MAX( MAXWRK, 2*N+N* 00332 $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) ) 00333 MAXWRK = MAX( MAXWRK, 2*N+N* 00334 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) ) 00335 ELSE IF( WNTQA ) THEN 00336 MAXWRK = MAX( MAXWRK, 2*N+N* 00337 $ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, N, -1 ) ) 00338 MAXWRK = MAX( MAXWRK, 2*N+M* 00339 $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) 00340 END IF 00341 END IF 00342 ELSE 00343 * 00344 * There is no complex work space needed for bidiagonal SVD 00345 * The real work space needed for bidiagonal SVD is BDSPAC 00346 * for computing singular values and singular vectors; BDSPAN 00347 * for computing singular values only. 00348 * BDSPAC = 5*M*M + 7*M 00349 * BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8)) 00350 * 00351 IF( N.GE.MNTHR1 ) THEN 00352 IF( WNTQN ) THEN 00353 * 00354 * Path 1t (N much larger than M, JOBZ='N') 00355 * 00356 MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, 00357 $ -1 ) 00358 MAXWRK = MAX( MAXWRK, 2*M+2*M* 00359 $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00360 MINWRK = 3*M 00361 ELSE IF( WNTQO ) THEN 00362 * 00363 * Path 2t (N much larger than M, JOBZ='O') 00364 * 00365 WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) 00366 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, 00367 $ N, M, -1 ) ) 00368 WRKBL = MAX( WRKBL, 2*M+2*M* 00369 $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00370 WRKBL = MAX( WRKBL, 2*M+M* 00371 $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) 00372 WRKBL = MAX( WRKBL, 2*M+M* 00373 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) 00374 MAXWRK = M*N + M*M + WRKBL 00375 MINWRK = 2*M*M + 3*M 00376 ELSE IF( WNTQS ) THEN 00377 * 00378 * Path 3t (N much larger than M, JOBZ='S') 00379 * 00380 WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) 00381 WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M, 00382 $ N, M, -1 ) ) 00383 WRKBL = MAX( WRKBL, 2*M+2*M* 00384 $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00385 WRKBL = MAX( WRKBL, 2*M+M* 00386 $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) 00387 WRKBL = MAX( WRKBL, 2*M+M* 00388 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) 00389 MAXWRK = M*M + WRKBL 00390 MINWRK = M*M + 3*M 00391 ELSE IF( WNTQA ) THEN 00392 * 00393 * Path 4t (N much larger than M, JOBZ='A') 00394 * 00395 WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) 00396 WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N, 00397 $ N, M, -1 ) ) 00398 WRKBL = MAX( WRKBL, 2*M+2*M* 00399 $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) ) 00400 WRKBL = MAX( WRKBL, 2*M+M* 00401 $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) ) 00402 WRKBL = MAX( WRKBL, 2*M+M* 00403 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) ) 00404 MAXWRK = M*M + WRKBL 00405 MINWRK = M*M + 2*M + N 00406 END IF 00407 ELSE IF( N.GE.MNTHR2 ) THEN 00408 * 00409 * Path 5t (N much larger than M, but not as much as MNTHR1) 00410 * 00411 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, 00412 $ -1, -1 ) 00413 MINWRK = 2*M + N 00414 IF( WNTQO ) THEN 00415 MAXWRK = MAX( MAXWRK, 2*M+M* 00416 $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) ) 00417 MAXWRK = MAX( MAXWRK, 2*M+M* 00418 $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) 00419 MAXWRK = MAXWRK + M*N 00420 MINWRK = MINWRK + M*M 00421 ELSE IF( WNTQS ) THEN 00422 MAXWRK = MAX( MAXWRK, 2*M+M* 00423 $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) ) 00424 MAXWRK = MAX( MAXWRK, 2*M+M* 00425 $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) 00426 ELSE IF( WNTQA ) THEN 00427 MAXWRK = MAX( MAXWRK, 2*M+N* 00428 $ ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) ) 00429 MAXWRK = MAX( MAXWRK, 2*M+M* 00430 $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) ) 00431 END IF 00432 ELSE 00433 * 00434 * Path 6t (N greater than M, but not much larger) 00435 * 00436 MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N, 00437 $ -1, -1 ) 00438 MINWRK = 2*M + N 00439 IF( WNTQO ) THEN 00440 MAXWRK = MAX( MAXWRK, 2*M+M* 00441 $ ILAENV( 1, 'ZUNMBR', 'PRC', M, N, M, -1 ) ) 00442 MAXWRK = MAX( MAXWRK, 2*M+M* 00443 $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 ) ) 00444 MAXWRK = MAXWRK + M*N 00445 MINWRK = MINWRK + M*M 00446 ELSE IF( WNTQS ) THEN 00447 MAXWRK = MAX( MAXWRK, 2*M+M* 00448 $ ILAENV( 1, 'ZUNGBR', 'PRC', M, N, M, -1 ) ) 00449 MAXWRK = MAX( MAXWRK, 2*M+M* 00450 $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) 00451 ELSE IF( WNTQA ) THEN 00452 MAXWRK = MAX( MAXWRK, 2*M+N* 00453 $ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, M, -1 ) ) 00454 MAXWRK = MAX( MAXWRK, 2*M+M* 00455 $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) ) 00456 END IF 00457 END IF 00458 END IF 00459 MAXWRK = MAX( MAXWRK, MINWRK ) 00460 END IF 00461 IF( INFO.EQ.0 ) THEN 00462 WORK( 1 ) = MAXWRK 00463 IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV ) 00464 $ INFO = -13 00465 END IF 00466 * 00467 * Quick returns 00468 * 00469 IF( INFO.NE.0 ) THEN 00470 CALL XERBLA( 'ZGESDD', -INFO ) 00471 RETURN 00472 END IF 00473 IF( LWORK.EQ.LQUERV ) 00474 $ RETURN 00475 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00476 RETURN 00477 END IF 00478 * 00479 * Get machine constants 00480 * 00481 EPS = DLAMCH( 'P' ) 00482 SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS 00483 BIGNUM = ONE / SMLNUM 00484 * 00485 * Scale A if max element outside range [SMLNUM,BIGNUM] 00486 * 00487 ANRM = ZLANGE( 'M', M, N, A, LDA, DUM ) 00488 ISCL = 0 00489 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00490 ISCL = 1 00491 CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR ) 00492 ELSE IF( ANRM.GT.BIGNUM ) THEN 00493 ISCL = 1 00494 CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR ) 00495 END IF 00496 * 00497 IF( M.GE.N ) THEN 00498 * 00499 * A has at least as many rows as columns. If A has sufficiently 00500 * more rows than columns, first reduce using the QR 00501 * decomposition (if sufficient workspace available) 00502 * 00503 IF( M.GE.MNTHR1 ) THEN 00504 * 00505 IF( WNTQN ) THEN 00506 * 00507 * Path 1 (M much larger than N, JOBZ='N') 00508 * No singular vectors to be computed 00509 * 00510 ITAU = 1 00511 NWORK = ITAU + N 00512 * 00513 * Compute A=Q*R 00514 * (CWorkspace: need 2*N, prefer N+N*NB) 00515 * (RWorkspace: need 0) 00516 * 00517 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00518 $ LWORK-NWORK+1, IERR ) 00519 * 00520 * Zero out below R 00521 * 00522 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00523 $ LDA ) 00524 IE = 1 00525 ITAUQ = 1 00526 ITAUP = ITAUQ + N 00527 NWORK = ITAUP + N 00528 * 00529 * Bidiagonalize R in A 00530 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00531 * (RWorkspace: need N) 00532 * 00533 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00534 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00535 $ IERR ) 00536 NRWORK = IE + N 00537 * 00538 * Perform bidiagonal SVD, compute singular values only 00539 * (CWorkspace: 0) 00540 * (RWorkspace: need BDSPAN) 00541 * 00542 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 00543 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 00544 * 00545 ELSE IF( WNTQO ) THEN 00546 * 00547 * Path 2 (M much larger than N, JOBZ='O') 00548 * N left singular vectors to be overwritten on A and 00549 * N right singular vectors to be computed in VT 00550 * 00551 IU = 1 00552 * 00553 * WORK(IU) is N by N 00554 * 00555 LDWRKU = N 00556 IR = IU + LDWRKU*N 00557 IF( LWORK.GE.M*N+N*N+3*N ) THEN 00558 * 00559 * WORK(IR) is M by N 00560 * 00561 LDWRKR = M 00562 ELSE 00563 LDWRKR = ( LWORK-N*N-3*N ) / N 00564 END IF 00565 ITAU = IR + LDWRKR*N 00566 NWORK = ITAU + N 00567 * 00568 * Compute A=Q*R 00569 * (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB) 00570 * (RWorkspace: 0) 00571 * 00572 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00573 $ LWORK-NWORK+1, IERR ) 00574 * 00575 * Copy R to WORK( IR ), zeroing out below it 00576 * 00577 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00578 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 00579 $ LDWRKR ) 00580 * 00581 * Generate Q in A 00582 * (CWorkspace: need 2*N, prefer N+N*NB) 00583 * (RWorkspace: 0) 00584 * 00585 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00586 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00587 IE = 1 00588 ITAUQ = ITAU 00589 ITAUP = ITAUQ + N 00590 NWORK = ITAUP + N 00591 * 00592 * Bidiagonalize R in WORK(IR) 00593 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB) 00594 * (RWorkspace: need N) 00595 * 00596 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 00597 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00598 $ LWORK-NWORK+1, IERR ) 00599 * 00600 * Perform bidiagonal SVD, computing left singular vectors 00601 * of R in WORK(IRU) and computing right singular vectors 00602 * of R in WORK(IRVT) 00603 * (CWorkspace: need 0) 00604 * (RWorkspace: need BDSPAC) 00605 * 00606 IRU = IE + N 00607 IRVT = IRU + N*N 00608 NRWORK = IRVT + N*N 00609 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00610 $ N, RWORK( IRVT ), N, DUM, IDUM, 00611 $ RWORK( NRWORK ), IWORK, INFO ) 00612 * 00613 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 00614 * Overwrite WORK(IU) by the left singular vectors of R 00615 * (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB) 00616 * (RWorkspace: 0) 00617 * 00618 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 00619 $ LDWRKU ) 00620 CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00621 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00622 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00623 * 00624 * Copy real matrix RWORK(IRVT) to complex matrix VT 00625 * Overwrite VT by the right singular vectors of R 00626 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 00627 * (RWorkspace: 0) 00628 * 00629 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00630 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 00631 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00632 $ LWORK-NWORK+1, IERR ) 00633 * 00634 * Multiply Q in A by left singular vectors of R in 00635 * WORK(IU), storing result in WORK(IR) and copying to A 00636 * (CWorkspace: need 2*N*N, prefer N*N+M*N) 00637 * (RWorkspace: 0) 00638 * 00639 DO 10 I = 1, M, LDWRKR 00640 CHUNK = MIN( M-I+1, LDWRKR ) 00641 CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ), 00642 $ LDA, WORK( IU ), LDWRKU, CZERO, 00643 $ WORK( IR ), LDWRKR ) 00644 CALL ZLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR, 00645 $ A( I, 1 ), LDA ) 00646 10 CONTINUE 00647 * 00648 ELSE IF( WNTQS ) THEN 00649 * 00650 * Path 3 (M much larger than N, JOBZ='S') 00651 * N left singular vectors to be computed in U and 00652 * N right singular vectors to be computed in VT 00653 * 00654 IR = 1 00655 * 00656 * WORK(IR) is N by N 00657 * 00658 LDWRKR = N 00659 ITAU = IR + LDWRKR*N 00660 NWORK = ITAU + N 00661 * 00662 * Compute A=Q*R 00663 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB) 00664 * (RWorkspace: 0) 00665 * 00666 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00667 $ LWORK-NWORK+1, IERR ) 00668 * 00669 * Copy R to WORK(IR), zeroing out below it 00670 * 00671 CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR ) 00672 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ), 00673 $ LDWRKR ) 00674 * 00675 * Generate Q in A 00676 * (CWorkspace: need 2*N, prefer N+N*NB) 00677 * (RWorkspace: 0) 00678 * 00679 CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ), 00680 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00681 IE = 1 00682 ITAUQ = ITAU 00683 ITAUP = ITAUQ + N 00684 NWORK = ITAUP + N 00685 * 00686 * Bidiagonalize R in WORK(IR) 00687 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB) 00688 * (RWorkspace: need N) 00689 * 00690 CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ), 00691 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 00692 $ LWORK-NWORK+1, IERR ) 00693 * 00694 * Perform bidiagonal SVD, computing left singular vectors 00695 * of bidiagonal matrix in RWORK(IRU) and computing right 00696 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00697 * (CWorkspace: need 0) 00698 * (RWorkspace: need BDSPAC) 00699 * 00700 IRU = IE + N 00701 IRVT = IRU + N*N 00702 NRWORK = IRVT + N*N 00703 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00704 $ N, RWORK( IRVT ), N, DUM, IDUM, 00705 $ RWORK( NRWORK ), IWORK, INFO ) 00706 * 00707 * Copy real matrix RWORK(IRU) to complex matrix U 00708 * Overwrite U by left singular vectors of R 00709 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00710 * (RWorkspace: 0) 00711 * 00712 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 00713 CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR, 00714 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 00715 $ LWORK-NWORK+1, IERR ) 00716 * 00717 * Copy real matrix RWORK(IRVT) to complex matrix VT 00718 * Overwrite VT by right singular vectors of R 00719 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00720 * (RWorkspace: 0) 00721 * 00722 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00723 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR, 00724 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00725 $ LWORK-NWORK+1, IERR ) 00726 * 00727 * Multiply Q in A by left singular vectors of R in 00728 * WORK(IR), storing result in U 00729 * (CWorkspace: need N*N) 00730 * (RWorkspace: 0) 00731 * 00732 CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR ) 00733 CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ), 00734 $ LDWRKR, CZERO, U, LDU ) 00735 * 00736 ELSE IF( WNTQA ) THEN 00737 * 00738 * Path 4 (M much larger than N, JOBZ='A') 00739 * M left singular vectors to be computed in U and 00740 * N right singular vectors to be computed in VT 00741 * 00742 IU = 1 00743 * 00744 * WORK(IU) is N by N 00745 * 00746 LDWRKU = N 00747 ITAU = IU + LDWRKU*N 00748 NWORK = ITAU + N 00749 * 00750 * Compute A=Q*R, copying result to U 00751 * (CWorkspace: need 2*N, prefer N+N*NB) 00752 * (RWorkspace: 0) 00753 * 00754 CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 00755 $ LWORK-NWORK+1, IERR ) 00756 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 00757 * 00758 * Generate Q in U 00759 * (CWorkspace: need N+M, prefer N+M*NB) 00760 * (RWorkspace: 0) 00761 * 00762 CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ), 00763 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00764 * 00765 * Produce R in A, zeroing out below it 00766 * 00767 CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ), 00768 $ LDA ) 00769 IE = 1 00770 ITAUQ = ITAU 00771 ITAUP = ITAUQ + N 00772 NWORK = ITAUP + N 00773 * 00774 * Bidiagonalize R in A 00775 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB) 00776 * (RWorkspace: need N) 00777 * 00778 CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00779 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00780 $ IERR ) 00781 IRU = IE + N 00782 IRVT = IRU + N*N 00783 NRWORK = IRVT + N*N 00784 * 00785 * Perform bidiagonal SVD, computing left singular vectors 00786 * of bidiagonal matrix in RWORK(IRU) and computing right 00787 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00788 * (CWorkspace: need 0) 00789 * (RWorkspace: need BDSPAC) 00790 * 00791 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00792 $ N, RWORK( IRVT ), N, DUM, IDUM, 00793 $ RWORK( NRWORK ), IWORK, INFO ) 00794 * 00795 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 00796 * Overwrite WORK(IU) by left singular vectors of R 00797 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB) 00798 * (RWorkspace: 0) 00799 * 00800 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 00801 $ LDWRKU ) 00802 CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA, 00803 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 00804 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00805 * 00806 * Copy real matrix RWORK(IRVT) to complex matrix VT 00807 * Overwrite VT by right singular vectors of R 00808 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 00809 * (RWorkspace: 0) 00810 * 00811 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 00812 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 00813 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 00814 $ LWORK-NWORK+1, IERR ) 00815 * 00816 * Multiply Q in U by left singular vectors of R in 00817 * WORK(IU), storing result in A 00818 * (CWorkspace: need N*N) 00819 * (RWorkspace: 0) 00820 * 00821 CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ), 00822 $ LDWRKU, CZERO, A, LDA ) 00823 * 00824 * Copy left singular vectors of A from A to U 00825 * 00826 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 00827 * 00828 END IF 00829 * 00830 ELSE IF( M.GE.MNTHR2 ) THEN 00831 * 00832 * MNTHR2 <= M < MNTHR1 00833 * 00834 * Path 5 (M much larger than N, but not as much as MNTHR1) 00835 * Reduce to bidiagonal form without QR decomposition, use 00836 * ZUNGBR and matrix multiplication to compute singular vectors 00837 * 00838 IE = 1 00839 NRWORK = IE + N 00840 ITAUQ = 1 00841 ITAUP = ITAUQ + N 00842 NWORK = ITAUP + N 00843 * 00844 * Bidiagonalize A 00845 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 00846 * (RWorkspace: need N) 00847 * 00848 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 00849 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 00850 $ IERR ) 00851 IF( WNTQN ) THEN 00852 * 00853 * Compute singular values only 00854 * (Cworkspace: 0) 00855 * (Rworkspace: need BDSPAN) 00856 * 00857 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 00858 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 00859 ELSE IF( WNTQO ) THEN 00860 IU = NWORK 00861 IRU = NRWORK 00862 IRVT = IRU + N*N 00863 NRWORK = IRVT + N*N 00864 * 00865 * Copy A to VT, generate P**H 00866 * (Cworkspace: need 2*N, prefer N+N*NB) 00867 * (Rworkspace: 0) 00868 * 00869 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00870 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00871 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00872 * 00873 * Generate Q in A 00874 * (CWorkspace: need 2*N, prefer N+N*NB) 00875 * (RWorkspace: 0) 00876 * 00877 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 00878 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00879 * 00880 IF( LWORK.GE.M*N+3*N ) THEN 00881 * 00882 * WORK( IU ) is M by N 00883 * 00884 LDWRKU = M 00885 ELSE 00886 * 00887 * WORK(IU) is LDWRKU by N 00888 * 00889 LDWRKU = ( LWORK-3*N ) / N 00890 END IF 00891 NWORK = IU + LDWRKU*N 00892 * 00893 * Perform bidiagonal SVD, computing left singular vectors 00894 * of bidiagonal matrix in RWORK(IRU) and computing right 00895 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00896 * (CWorkspace: need 0) 00897 * (RWorkspace: need BDSPAC) 00898 * 00899 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00900 $ N, RWORK( IRVT ), N, DUM, IDUM, 00901 $ RWORK( NRWORK ), IWORK, INFO ) 00902 * 00903 * Multiply real matrix RWORK(IRVT) by P**H in VT, 00904 * storing the result in WORK(IU), copying to VT 00905 * (Cworkspace: need 0) 00906 * (Rworkspace: need 3*N*N) 00907 * 00908 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, 00909 $ WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 00910 CALL ZLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT ) 00911 * 00912 * Multiply Q in A by real matrix RWORK(IRU), storing the 00913 * result in WORK(IU), copying to A 00914 * (CWorkspace: need N*N, prefer M*N) 00915 * (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 00916 * 00917 NRWORK = IRVT 00918 DO 20 I = 1, M, LDWRKU 00919 CHUNK = MIN( M-I+1, LDWRKU ) 00920 CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ), 00921 $ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) ) 00922 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 00923 $ A( I, 1 ), LDA ) 00924 20 CONTINUE 00925 * 00926 ELSE IF( WNTQS ) THEN 00927 * 00928 * Copy A to VT, generate P**H 00929 * (Cworkspace: need 2*N, prefer N+N*NB) 00930 * (Rworkspace: 0) 00931 * 00932 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00933 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00934 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00935 * 00936 * Copy A to U, generate Q 00937 * (Cworkspace: need 2*N, prefer N+N*NB) 00938 * (Rworkspace: 0) 00939 * 00940 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 00941 CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ), 00942 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00943 * 00944 * Perform bidiagonal SVD, computing left singular vectors 00945 * of bidiagonal matrix in RWORK(IRU) and computing right 00946 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00947 * (CWorkspace: need 0) 00948 * (RWorkspace: need BDSPAC) 00949 * 00950 IRU = NRWORK 00951 IRVT = IRU + N*N 00952 NRWORK = IRVT + N*N 00953 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 00954 $ N, RWORK( IRVT ), N, DUM, IDUM, 00955 $ RWORK( NRWORK ), IWORK, INFO ) 00956 * 00957 * Multiply real matrix RWORK(IRVT) by P**H in VT, 00958 * storing the result in A, copying to VT 00959 * (Cworkspace: need 0) 00960 * (Rworkspace: need 3*N*N) 00961 * 00962 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 00963 $ RWORK( NRWORK ) ) 00964 CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT ) 00965 * 00966 * Multiply Q in U by real matrix RWORK(IRU), storing the 00967 * result in A, copying to U 00968 * (CWorkspace: need 0) 00969 * (Rworkspace: need N*N+2*M*N) 00970 * 00971 NRWORK = IRVT 00972 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 00973 $ RWORK( NRWORK ) ) 00974 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 00975 ELSE 00976 * 00977 * Copy A to VT, generate P**H 00978 * (Cworkspace: need 2*N, prefer N+N*NB) 00979 * (Rworkspace: 0) 00980 * 00981 CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT ) 00982 CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ), 00983 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00984 * 00985 * Copy A to U, generate Q 00986 * (Cworkspace: need 2*N, prefer N+N*NB) 00987 * (Rworkspace: 0) 00988 * 00989 CALL ZLACPY( 'L', M, N, A, LDA, U, LDU ) 00990 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 00991 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 00992 * 00993 * Perform bidiagonal SVD, computing left singular vectors 00994 * of bidiagonal matrix in RWORK(IRU) and computing right 00995 * singular vectors of bidiagonal matrix in RWORK(IRVT) 00996 * (CWorkspace: need 0) 00997 * (RWorkspace: need BDSPAC) 00998 * 00999 IRU = NRWORK 01000 IRVT = IRU + N*N 01001 NRWORK = IRVT + N*N 01002 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01003 $ N, RWORK( IRVT ), N, DUM, IDUM, 01004 $ RWORK( NRWORK ), IWORK, INFO ) 01005 * 01006 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01007 * storing the result in A, copying to VT 01008 * (Cworkspace: need 0) 01009 * (Rworkspace: need 3*N*N) 01010 * 01011 CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA, 01012 $ RWORK( NRWORK ) ) 01013 CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT ) 01014 * 01015 * Multiply Q in U by real matrix RWORK(IRU), storing the 01016 * result in A, copying to U 01017 * (CWorkspace: 0) 01018 * (Rworkspace: need 3*N*N) 01019 * 01020 NRWORK = IRVT 01021 CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA, 01022 $ RWORK( NRWORK ) ) 01023 CALL ZLACPY( 'F', M, N, A, LDA, U, LDU ) 01024 END IF 01025 * 01026 ELSE 01027 * 01028 * M .LT. MNTHR2 01029 * 01030 * Path 6 (M at least N, but not much larger) 01031 * Reduce to bidiagonal form without QR decomposition 01032 * Use ZUNMBR to compute singular vectors 01033 * 01034 IE = 1 01035 NRWORK = IE + N 01036 ITAUQ = 1 01037 ITAUP = ITAUQ + N 01038 NWORK = ITAUP + N 01039 * 01040 * Bidiagonalize A 01041 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB) 01042 * (RWorkspace: need N) 01043 * 01044 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01045 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01046 $ IERR ) 01047 IF( WNTQN ) THEN 01048 * 01049 * Compute singular values only 01050 * (Cworkspace: 0) 01051 * (Rworkspace: need BDSPAN) 01052 * 01053 CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1, 01054 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01055 ELSE IF( WNTQO ) THEN 01056 IU = NWORK 01057 IRU = NRWORK 01058 IRVT = IRU + N*N 01059 NRWORK = IRVT + N*N 01060 IF( LWORK.GE.M*N+3*N ) THEN 01061 * 01062 * WORK( IU ) is M by N 01063 * 01064 LDWRKU = M 01065 ELSE 01066 * 01067 * WORK( IU ) is LDWRKU by N 01068 * 01069 LDWRKU = ( LWORK-3*N ) / N 01070 END IF 01071 NWORK = IU + LDWRKU*N 01072 * 01073 * Perform bidiagonal SVD, computing left singular vectors 01074 * of bidiagonal matrix in RWORK(IRU) and computing right 01075 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01076 * (CWorkspace: need 0) 01077 * (RWorkspace: need BDSPAC) 01078 * 01079 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01080 $ N, RWORK( IRVT ), N, DUM, IDUM, 01081 $ RWORK( NRWORK ), IWORK, INFO ) 01082 * 01083 * Copy real matrix RWORK(IRVT) to complex matrix VT 01084 * Overwrite VT by right singular vectors of A 01085 * (Cworkspace: need 2*N, prefer N+N*NB) 01086 * (Rworkspace: need 0) 01087 * 01088 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01089 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01090 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01091 $ LWORK-NWORK+1, IERR ) 01092 * 01093 IF( LWORK.GE.M*N+3*N ) THEN 01094 * 01095 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 01096 * Overwrite WORK(IU) by left singular vectors of A, copying 01097 * to A 01098 * (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB) 01099 * (Rworkspace: need 0) 01100 * 01101 CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ), 01102 $ LDWRKU ) 01103 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ), 01104 $ LDWRKU ) 01105 CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 01106 $ WORK( ITAUQ ), WORK( IU ), LDWRKU, 01107 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01108 CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA ) 01109 ELSE 01110 * 01111 * Generate Q in A 01112 * (Cworkspace: need 2*N, prefer N+N*NB) 01113 * (Rworkspace: need 0) 01114 * 01115 CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ), 01116 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01117 * 01118 * Multiply Q in A by real matrix RWORK(IRU), storing the 01119 * result in WORK(IU), copying to A 01120 * (CWorkspace: need N*N, prefer M*N) 01121 * (Rworkspace: need 3*N*N, prefer N*N+2*M*N) 01122 * 01123 NRWORK = IRVT 01124 DO 30 I = 1, M, LDWRKU 01125 CHUNK = MIN( M-I+1, LDWRKU ) 01126 CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, 01127 $ RWORK( IRU ), N, WORK( IU ), LDWRKU, 01128 $ RWORK( NRWORK ) ) 01129 CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU, 01130 $ A( I, 1 ), LDA ) 01131 30 CONTINUE 01132 END IF 01133 * 01134 ELSE IF( WNTQS ) THEN 01135 * 01136 * Perform bidiagonal SVD, computing left singular vectors 01137 * of bidiagonal matrix in RWORK(IRU) and computing right 01138 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01139 * (CWorkspace: need 0) 01140 * (RWorkspace: need BDSPAC) 01141 * 01142 IRU = NRWORK 01143 IRVT = IRU + N*N 01144 NRWORK = IRVT + N*N 01145 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01146 $ N, RWORK( IRVT ), N, DUM, IDUM, 01147 $ RWORK( NRWORK ), IWORK, INFO ) 01148 * 01149 * Copy real matrix RWORK(IRU) to complex matrix U 01150 * Overwrite U by left singular vectors of A 01151 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01152 * (RWorkspace: 0) 01153 * 01154 CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU ) 01155 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 01156 CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA, 01157 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01158 $ LWORK-NWORK+1, IERR ) 01159 * 01160 * Copy real matrix RWORK(IRVT) to complex matrix VT 01161 * Overwrite VT by right singular vectors of A 01162 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01163 * (RWorkspace: 0) 01164 * 01165 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01166 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01167 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01168 $ LWORK-NWORK+1, IERR ) 01169 ELSE 01170 * 01171 * Perform bidiagonal SVD, computing left singular vectors 01172 * of bidiagonal matrix in RWORK(IRU) and computing right 01173 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01174 * (CWorkspace: need 0) 01175 * (RWorkspace: need BDSPAC) 01176 * 01177 IRU = NRWORK 01178 IRVT = IRU + N*N 01179 NRWORK = IRVT + N*N 01180 CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ), 01181 $ N, RWORK( IRVT ), N, DUM, IDUM, 01182 $ RWORK( NRWORK ), IWORK, INFO ) 01183 * 01184 * Set the right corner of U to identity matrix 01185 * 01186 CALL ZLASET( 'F', M, M, CZERO, CZERO, U, LDU ) 01187 IF( M.GT.N ) THEN 01188 CALL ZLASET( 'F', M-N, M-N, CZERO, CONE, 01189 $ U( N+1, N+1 ), LDU ) 01190 END IF 01191 * 01192 * Copy real matrix RWORK(IRU) to complex matrix U 01193 * Overwrite U by left singular vectors of A 01194 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB) 01195 * (RWorkspace: 0) 01196 * 01197 CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU ) 01198 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01199 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01200 $ LWORK-NWORK+1, IERR ) 01201 * 01202 * Copy real matrix RWORK(IRVT) to complex matrix VT 01203 * Overwrite VT by right singular vectors of A 01204 * (CWorkspace: need 3*N, prefer 2*N+N*NB) 01205 * (RWorkspace: 0) 01206 * 01207 CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT ) 01208 CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA, 01209 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01210 $ LWORK-NWORK+1, IERR ) 01211 END IF 01212 * 01213 END IF 01214 * 01215 ELSE 01216 * 01217 * A has more columns than rows. If A has sufficiently more 01218 * columns than rows, first reduce using the LQ decomposition (if 01219 * sufficient workspace available) 01220 * 01221 IF( N.GE.MNTHR1 ) THEN 01222 * 01223 IF( WNTQN ) THEN 01224 * 01225 * Path 1t (N much larger than M, JOBZ='N') 01226 * No singular vectors to be computed 01227 * 01228 ITAU = 1 01229 NWORK = ITAU + M 01230 * 01231 * Compute A=L*Q 01232 * (CWorkspace: need 2*M, prefer M+M*NB) 01233 * (RWorkspace: 0) 01234 * 01235 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01236 $ LWORK-NWORK+1, IERR ) 01237 * 01238 * Zero out above L 01239 * 01240 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 01241 $ LDA ) 01242 IE = 1 01243 ITAUQ = 1 01244 ITAUP = ITAUQ + M 01245 NWORK = ITAUP + M 01246 * 01247 * Bidiagonalize L in A 01248 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB) 01249 * (RWorkspace: need M) 01250 * 01251 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01252 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01253 $ IERR ) 01254 NRWORK = IE + M 01255 * 01256 * Perform bidiagonal SVD, compute singular values only 01257 * (CWorkspace: 0) 01258 * (RWorkspace: need BDSPAN) 01259 * 01260 CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01261 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01262 * 01263 ELSE IF( WNTQO ) THEN 01264 * 01265 * Path 2t (N much larger than M, JOBZ='O') 01266 * M right singular vectors to be overwritten on A and 01267 * M left singular vectors to be computed in U 01268 * 01269 IVT = 1 01270 LDWKVT = M 01271 * 01272 * WORK(IVT) is M by M 01273 * 01274 IL = IVT + LDWKVT*M 01275 IF( LWORK.GE.M*N+M*M+3*M ) THEN 01276 * 01277 * WORK(IL) M by N 01278 * 01279 LDWRKL = M 01280 CHUNK = N 01281 ELSE 01282 * 01283 * WORK(IL) is M by CHUNK 01284 * 01285 LDWRKL = M 01286 CHUNK = ( LWORK-M*M-3*M ) / M 01287 END IF 01288 ITAU = IL + LDWRKL*CHUNK 01289 NWORK = ITAU + M 01290 * 01291 * Compute A=L*Q 01292 * (CWorkspace: need 2*M, prefer M+M*NB) 01293 * (RWorkspace: 0) 01294 * 01295 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01296 $ LWORK-NWORK+1, IERR ) 01297 * 01298 * Copy L to WORK(IL), zeroing about above it 01299 * 01300 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01301 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 01302 $ WORK( IL+LDWRKL ), LDWRKL ) 01303 * 01304 * Generate Q in A 01305 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 01306 * (RWorkspace: 0) 01307 * 01308 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 01309 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01310 IE = 1 01311 ITAUQ = ITAU 01312 ITAUP = ITAUQ + M 01313 NWORK = ITAUP + M 01314 * 01315 * Bidiagonalize L in WORK(IL) 01316 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01317 * (RWorkspace: need M) 01318 * 01319 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 01320 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01321 $ LWORK-NWORK+1, IERR ) 01322 * 01323 * Perform bidiagonal SVD, computing left singular vectors 01324 * of bidiagonal matrix in RWORK(IRU) and computing right 01325 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01326 * (CWorkspace: need 0) 01327 * (RWorkspace: need BDSPAC) 01328 * 01329 IRU = IE + M 01330 IRVT = IRU + M*M 01331 NRWORK = IRVT + M*M 01332 CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01333 $ M, RWORK( IRVT ), M, DUM, IDUM, 01334 $ RWORK( NRWORK ), IWORK, INFO ) 01335 * 01336 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU) 01337 * Overwrite WORK(IU) by the left singular vectors of L 01338 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 01339 * (RWorkspace: 0) 01340 * 01341 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01342 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01343 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01344 $ LWORK-NWORK+1, IERR ) 01345 * 01346 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01347 * Overwrite WORK(IVT) by the right singular vectors of L 01348 * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB) 01349 * (RWorkspace: 0) 01350 * 01351 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01352 $ LDWKVT ) 01353 CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 01354 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01355 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01356 * 01357 * Multiply right singular vectors of L in WORK(IL) by Q 01358 * in A, storing result in WORK(IL) and copying to A 01359 * (CWorkspace: need 2*M*M, prefer M*M+M*N)) 01360 * (RWorkspace: 0) 01361 * 01362 DO 40 I = 1, N, CHUNK 01363 BLK = MIN( N-I+1, CHUNK ) 01364 CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M, 01365 $ A( 1, I ), LDA, CZERO, WORK( IL ), 01366 $ LDWRKL ) 01367 CALL ZLACPY( 'F', M, BLK, WORK( IL ), LDWRKL, 01368 $ A( 1, I ), LDA ) 01369 40 CONTINUE 01370 * 01371 ELSE IF( WNTQS ) THEN 01372 * 01373 * Path 3t (N much larger than M, JOBZ='S') 01374 * M right singular vectors to be computed in VT and 01375 * M left singular vectors to be computed in U 01376 * 01377 IL = 1 01378 * 01379 * WORK(IL) is M by M 01380 * 01381 LDWRKL = M 01382 ITAU = IL + LDWRKL*M 01383 NWORK = ITAU + M 01384 * 01385 * Compute A=L*Q 01386 * (CWorkspace: need 2*M, prefer M+M*NB) 01387 * (RWorkspace: 0) 01388 * 01389 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01390 $ LWORK-NWORK+1, IERR ) 01391 * 01392 * Copy L to WORK(IL), zeroing out above it 01393 * 01394 CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL ) 01395 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, 01396 $ WORK( IL+LDWRKL ), LDWRKL ) 01397 * 01398 * Generate Q in A 01399 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB) 01400 * (RWorkspace: 0) 01401 * 01402 CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ), 01403 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01404 IE = 1 01405 ITAUQ = ITAU 01406 ITAUP = ITAUQ + M 01407 NWORK = ITAUP + M 01408 * 01409 * Bidiagonalize L in WORK(IL) 01410 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01411 * (RWorkspace: need M) 01412 * 01413 CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ), 01414 $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ), 01415 $ LWORK-NWORK+1, IERR ) 01416 * 01417 * Perform bidiagonal SVD, computing left singular vectors 01418 * of bidiagonal matrix in RWORK(IRU) and computing right 01419 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01420 * (CWorkspace: need 0) 01421 * (RWorkspace: need BDSPAC) 01422 * 01423 IRU = IE + M 01424 IRVT = IRU + M*M 01425 NRWORK = IRVT + M*M 01426 CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01427 $ M, RWORK( IRVT ), M, DUM, IDUM, 01428 $ RWORK( NRWORK ), IWORK, INFO ) 01429 * 01430 * Copy real matrix RWORK(IRU) to complex matrix U 01431 * Overwrite U by left singular vectors of L 01432 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01433 * (RWorkspace: 0) 01434 * 01435 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01436 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL, 01437 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01438 $ LWORK-NWORK+1, IERR ) 01439 * 01440 * Copy real matrix RWORK(IRVT) to complex matrix VT 01441 * Overwrite VT by left singular vectors of L 01442 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01443 * (RWorkspace: 0) 01444 * 01445 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 01446 CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL, 01447 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01448 $ LWORK-NWORK+1, IERR ) 01449 * 01450 * Copy VT to WORK(IL), multiply right singular vectors of L 01451 * in WORK(IL) by Q in A, storing result in VT 01452 * (CWorkspace: need M*M) 01453 * (RWorkspace: 0) 01454 * 01455 CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL ) 01456 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL, 01457 $ A, LDA, CZERO, VT, LDVT ) 01458 * 01459 ELSE IF( WNTQA ) THEN 01460 * 01461 * Path 9t (N much larger than M, JOBZ='A') 01462 * N right singular vectors to be computed in VT and 01463 * M left singular vectors to be computed in U 01464 * 01465 IVT = 1 01466 * 01467 * WORK(IVT) is M by M 01468 * 01469 LDWKVT = M 01470 ITAU = IVT + LDWKVT*M 01471 NWORK = ITAU + M 01472 * 01473 * Compute A=L*Q, copying result to VT 01474 * (CWorkspace: need 2*M, prefer M+M*NB) 01475 * (RWorkspace: 0) 01476 * 01477 CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ), 01478 $ LWORK-NWORK+1, IERR ) 01479 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01480 * 01481 * Generate Q in VT 01482 * (CWorkspace: need M+N, prefer M+N*NB) 01483 * (RWorkspace: 0) 01484 * 01485 CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ), 01486 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01487 * 01488 * Produce L in A, zeroing out above it 01489 * 01490 CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ), 01491 $ LDA ) 01492 IE = 1 01493 ITAUQ = ITAU 01494 ITAUP = ITAUQ + M 01495 NWORK = ITAUP + M 01496 * 01497 * Bidiagonalize L in A 01498 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB) 01499 * (RWorkspace: need M) 01500 * 01501 CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01502 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01503 $ IERR ) 01504 * 01505 * Perform bidiagonal SVD, computing left singular vectors 01506 * of bidiagonal matrix in RWORK(IRU) and computing right 01507 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01508 * (CWorkspace: need 0) 01509 * (RWorkspace: need BDSPAC) 01510 * 01511 IRU = IE + M 01512 IRVT = IRU + M*M 01513 NRWORK = IRVT + M*M 01514 CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01515 $ M, RWORK( IRVT ), M, DUM, IDUM, 01516 $ RWORK( NRWORK ), IWORK, INFO ) 01517 * 01518 * Copy real matrix RWORK(IRU) to complex matrix U 01519 * Overwrite U by left singular vectors of L 01520 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01521 * (RWorkspace: 0) 01522 * 01523 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01524 CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA, 01525 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01526 $ LWORK-NWORK+1, IERR ) 01527 * 01528 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01529 * Overwrite WORK(IVT) by right singular vectors of L 01530 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB) 01531 * (RWorkspace: 0) 01532 * 01533 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01534 $ LDWKVT ) 01535 CALL ZUNMBR( 'P', 'R', 'C', M, M, M, A, LDA, 01536 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01537 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01538 * 01539 * Multiply right singular vectors of L in WORK(IVT) by 01540 * Q in VT, storing result in A 01541 * (CWorkspace: need M*M) 01542 * (RWorkspace: 0) 01543 * 01544 CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT, 01545 $ VT, LDVT, CZERO, A, LDA ) 01546 * 01547 * Copy right singular vectors of A from A to VT 01548 * 01549 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01550 * 01551 END IF 01552 * 01553 ELSE IF( N.GE.MNTHR2 ) THEN 01554 * 01555 * MNTHR2 <= N < MNTHR1 01556 * 01557 * Path 5t (N much larger than M, but not as much as MNTHR1) 01558 * Reduce to bidiagonal form without QR decomposition, use 01559 * ZUNGBR and matrix multiplication to compute singular vectors 01560 * 01561 * 01562 IE = 1 01563 NRWORK = IE + M 01564 ITAUQ = 1 01565 ITAUP = ITAUQ + M 01566 NWORK = ITAUP + M 01567 * 01568 * Bidiagonalize A 01569 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 01570 * (RWorkspace: M) 01571 * 01572 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01573 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01574 $ IERR ) 01575 * 01576 IF( WNTQN ) THEN 01577 * 01578 * Compute singular values only 01579 * (Cworkspace: 0) 01580 * (Rworkspace: need BDSPAN) 01581 * 01582 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01583 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01584 ELSE IF( WNTQO ) THEN 01585 IRVT = NRWORK 01586 IRU = IRVT + M*M 01587 NRWORK = IRU + M*M 01588 IVT = NWORK 01589 * 01590 * Copy A to U, generate Q 01591 * (Cworkspace: need 2*M, prefer M+M*NB) 01592 * (Rworkspace: 0) 01593 * 01594 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 01595 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01596 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01597 * 01598 * Generate P**H in A 01599 * (Cworkspace: need 2*M, prefer M+M*NB) 01600 * (Rworkspace: 0) 01601 * 01602 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 01603 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01604 * 01605 LDWKVT = M 01606 IF( LWORK.GE.M*N+3*M ) THEN 01607 * 01608 * WORK( IVT ) is M by N 01609 * 01610 NWORK = IVT + LDWKVT*N 01611 CHUNK = N 01612 ELSE 01613 * 01614 * WORK( IVT ) is M by CHUNK 01615 * 01616 CHUNK = ( LWORK-3*M ) / M 01617 NWORK = IVT + LDWKVT*CHUNK 01618 END IF 01619 * 01620 * Perform bidiagonal SVD, computing left singular vectors 01621 * of bidiagonal matrix in RWORK(IRU) and computing right 01622 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01623 * (CWorkspace: need 0) 01624 * (RWorkspace: need BDSPAC) 01625 * 01626 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01627 $ M, RWORK( IRVT ), M, DUM, IDUM, 01628 $ RWORK( NRWORK ), IWORK, INFO ) 01629 * 01630 * Multiply Q in U by real matrix RWORK(IRVT) 01631 * storing the result in WORK(IVT), copying to U 01632 * (Cworkspace: need 0) 01633 * (Rworkspace: need 2*M*M) 01634 * 01635 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ), 01636 $ LDWKVT, RWORK( NRWORK ) ) 01637 CALL ZLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU ) 01638 * 01639 * Multiply RWORK(IRVT) by P**H in A, storing the 01640 * result in WORK(IVT), copying to A 01641 * (CWorkspace: need M*M, prefer M*N) 01642 * (Rworkspace: need 2*M*M, prefer 2*M*N) 01643 * 01644 NRWORK = IRU 01645 DO 50 I = 1, N, CHUNK 01646 BLK = MIN( N-I+1, CHUNK ) 01647 CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA, 01648 $ WORK( IVT ), LDWKVT, RWORK( NRWORK ) ) 01649 CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 01650 $ A( 1, I ), LDA ) 01651 50 CONTINUE 01652 ELSE IF( WNTQS ) THEN 01653 * 01654 * Copy A to U, generate Q 01655 * (Cworkspace: need 2*M, prefer M+M*NB) 01656 * (Rworkspace: 0) 01657 * 01658 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 01659 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01660 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01661 * 01662 * Copy A to VT, generate P**H 01663 * (Cworkspace: need 2*M, prefer M+M*NB) 01664 * (Rworkspace: 0) 01665 * 01666 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01667 CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ), 01668 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01669 * 01670 * Perform bidiagonal SVD, computing left singular vectors 01671 * of bidiagonal matrix in RWORK(IRU) and computing right 01672 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01673 * (CWorkspace: need 0) 01674 * (RWorkspace: need BDSPAC) 01675 * 01676 IRVT = NRWORK 01677 IRU = IRVT + M*M 01678 NRWORK = IRU + M*M 01679 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01680 $ M, RWORK( IRVT ), M, DUM, IDUM, 01681 $ RWORK( NRWORK ), IWORK, INFO ) 01682 * 01683 * Multiply Q in U by real matrix RWORK(IRU), storing the 01684 * result in A, copying to U 01685 * (CWorkspace: need 0) 01686 * (Rworkspace: need 3*M*M) 01687 * 01688 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 01689 $ RWORK( NRWORK ) ) 01690 CALL ZLACPY( 'F', M, M, A, LDA, U, LDU ) 01691 * 01692 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01693 * storing the result in A, copying to VT 01694 * (Cworkspace: need 0) 01695 * (Rworkspace: need M*M+2*M*N) 01696 * 01697 NRWORK = IRU 01698 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 01699 $ RWORK( NRWORK ) ) 01700 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01701 ELSE 01702 * 01703 * Copy A to U, generate Q 01704 * (Cworkspace: need 2*M, prefer M+M*NB) 01705 * (Rworkspace: 0) 01706 * 01707 CALL ZLACPY( 'L', M, M, A, LDA, U, LDU ) 01708 CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ), 01709 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01710 * 01711 * Copy A to VT, generate P**H 01712 * (Cworkspace: need 2*M, prefer M+M*NB) 01713 * (Rworkspace: 0) 01714 * 01715 CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT ) 01716 CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ), 01717 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01718 * 01719 * Perform bidiagonal SVD, computing left singular vectors 01720 * of bidiagonal matrix in RWORK(IRU) and computing right 01721 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01722 * (CWorkspace: need 0) 01723 * (RWorkspace: need BDSPAC) 01724 * 01725 IRVT = NRWORK 01726 IRU = IRVT + M*M 01727 NRWORK = IRU + M*M 01728 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01729 $ M, RWORK( IRVT ), M, DUM, IDUM, 01730 $ RWORK( NRWORK ), IWORK, INFO ) 01731 * 01732 * Multiply Q in U by real matrix RWORK(IRU), storing the 01733 * result in A, copying to U 01734 * (CWorkspace: need 0) 01735 * (Rworkspace: need 3*M*M) 01736 * 01737 CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA, 01738 $ RWORK( NRWORK ) ) 01739 CALL ZLACPY( 'F', M, M, A, LDA, U, LDU ) 01740 * 01741 * Multiply real matrix RWORK(IRVT) by P**H in VT, 01742 * storing the result in A, copying to VT 01743 * (Cworkspace: need 0) 01744 * (Rworkspace: need M*M+2*M*N) 01745 * 01746 CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA, 01747 $ RWORK( NRWORK ) ) 01748 CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT ) 01749 END IF 01750 * 01751 ELSE 01752 * 01753 * N .LT. MNTHR2 01754 * 01755 * Path 6t (N greater than M, but not much larger) 01756 * Reduce to bidiagonal form without LQ decomposition 01757 * Use ZUNMBR to compute singular vectors 01758 * 01759 IE = 1 01760 NRWORK = IE + M 01761 ITAUQ = 1 01762 ITAUP = ITAUQ + M 01763 NWORK = ITAUP + M 01764 * 01765 * Bidiagonalize A 01766 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB) 01767 * (RWorkspace: M) 01768 * 01769 CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ), 01770 $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1, 01771 $ IERR ) 01772 IF( WNTQN ) THEN 01773 * 01774 * Compute singular values only 01775 * (Cworkspace: 0) 01776 * (Rworkspace: need BDSPAN) 01777 * 01778 CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1, 01779 $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO ) 01780 ELSE IF( WNTQO ) THEN 01781 LDWKVT = M 01782 IVT = NWORK 01783 IF( LWORK.GE.M*N+3*M ) THEN 01784 * 01785 * WORK( IVT ) is M by N 01786 * 01787 CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ), 01788 $ LDWKVT ) 01789 NWORK = IVT + LDWKVT*N 01790 ELSE 01791 * 01792 * WORK( IVT ) is M by CHUNK 01793 * 01794 CHUNK = ( LWORK-3*M ) / M 01795 NWORK = IVT + LDWKVT*CHUNK 01796 END IF 01797 * 01798 * Perform bidiagonal SVD, computing left singular vectors 01799 * of bidiagonal matrix in RWORK(IRU) and computing right 01800 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01801 * (CWorkspace: need 0) 01802 * (RWorkspace: need BDSPAC) 01803 * 01804 IRVT = NRWORK 01805 IRU = IRVT + M*M 01806 NRWORK = IRU + M*M 01807 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01808 $ M, RWORK( IRVT ), M, DUM, IDUM, 01809 $ RWORK( NRWORK ), IWORK, INFO ) 01810 * 01811 * Copy real matrix RWORK(IRU) to complex matrix U 01812 * Overwrite U by left singular vectors of A 01813 * (Cworkspace: need 2*M, prefer M+M*NB) 01814 * (Rworkspace: need 0) 01815 * 01816 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01817 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01818 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01819 $ LWORK-NWORK+1, IERR ) 01820 * 01821 IF( LWORK.GE.M*N+3*M ) THEN 01822 * 01823 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT) 01824 * Overwrite WORK(IVT) by right singular vectors of A, 01825 * copying to A 01826 * (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB) 01827 * (Rworkspace: need 0) 01828 * 01829 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ), 01830 $ LDWKVT ) 01831 CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 01832 $ WORK( ITAUP ), WORK( IVT ), LDWKVT, 01833 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01834 CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA ) 01835 ELSE 01836 * 01837 * Generate P**H in A 01838 * (Cworkspace: need 2*M, prefer M+M*NB) 01839 * (Rworkspace: need 0) 01840 * 01841 CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ), 01842 $ WORK( NWORK ), LWORK-NWORK+1, IERR ) 01843 * 01844 * Multiply Q in A by real matrix RWORK(IRU), storing the 01845 * result in WORK(IU), copying to A 01846 * (CWorkspace: need M*M, prefer M*N) 01847 * (Rworkspace: need 3*M*M, prefer M*M+2*M*N) 01848 * 01849 NRWORK = IRU 01850 DO 60 I = 1, N, CHUNK 01851 BLK = MIN( N-I+1, CHUNK ) 01852 CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), 01853 $ LDA, WORK( IVT ), LDWKVT, 01854 $ RWORK( NRWORK ) ) 01855 CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT, 01856 $ A( 1, I ), LDA ) 01857 60 CONTINUE 01858 END IF 01859 ELSE IF( WNTQS ) THEN 01860 * 01861 * Perform bidiagonal SVD, computing left singular vectors 01862 * of bidiagonal matrix in RWORK(IRU) and computing right 01863 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01864 * (CWorkspace: need 0) 01865 * (RWorkspace: need BDSPAC) 01866 * 01867 IRVT = NRWORK 01868 IRU = IRVT + M*M 01869 NRWORK = IRU + M*M 01870 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01871 $ M, RWORK( IRVT ), M, DUM, IDUM, 01872 $ RWORK( NRWORK ), IWORK, INFO ) 01873 * 01874 * Copy real matrix RWORK(IRU) to complex matrix U 01875 * Overwrite U by left singular vectors of A 01876 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01877 * (RWorkspace: M*M) 01878 * 01879 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01880 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01881 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01882 $ LWORK-NWORK+1, IERR ) 01883 * 01884 * Copy real matrix RWORK(IRVT) to complex matrix VT 01885 * Overwrite VT by right singular vectors of A 01886 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01887 * (RWorkspace: M*M) 01888 * 01889 CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT ) 01890 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 01891 CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA, 01892 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01893 $ LWORK-NWORK+1, IERR ) 01894 ELSE 01895 * 01896 * Perform bidiagonal SVD, computing left singular vectors 01897 * of bidiagonal matrix in RWORK(IRU) and computing right 01898 * singular vectors of bidiagonal matrix in RWORK(IRVT) 01899 * (CWorkspace: need 0) 01900 * (RWorkspace: need BDSPAC) 01901 * 01902 IRVT = NRWORK 01903 IRU = IRVT + M*M 01904 NRWORK = IRU + M*M 01905 * 01906 CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ), 01907 $ M, RWORK( IRVT ), M, DUM, IDUM, 01908 $ RWORK( NRWORK ), IWORK, INFO ) 01909 * 01910 * Copy real matrix RWORK(IRU) to complex matrix U 01911 * Overwrite U by left singular vectors of A 01912 * (CWorkspace: need 3*M, prefer 2*M+M*NB) 01913 * (RWorkspace: M*M) 01914 * 01915 CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU ) 01916 CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA, 01917 $ WORK( ITAUQ ), U, LDU, WORK( NWORK ), 01918 $ LWORK-NWORK+1, IERR ) 01919 * 01920 * Set all of VT to identity matrix 01921 * 01922 CALL ZLASET( 'F', N, N, CZERO, CONE, VT, LDVT ) 01923 * 01924 * Copy real matrix RWORK(IRVT) to complex matrix VT 01925 * Overwrite VT by right singular vectors of A 01926 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB) 01927 * (RWorkspace: M*M) 01928 * 01929 CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT ) 01930 CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA, 01931 $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ), 01932 $ LWORK-NWORK+1, IERR ) 01933 END IF 01934 * 01935 END IF 01936 * 01937 END IF 01938 * 01939 * Undo scaling if necessary 01940 * 01941 IF( ISCL.EQ.1 ) THEN 01942 IF( ANRM.GT.BIGNUM ) 01943 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN, 01944 $ IERR ) 01945 IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM ) 01946 $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1, 01947 $ RWORK( IE ), MINMN, IERR ) 01948 IF( ANRM.LT.SMLNUM ) 01949 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN, 01950 $ IERR ) 01951 IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM ) 01952 $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1, 01953 $ RWORK( IE ), MINMN, IERR ) 01954 END IF 01955 * 01956 * Return optimal workspace in WORK(1) 01957 * 01958 WORK( 1 ) = MAXWRK 01959 * 01960 RETURN 01961 * 01962 * End of ZGESDD 01963 * 01964 END
1.7.3 
