![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTGSJA 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTGSJA + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsja.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsja.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsja.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 00022 * LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 00023 * Q, LDQ, WORK, NCYCLE, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBQ, JOBU, JOBV 00027 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 00028 * $ NCYCLE, P 00029 * DOUBLE PRECISION TOLA, TOLB 00030 * .. 00031 * .. Array Arguments .. 00032 * DOUBLE PRECISION ALPHA( * ), BETA( * ) 00033 * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00034 * $ U( LDU, * ), V( LDV, * ), WORK( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZTGSJA computes the generalized singular value decomposition (GSVD) 00044 *> of two complex upper triangular (or trapezoidal) matrices A and B. 00045 *> 00046 *> On entry, it is assumed that matrices A and B have the following 00047 *> forms, which may be obtained by the preprocessing subroutine ZGGSVP 00048 *> from a general M-by-N matrix A and P-by-N matrix B: 00049 *> 00050 *> N-K-L K L 00051 *> A = K ( 0 A12 A13 ) if M-K-L >= 0; 00052 *> L ( 0 0 A23 ) 00053 *> M-K-L ( 0 0 0 ) 00054 *> 00055 *> N-K-L K L 00056 *> A = K ( 0 A12 A13 ) if M-K-L < 0; 00057 *> M-K ( 0 0 A23 ) 00058 *> 00059 *> N-K-L K L 00060 *> B = L ( 0 0 B13 ) 00061 *> P-L ( 0 0 0 ) 00062 *> 00063 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular 00064 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0, 00065 *> otherwise A23 is (M-K)-by-L upper trapezoidal. 00066 *> 00067 *> On exit, 00068 *> 00069 *> U**H *A*Q = D1*( 0 R ), V**H *B*Q = D2*( 0 R ), 00070 *> 00071 *> where U, V and Q are unitary matrices. 00072 *> R is a nonsingular upper triangular matrix, and D1 00073 *> and D2 are ``diagonal'' matrices, which are of the following 00074 *> structures: 00075 *> 00076 *> If M-K-L >= 0, 00077 *> 00078 *> K L 00079 *> D1 = K ( I 0 ) 00080 *> L ( 0 C ) 00081 *> M-K-L ( 0 0 ) 00082 *> 00083 *> K L 00084 *> D2 = L ( 0 S ) 00085 *> P-L ( 0 0 ) 00086 *> 00087 *> N-K-L K L 00088 *> ( 0 R ) = K ( 0 R11 R12 ) K 00089 *> L ( 0 0 R22 ) L 00090 *> 00091 *> where 00092 *> 00093 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), 00094 *> S = diag( BETA(K+1), ... , BETA(K+L) ), 00095 *> C**2 + S**2 = I. 00096 *> 00097 *> R is stored in A(1:K+L,N-K-L+1:N) on exit. 00098 *> 00099 *> If M-K-L < 0, 00100 *> 00101 *> K M-K K+L-M 00102 *> D1 = K ( I 0 0 ) 00103 *> M-K ( 0 C 0 ) 00104 *> 00105 *> K M-K K+L-M 00106 *> D2 = M-K ( 0 S 0 ) 00107 *> K+L-M ( 0 0 I ) 00108 *> P-L ( 0 0 0 ) 00109 *> 00110 *> N-K-L K M-K K+L-M 00111 *> ( 0 R ) = K ( 0 R11 R12 R13 ) 00112 *> M-K ( 0 0 R22 R23 ) 00113 *> K+L-M ( 0 0 0 R33 ) 00114 *> 00115 *> where 00116 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ), 00117 *> S = diag( BETA(K+1), ... , BETA(M) ), 00118 *> C**2 + S**2 = I. 00119 *> 00120 *> R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored 00121 *> ( 0 R22 R23 ) 00122 *> in B(M-K+1:L,N+M-K-L+1:N) on exit. 00123 *> 00124 *> The computation of the unitary transformation matrices U, V or Q 00125 *> is optional. These matrices may either be formed explicitly, or they 00126 *> may be postmultiplied into input matrices U1, V1, or Q1. 00127 *> \endverbatim 00128 * 00129 * Arguments: 00130 * ========== 00131 * 00132 *> \param[in] JOBU 00133 *> \verbatim 00134 *> JOBU is CHARACTER*1 00135 *> = 'U': U must contain a unitary matrix U1 on entry, and 00136 *> the product U1*U is returned; 00137 *> = 'I': U is initialized to the unit matrix, and the 00138 *> unitary matrix U is returned; 00139 *> = 'N': U is not computed. 00140 *> \endverbatim 00141 *> 00142 *> \param[in] JOBV 00143 *> \verbatim 00144 *> JOBV is CHARACTER*1 00145 *> = 'V': V must contain a unitary matrix V1 on entry, and 00146 *> the product V1*V is returned; 00147 *> = 'I': V is initialized to the unit matrix, and the 00148 *> unitary matrix V is returned; 00149 *> = 'N': V is not computed. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] JOBQ 00153 *> \verbatim 00154 *> JOBQ is CHARACTER*1 00155 *> = 'Q': Q must contain a unitary matrix Q1 on entry, and 00156 *> the product Q1*Q is returned; 00157 *> = 'I': Q is initialized to the unit matrix, and the 00158 *> unitary matrix Q is returned; 00159 *> = 'N': Q is not computed. 00160 *> \endverbatim 00161 *> 00162 *> \param[in] M 00163 *> \verbatim 00164 *> M is INTEGER 00165 *> The number of rows of the matrix A. M >= 0. 00166 *> \endverbatim 00167 *> 00168 *> \param[in] P 00169 *> \verbatim 00170 *> P is INTEGER 00171 *> The number of rows of the matrix B. P >= 0. 00172 *> \endverbatim 00173 *> 00174 *> \param[in] N 00175 *> \verbatim 00176 *> N is INTEGER 00177 *> The number of columns of the matrices A and B. N >= 0. 00178 *> \endverbatim 00179 *> 00180 *> \param[in] K 00181 *> \verbatim 00182 *> K is INTEGER 00183 *> \endverbatim 00184 *> 00185 *> \param[in] L 00186 *> \verbatim 00187 *> L is INTEGER 00188 *> 00189 *> K and L specify the subblocks in the input matrices A and B: 00190 *> A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N-L+1:N) 00191 *> of A and B, whose GSVD is going to be computed by ZTGSJA. 00192 *> See Further Details. 00193 *> \endverbatim 00194 *> 00195 *> \param[in,out] A 00196 *> \verbatim 00197 *> A is COMPLEX*16 array, dimension (LDA,N) 00198 *> On entry, the M-by-N matrix A. 00199 *> On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular 00200 *> matrix R or part of R. See Purpose for details. 00201 *> \endverbatim 00202 *> 00203 *> \param[in] LDA 00204 *> \verbatim 00205 *> LDA is INTEGER 00206 *> The leading dimension of the array A. LDA >= max(1,M). 00207 *> \endverbatim 00208 *> 00209 *> \param[in,out] B 00210 *> \verbatim 00211 *> B is COMPLEX*16 array, dimension (LDB,N) 00212 *> On entry, the P-by-N matrix B. 00213 *> On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains 00214 *> a part of R. See Purpose for details. 00215 *> \endverbatim 00216 *> 00217 *> \param[in] LDB 00218 *> \verbatim 00219 *> LDB is INTEGER 00220 *> The leading dimension of the array B. LDB >= max(1,P). 00221 *> \endverbatim 00222 *> 00223 *> \param[in] TOLA 00224 *> \verbatim 00225 *> TOLA is DOUBLE PRECISION 00226 *> \endverbatim 00227 *> 00228 *> \param[in] TOLB 00229 *> \verbatim 00230 *> TOLB is DOUBLE PRECISION 00231 *> 00232 *> TOLA and TOLB are the convergence criteria for the Jacobi- 00233 *> Kogbetliantz iteration procedure. Generally, they are the 00234 *> same as used in the preprocessing step, say 00235 *> TOLA = MAX(M,N)*norm(A)*MAZHEPS, 00236 *> TOLB = MAX(P,N)*norm(B)*MAZHEPS. 00237 *> \endverbatim 00238 *> 00239 *> \param[out] ALPHA 00240 *> \verbatim 00241 *> ALPHA is DOUBLE PRECISION array, dimension (N) 00242 *> \endverbatim 00243 *> 00244 *> \param[out] BETA 00245 *> \verbatim 00246 *> BETA is DOUBLE PRECISION array, dimension (N) 00247 *> 00248 *> On exit, ALPHA and BETA contain the generalized singular 00249 *> value pairs of A and B; 00250 *> ALPHA(1:K) = 1, 00251 *> BETA(1:K) = 0, 00252 *> and if M-K-L >= 0, 00253 *> ALPHA(K+1:K+L) = diag(C), 00254 *> BETA(K+1:K+L) = diag(S), 00255 *> or if M-K-L < 0, 00256 *> ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0 00257 *> BETA(K+1:M) = S, BETA(M+1:K+L) = 1. 00258 *> Furthermore, if K+L < N, 00259 *> ALPHA(K+L+1:N) = 0 and 00260 *> BETA(K+L+1:N) = 0. 00261 *> \endverbatim 00262 *> 00263 *> \param[in,out] U 00264 *> \verbatim 00265 *> U is COMPLEX*16 array, dimension (LDU,M) 00266 *> On entry, if JOBU = 'U', U must contain a matrix U1 (usually 00267 *> the unitary matrix returned by ZGGSVP). 00268 *> On exit, 00269 *> if JOBU = 'I', U contains the unitary matrix U; 00270 *> if JOBU = 'U', U contains the product U1*U. 00271 *> If JOBU = 'N', U is not referenced. 00272 *> \endverbatim 00273 *> 00274 *> \param[in] LDU 00275 *> \verbatim 00276 *> LDU is INTEGER 00277 *> The leading dimension of the array U. LDU >= max(1,M) if 00278 *> JOBU = 'U'; LDU >= 1 otherwise. 00279 *> \endverbatim 00280 *> 00281 *> \param[in,out] V 00282 *> \verbatim 00283 *> V is COMPLEX*16 array, dimension (LDV,P) 00284 *> On entry, if JOBV = 'V', V must contain a matrix V1 (usually 00285 *> the unitary matrix returned by ZGGSVP). 00286 *> On exit, 00287 *> if JOBV = 'I', V contains the unitary matrix V; 00288 *> if JOBV = 'V', V contains the product V1*V. 00289 *> If JOBV = 'N', V is not referenced. 00290 *> \endverbatim 00291 *> 00292 *> \param[in] LDV 00293 *> \verbatim 00294 *> LDV is INTEGER 00295 *> The leading dimension of the array V. LDV >= max(1,P) if 00296 *> JOBV = 'V'; LDV >= 1 otherwise. 00297 *> \endverbatim 00298 *> 00299 *> \param[in,out] Q 00300 *> \verbatim 00301 *> Q is COMPLEX*16 array, dimension (LDQ,N) 00302 *> On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually 00303 *> the unitary matrix returned by ZGGSVP). 00304 *> On exit, 00305 *> if JOBQ = 'I', Q contains the unitary matrix Q; 00306 *> if JOBQ = 'Q', Q contains the product Q1*Q. 00307 *> If JOBQ = 'N', Q is not referenced. 00308 *> \endverbatim 00309 *> 00310 *> \param[in] LDQ 00311 *> \verbatim 00312 *> LDQ is INTEGER 00313 *> The leading dimension of the array Q. LDQ >= max(1,N) if 00314 *> JOBQ = 'Q'; LDQ >= 1 otherwise. 00315 *> \endverbatim 00316 *> 00317 *> \param[out] WORK 00318 *> \verbatim 00319 *> WORK is COMPLEX*16 array, dimension (2*N) 00320 *> \endverbatim 00321 *> 00322 *> \param[out] NCYCLE 00323 *> \verbatim 00324 *> NCYCLE is INTEGER 00325 *> The number of cycles required for convergence. 00326 *> \endverbatim 00327 *> 00328 *> \param[out] INFO 00329 *> \verbatim 00330 *> INFO is INTEGER 00331 *> = 0: successful exit 00332 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00333 *> = 1: the procedure does not converge after MAXIT cycles. 00334 *> \endverbatim 00335 * 00336 *> \par Internal Parameters: 00337 * ========================= 00338 *> 00339 *> \verbatim 00340 *> MAXIT INTEGER 00341 *> MAXIT specifies the total loops that the iterative procedure 00342 *> may take. If after MAXIT cycles, the routine fails to 00343 *> converge, we return INFO = 1. 00344 *> \endverbatim 00345 * 00346 * Authors: 00347 * ======== 00348 * 00349 *> \author Univ. of Tennessee 00350 *> \author Univ. of California Berkeley 00351 *> \author Univ. of Colorado Denver 00352 *> \author NAG Ltd. 00353 * 00354 *> \date November 2011 00355 * 00356 *> \ingroup complex16OTHERcomputational 00357 * 00358 *> \par Further Details: 00359 * ===================== 00360 *> 00361 *> \verbatim 00362 *> 00363 *> ZTGSJA essentially uses a variant of Kogbetliantz algorithm to reduce 00364 *> min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L 00365 *> matrix B13 to the form: 00366 *> 00367 *> U1**H *A13*Q1 = C1*R1; V1**H *B13*Q1 = S1*R1, 00368 *> 00369 *> where U1, V1 and Q1 are unitary matrix. 00370 *> C1 and S1 are diagonal matrices satisfying 00371 *> 00372 *> C1**2 + S1**2 = I, 00373 *> 00374 *> and R1 is an L-by-L nonsingular upper triangular matrix. 00375 *> \endverbatim 00376 *> 00377 * ===================================================================== 00378 SUBROUTINE ZTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, 00379 $ LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, 00380 $ Q, LDQ, WORK, NCYCLE, INFO ) 00381 * 00382 * -- LAPACK computational routine (version 3.4.0) -- 00383 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00384 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00385 * November 2011 00386 * 00387 * .. Scalar Arguments .. 00388 CHARACTER JOBQ, JOBU, JOBV 00389 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, 00390 $ NCYCLE, P 00391 DOUBLE PRECISION TOLA, TOLB 00392 * .. 00393 * .. Array Arguments .. 00394 DOUBLE PRECISION ALPHA( * ), BETA( * ) 00395 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00396 $ U( LDU, * ), V( LDV, * ), WORK( * ) 00397 * .. 00398 * 00399 * ===================================================================== 00400 * 00401 * .. Parameters .. 00402 INTEGER MAXIT 00403 PARAMETER ( MAXIT = 40 ) 00404 DOUBLE PRECISION ZERO, ONE 00405 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00406 COMPLEX*16 CZERO, CONE 00407 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00408 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00409 * .. 00410 * .. Local Scalars .. 00411 * 00412 LOGICAL INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV 00413 INTEGER I, J, KCYCLE 00414 DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV, ERROR, GAMMA, 00415 $ RWK, SSMIN 00416 COMPLEX*16 A2, B2, SNQ, SNU, SNV 00417 * .. 00418 * .. External Functions .. 00419 LOGICAL LSAME 00420 EXTERNAL LSAME 00421 * .. 00422 * .. External Subroutines .. 00423 EXTERNAL DLARTG, XERBLA, ZCOPY, ZDSCAL, ZLAGS2, ZLAPLL, 00424 $ ZLASET, ZROT 00425 * .. 00426 * .. Intrinsic Functions .. 00427 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN 00428 * .. 00429 * .. Executable Statements .. 00430 * 00431 * Decode and test the input parameters 00432 * 00433 INITU = LSAME( JOBU, 'I' ) 00434 WANTU = INITU .OR. LSAME( JOBU, 'U' ) 00435 * 00436 INITV = LSAME( JOBV, 'I' ) 00437 WANTV = INITV .OR. LSAME( JOBV, 'V' ) 00438 * 00439 INITQ = LSAME( JOBQ, 'I' ) 00440 WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' ) 00441 * 00442 INFO = 0 00443 IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN 00444 INFO = -1 00445 ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00446 INFO = -2 00447 ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN 00448 INFO = -3 00449 ELSE IF( M.LT.0 ) THEN 00450 INFO = -4 00451 ELSE IF( P.LT.0 ) THEN 00452 INFO = -5 00453 ELSE IF( N.LT.0 ) THEN 00454 INFO = -6 00455 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00456 INFO = -10 00457 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00458 INFO = -12 00459 ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN 00460 INFO = -18 00461 ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN 00462 INFO = -20 00463 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN 00464 INFO = -22 00465 END IF 00466 IF( INFO.NE.0 ) THEN 00467 CALL XERBLA( 'ZTGSJA', -INFO ) 00468 RETURN 00469 END IF 00470 * 00471 * Initialize U, V and Q, if necessary 00472 * 00473 IF( INITU ) 00474 $ CALL ZLASET( 'Full', M, M, CZERO, CONE, U, LDU ) 00475 IF( INITV ) 00476 $ CALL ZLASET( 'Full', P, P, CZERO, CONE, V, LDV ) 00477 IF( INITQ ) 00478 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ ) 00479 * 00480 * Loop until convergence 00481 * 00482 UPPER = .FALSE. 00483 DO 40 KCYCLE = 1, MAXIT 00484 * 00485 UPPER = .NOT.UPPER 00486 * 00487 DO 20 I = 1, L - 1 00488 DO 10 J = I + 1, L 00489 * 00490 A1 = ZERO 00491 A2 = CZERO 00492 A3 = ZERO 00493 IF( K+I.LE.M ) 00494 $ A1 = DBLE( A( K+I, N-L+I ) ) 00495 IF( K+J.LE.M ) 00496 $ A3 = DBLE( A( K+J, N-L+J ) ) 00497 * 00498 B1 = DBLE( B( I, N-L+I ) ) 00499 B3 = DBLE( B( J, N-L+J ) ) 00500 * 00501 IF( UPPER ) THEN 00502 IF( K+I.LE.M ) 00503 $ A2 = A( K+I, N-L+J ) 00504 B2 = B( I, N-L+J ) 00505 ELSE 00506 IF( K+J.LE.M ) 00507 $ A2 = A( K+J, N-L+I ) 00508 B2 = B( J, N-L+I ) 00509 END IF 00510 * 00511 CALL ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, 00512 $ CSV, SNV, CSQ, SNQ ) 00513 * 00514 * Update (K+I)-th and (K+J)-th rows of matrix A: U**H *A 00515 * 00516 IF( K+J.LE.M ) 00517 $ CALL ZROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ), 00518 $ LDA, CSU, DCONJG( SNU ) ) 00519 * 00520 * Update I-th and J-th rows of matrix B: V**H *B 00521 * 00522 CALL ZROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB, 00523 $ CSV, DCONJG( SNV ) ) 00524 * 00525 * Update (N-L+I)-th and (N-L+J)-th columns of matrices 00526 * A and B: A*Q and B*Q 00527 * 00528 CALL ZROT( MIN( K+L, M ), A( 1, N-L+J ), 1, 00529 $ A( 1, N-L+I ), 1, CSQ, SNQ ) 00530 * 00531 CALL ZROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ, 00532 $ SNQ ) 00533 * 00534 IF( UPPER ) THEN 00535 IF( K+I.LE.M ) 00536 $ A( K+I, N-L+J ) = CZERO 00537 B( I, N-L+J ) = CZERO 00538 ELSE 00539 IF( K+J.LE.M ) 00540 $ A( K+J, N-L+I ) = CZERO 00541 B( J, N-L+I ) = CZERO 00542 END IF 00543 * 00544 * Ensure that the diagonal elements of A and B are real. 00545 * 00546 IF( K+I.LE.M ) 00547 $ A( K+I, N-L+I ) = DBLE( A( K+I, N-L+I ) ) 00548 IF( K+J.LE.M ) 00549 $ A( K+J, N-L+J ) = DBLE( A( K+J, N-L+J ) ) 00550 B( I, N-L+I ) = DBLE( B( I, N-L+I ) ) 00551 B( J, N-L+J ) = DBLE( B( J, N-L+J ) ) 00552 * 00553 * Update unitary matrices U, V, Q, if desired. 00554 * 00555 IF( WANTU .AND. K+J.LE.M ) 00556 $ CALL ZROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU, 00557 $ SNU ) 00558 * 00559 IF( WANTV ) 00560 $ CALL ZROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV ) 00561 * 00562 IF( WANTQ ) 00563 $ CALL ZROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ, 00564 $ SNQ ) 00565 * 00566 10 CONTINUE 00567 20 CONTINUE 00568 * 00569 IF( .NOT.UPPER ) THEN 00570 * 00571 * The matrices A13 and B13 were lower triangular at the start 00572 * of the cycle, and are now upper triangular. 00573 * 00574 * Convergence test: test the parallelism of the corresponding 00575 * rows of A and B. 00576 * 00577 ERROR = ZERO 00578 DO 30 I = 1, MIN( L, M-K ) 00579 CALL ZCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 ) 00580 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 ) 00581 CALL ZLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN ) 00582 ERROR = MAX( ERROR, SSMIN ) 00583 30 CONTINUE 00584 * 00585 IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) ) 00586 $ GO TO 50 00587 END IF 00588 * 00589 * End of cycle loop 00590 * 00591 40 CONTINUE 00592 * 00593 * The algorithm has not converged after MAXIT cycles. 00594 * 00595 INFO = 1 00596 GO TO 100 00597 * 00598 50 CONTINUE 00599 * 00600 * If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged. 00601 * Compute the generalized singular value pairs (ALPHA, BETA), and 00602 * set the triangular matrix R to array A. 00603 * 00604 DO 60 I = 1, K 00605 ALPHA( I ) = ONE 00606 BETA( I ) = ZERO 00607 60 CONTINUE 00608 * 00609 DO 70 I = 1, MIN( L, M-K ) 00610 * 00611 A1 = DBLE( A( K+I, N-L+I ) ) 00612 B1 = DBLE( B( I, N-L+I ) ) 00613 * 00614 IF( A1.NE.ZERO ) THEN 00615 GAMMA = B1 / A1 00616 * 00617 IF( GAMMA.LT.ZERO ) THEN 00618 CALL ZDSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB ) 00619 IF( WANTV ) 00620 $ CALL ZDSCAL( P, -ONE, V( 1, I ), 1 ) 00621 END IF 00622 * 00623 CALL DLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ), 00624 $ RWK ) 00625 * 00626 IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN 00627 CALL ZDSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ), 00628 $ LDA ) 00629 ELSE 00630 CALL ZDSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ), 00631 $ LDB ) 00632 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00633 $ LDA ) 00634 END IF 00635 * 00636 ELSE 00637 * 00638 ALPHA( K+I ) = ZERO 00639 BETA( K+I ) = ONE 00640 CALL ZCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ), 00641 $ LDA ) 00642 END IF 00643 70 CONTINUE 00644 * 00645 * Post-assignment 00646 * 00647 DO 80 I = M + 1, K + L 00648 ALPHA( I ) = ZERO 00649 BETA( I ) = ONE 00650 80 CONTINUE 00651 * 00652 IF( K+L.LT.N ) THEN 00653 DO 90 I = K + L + 1, N 00654 ALPHA( I ) = ZERO 00655 BETA( I ) = ZERO 00656 90 CONTINUE 00657 END IF 00658 * 00659 100 CONTINUE 00660 NCYCLE = KCYCLE 00661 * 00662 RETURN 00663 * 00664 * End of ZTGSJA 00665 * 00666 END