|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 00002 $ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, 00003 $ IWORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * January 2007 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER TRANS 00012 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, 00013 $ LWORK, M, N 00014 DOUBLE PRECISION DIF, SCALE 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IWORK( * ) 00018 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), 00019 $ D( LDD, * ), E( LDE, * ), F( LDF, * ), 00020 $ WORK( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZTGSYL solves the generalized Sylvester equation: 00027 * 00028 * A * R - L * B = scale * C (1) 00029 * D * R - L * E = scale * F 00030 * 00031 * where R and L are unknown m-by-n matrices, (A, D), (B, E) and 00032 * (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n, 00033 * respectively, with complex entries. A, B, D and E are upper 00034 * triangular (i.e., (A,D) and (B,E) in generalized Schur form). 00035 * 00036 * The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 00037 * is an output scaling factor chosen to avoid overflow. 00038 * 00039 * In matrix notation (1) is equivalent to solve Zx = scale*b, where Z 00040 * is defined as 00041 * 00042 * Z = [ kron(In, A) -kron(B', Im) ] (2) 00043 * [ kron(In, D) -kron(E', Im) ], 00044 * 00045 * Here Ix is the identity matrix of size x and X' is the conjugate 00046 * transpose of X. Kron(X, Y) is the Kronecker product between the 00047 * matrices X and Y. 00048 * 00049 * If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b 00050 * is solved for, which is equivalent to solve for R and L in 00051 * 00052 * A' * R + D' * L = scale * C (3) 00053 * R * B' + L * E' = scale * -F 00054 * 00055 * This case (TRANS = 'C') is used to compute an one-norm-based estimate 00056 * of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D) 00057 * and (B,E), using ZLACON. 00058 * 00059 * If IJOB >= 1, ZTGSYL computes a Frobenius norm-based estimate of 00060 * Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the 00061 * reciprocal of the smallest singular value of Z. 00062 * 00063 * This is a level-3 BLAS algorithm. 00064 * 00065 * Arguments 00066 * ========= 00067 * 00068 * TRANS (input) CHARACTER*1 00069 * = 'N': solve the generalized sylvester equation (1). 00070 * = 'C': solve the "conjugate transposed" system (3). 00071 * 00072 * IJOB (input) INTEGER 00073 * Specifies what kind of functionality to be performed. 00074 * =0: solve (1) only. 00075 * =1: The functionality of 0 and 3. 00076 * =2: The functionality of 0 and 4. 00077 * =3: Only an estimate of Dif[(A,D), (B,E)] is computed. 00078 * (look ahead strategy is used). 00079 * =4: Only an estimate of Dif[(A,D), (B,E)] is computed. 00080 * (ZGECON on sub-systems is used). 00081 * Not referenced if TRANS = 'C'. 00082 * 00083 * M (input) INTEGER 00084 * The order of the matrices A and D, and the row dimension of 00085 * the matrices C, F, R and L. 00086 * 00087 * N (input) INTEGER 00088 * The order of the matrices B and E, and the column dimension 00089 * of the matrices C, F, R and L. 00090 * 00091 * A (input) COMPLEX*16 array, dimension (LDA, M) 00092 * The upper triangular matrix A. 00093 * 00094 * LDA (input) INTEGER 00095 * The leading dimension of the array A. LDA >= max(1, M). 00096 * 00097 * B (input) COMPLEX*16 array, dimension (LDB, N) 00098 * The upper triangular matrix B. 00099 * 00100 * LDB (input) INTEGER 00101 * The leading dimension of the array B. LDB >= max(1, N). 00102 * 00103 * C (input/output) COMPLEX*16 array, dimension (LDC, N) 00104 * On entry, C contains the right-hand-side of the first matrix 00105 * equation in (1) or (3). 00106 * On exit, if IJOB = 0, 1 or 2, C has been overwritten by 00107 * the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, 00108 * the solution achieved during the computation of the 00109 * Dif-estimate. 00110 * 00111 * LDC (input) INTEGER 00112 * The leading dimension of the array C. LDC >= max(1, M). 00113 * 00114 * D (input) COMPLEX*16 array, dimension (LDD, M) 00115 * The upper triangular matrix D. 00116 * 00117 * LDD (input) INTEGER 00118 * The leading dimension of the array D. LDD >= max(1, M). 00119 * 00120 * E (input) COMPLEX*16 array, dimension (LDE, N) 00121 * The upper triangular matrix E. 00122 * 00123 * LDE (input) INTEGER 00124 * The leading dimension of the array E. LDE >= max(1, N). 00125 * 00126 * F (input/output) COMPLEX*16 array, dimension (LDF, N) 00127 * On entry, F contains the right-hand-side of the second matrix 00128 * equation in (1) or (3). 00129 * On exit, if IJOB = 0, 1 or 2, F has been overwritten by 00130 * the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, 00131 * the solution achieved during the computation of the 00132 * Dif-estimate. 00133 * 00134 * LDF (input) INTEGER 00135 * The leading dimension of the array F. LDF >= max(1, M). 00136 * 00137 * DIF (output) DOUBLE PRECISION 00138 * On exit DIF is the reciprocal of a lower bound of the 00139 * reciprocal of the Dif-function, i.e. DIF is an upper bound of 00140 * Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2). 00141 * IF IJOB = 0 or TRANS = 'C', DIF is not referenced. 00142 * 00143 * SCALE (output) DOUBLE PRECISION 00144 * On exit SCALE is the scaling factor in (1) or (3). 00145 * If 0 < SCALE < 1, C and F hold the solutions R and L, resp., 00146 * to a slightly perturbed system but the input matrices A, B, 00147 * D and E have not been changed. If SCALE = 0, R and L will 00148 * hold the solutions to the homogenious system with C = F = 0. 00149 * 00150 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00151 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00152 * 00153 * LWORK (input) INTEGER 00154 * The dimension of the array WORK. LWORK > = 1. 00155 * If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). 00156 * 00157 * If LWORK = -1, then a workspace query is assumed; the routine 00158 * only calculates the optimal size of the WORK array, returns 00159 * this value as the first entry of the WORK array, and no error 00160 * message related to LWORK is issued by XERBLA. 00161 * 00162 * IWORK (workspace) INTEGER array, dimension (M+N+2) 00163 * 00164 * INFO (output) INTEGER 00165 * =0: successful exit 00166 * <0: If INFO = -i, the i-th argument had an illegal value. 00167 * >0: (A, D) and (B, E) have common or very close 00168 * eigenvalues. 00169 * 00170 * Further Details 00171 * =============== 00172 * 00173 * Based on contributions by 00174 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00175 * Umea University, S-901 87 Umea, Sweden. 00176 * 00177 * [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software 00178 * for Solving the Generalized Sylvester Equation and Estimating the 00179 * Separation between Regular Matrix Pairs, Report UMINF - 93.23, 00180 * Department of Computing Science, Umea University, S-901 87 Umea, 00181 * Sweden, December 1993, Revised April 1994, Also as LAPACK Working 00182 * Note 75. To appear in ACM Trans. on Math. Software, Vol 22, 00183 * No 1, 1996. 00184 * 00185 * [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester 00186 * Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. 00187 * Appl., 15(4):1045-1060, 1994. 00188 * 00189 * [3] B. Kagstrom and L. Westin, Generalized Schur Methods with 00190 * Condition Estimators for Solving the Generalized Sylvester 00191 * Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, 00192 * July 1989, pp 745-751. 00193 * 00194 * ===================================================================== 00195 * Replaced various illegal calls to CCOPY by calls to CLASET. 00196 * Sven Hammarling, 1/5/02. 00197 * 00198 * .. Parameters .. 00199 DOUBLE PRECISION ZERO, ONE 00200 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00201 COMPLEX*16 CZERO 00202 PARAMETER ( CZERO = (0.0D+0, 0.0D+0) ) 00203 * .. 00204 * .. Local Scalars .. 00205 LOGICAL LQUERY, NOTRAN 00206 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K, 00207 $ LINFO, LWMIN, MB, NB, P, PQ, Q 00208 DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC 00209 * .. 00210 * .. External Functions .. 00211 LOGICAL LSAME 00212 INTEGER ILAENV 00213 EXTERNAL LSAME, ILAENV 00214 * .. 00215 * .. External Subroutines .. 00216 EXTERNAL XERBLA, ZGEMM, ZLACPY, ZLASET, ZSCAL, ZTGSY2 00217 * .. 00218 * .. Intrinsic Functions .. 00219 INTRINSIC DBLE, DCMPLX, MAX, SQRT 00220 * .. 00221 * .. Executable Statements .. 00222 * 00223 * Decode and test input parameters 00224 * 00225 INFO = 0 00226 NOTRAN = LSAME( TRANS, 'N' ) 00227 LQUERY = ( LWORK.EQ.-1 ) 00228 * 00229 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00230 INFO = -1 00231 ELSE IF( NOTRAN ) THEN 00232 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN 00233 INFO = -2 00234 END IF 00235 END IF 00236 IF( INFO.EQ.0 ) THEN 00237 IF( M.LE.0 ) THEN 00238 INFO = -3 00239 ELSE IF( N.LE.0 ) THEN 00240 INFO = -4 00241 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00242 INFO = -6 00243 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00244 INFO = -8 00245 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00246 INFO = -10 00247 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 00248 INFO = -12 00249 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 00250 INFO = -14 00251 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 00252 INFO = -16 00253 END IF 00254 END IF 00255 * 00256 IF( INFO.EQ.0 ) THEN 00257 IF( NOTRAN ) THEN 00258 IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN 00259 LWMIN = MAX( 1, 2*M*N ) 00260 ELSE 00261 LWMIN = 1 00262 END IF 00263 ELSE 00264 LWMIN = 1 00265 END IF 00266 WORK( 1 ) = LWMIN 00267 * 00268 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00269 INFO = -20 00270 END IF 00271 END IF 00272 * 00273 IF( INFO.NE.0 ) THEN 00274 CALL XERBLA( 'ZTGSYL', -INFO ) 00275 RETURN 00276 ELSE IF( LQUERY ) THEN 00277 RETURN 00278 END IF 00279 * 00280 * Quick return if possible 00281 * 00282 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00283 SCALE = 1 00284 IF( NOTRAN ) THEN 00285 IF( IJOB.NE.0 ) THEN 00286 DIF = 0 00287 END IF 00288 END IF 00289 RETURN 00290 END IF 00291 * 00292 * Determine optimal block sizes MB and NB 00293 * 00294 MB = ILAENV( 2, 'ZTGSYL', TRANS, M, N, -1, -1 ) 00295 NB = ILAENV( 5, 'ZTGSYL', TRANS, M, N, -1, -1 ) 00296 * 00297 ISOLVE = 1 00298 IFUNC = 0 00299 IF( NOTRAN ) THEN 00300 IF( IJOB.GE.3 ) THEN 00301 IFUNC = IJOB - 2 00302 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00303 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 00304 ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN 00305 ISOLVE = 2 00306 END IF 00307 END IF 00308 * 00309 IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) ) 00310 $ THEN 00311 * 00312 * Use unblocked Level 2 solver 00313 * 00314 DO 30 IROUND = 1, ISOLVE 00315 * 00316 SCALE = ONE 00317 DSCALE = ZERO 00318 DSUM = ONE 00319 PQ = M*N 00320 CALL ZTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D, 00321 $ LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE, 00322 $ INFO ) 00323 IF( DSCALE.NE.ZERO ) THEN 00324 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN 00325 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 00326 ELSE 00327 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) 00328 END IF 00329 END IF 00330 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN 00331 IF( NOTRAN ) THEN 00332 IFUNC = IJOB 00333 END IF 00334 SCALE2 = SCALE 00335 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M ) 00336 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 00337 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00338 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 00339 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 00340 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC ) 00341 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) 00342 SCALE = SCALE2 00343 END IF 00344 30 CONTINUE 00345 * 00346 RETURN 00347 * 00348 END IF 00349 * 00350 * Determine block structure of A 00351 * 00352 P = 0 00353 I = 1 00354 40 CONTINUE 00355 IF( I.GT.M ) 00356 $ GO TO 50 00357 P = P + 1 00358 IWORK( P ) = I 00359 I = I + MB 00360 IF( I.GE.M ) 00361 $ GO TO 50 00362 GO TO 40 00363 50 CONTINUE 00364 IWORK( P+1 ) = M + 1 00365 IF( IWORK( P ).EQ.IWORK( P+1 ) ) 00366 $ P = P - 1 00367 * 00368 * Determine block structure of B 00369 * 00370 Q = P + 1 00371 J = 1 00372 60 CONTINUE 00373 IF( J.GT.N ) 00374 $ GO TO 70 00375 * 00376 Q = Q + 1 00377 IWORK( Q ) = J 00378 J = J + NB 00379 IF( J.GE.N ) 00380 $ GO TO 70 00381 GO TO 60 00382 * 00383 70 CONTINUE 00384 IWORK( Q+1 ) = N + 1 00385 IF( IWORK( Q ).EQ.IWORK( Q+1 ) ) 00386 $ Q = Q - 1 00387 * 00388 IF( NOTRAN ) THEN 00389 DO 150 IROUND = 1, ISOLVE 00390 * 00391 * Solve (I, J) - subsystem 00392 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 00393 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 00394 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q 00395 * 00396 PQ = 0 00397 SCALE = ONE 00398 DSCALE = ZERO 00399 DSUM = ONE 00400 DO 130 J = P + 2, Q 00401 JS = IWORK( J ) 00402 JE = IWORK( J+1 ) - 1 00403 NB = JE - JS + 1 00404 DO 120 I = P, 1, -1 00405 IS = IWORK( I ) 00406 IE = IWORK( I+1 ) - 1 00407 MB = IE - IS + 1 00408 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA, 00409 $ B( JS, JS ), LDB, C( IS, JS ), LDC, 00410 $ D( IS, IS ), LDD, E( JS, JS ), LDE, 00411 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE, 00412 $ LINFO ) 00413 IF( LINFO.GT.0 ) 00414 $ INFO = LINFO 00415 PQ = PQ + MB*NB 00416 IF( SCALOC.NE.ONE ) THEN 00417 DO 80 K = 1, JS - 1 00418 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00419 $ C( 1, K ), 1 ) 00420 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00421 $ F( 1, K ), 1 ) 00422 80 CONTINUE 00423 DO 90 K = JS, JE 00424 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 00425 $ C( 1, K ), 1 ) 00426 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 00427 $ F( 1, K ), 1 ) 00428 90 CONTINUE 00429 DO 100 K = JS, JE 00430 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 00431 $ C( IE+1, K ), 1 ) 00432 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 00433 $ F( IE+1, K ), 1 ) 00434 100 CONTINUE 00435 DO 110 K = JE + 1, N 00436 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00437 $ C( 1, K ), 1 ) 00438 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), 00439 $ F( 1, K ), 1 ) 00440 110 CONTINUE 00441 SCALE = SCALE*SCALOC 00442 END IF 00443 * 00444 * Substitute R(I,J) and L(I,J) into remaining equation. 00445 * 00446 IF( I.GT.1 ) THEN 00447 CALL ZGEMM( 'N', 'N', IS-1, NB, MB, 00448 $ DCMPLX( -ONE, ZERO ), A( 1, IS ), LDA, 00449 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 00450 $ C( 1, JS ), LDC ) 00451 CALL ZGEMM( 'N', 'N', IS-1, NB, MB, 00452 $ DCMPLX( -ONE, ZERO ), D( 1, IS ), LDD, 00453 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 00454 $ F( 1, JS ), LDF ) 00455 END IF 00456 IF( J.LT.Q ) THEN 00457 CALL ZGEMM( 'N', 'N', MB, N-JE, NB, 00458 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00459 $ B( JS, JE+1 ), LDB, 00460 $ DCMPLX( ONE, ZERO ), C( IS, JE+1 ), 00461 $ LDC ) 00462 CALL ZGEMM( 'N', 'N', MB, N-JE, NB, 00463 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00464 $ E( JS, JE+1 ), LDE, 00465 $ DCMPLX( ONE, ZERO ), F( IS, JE+1 ), 00466 $ LDF ) 00467 END IF 00468 120 CONTINUE 00469 130 CONTINUE 00470 IF( DSCALE.NE.ZERO ) THEN 00471 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN 00472 DIF = SQRT( DBLE( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 00473 ELSE 00474 DIF = SQRT( DBLE( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) 00475 END IF 00476 END IF 00477 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN 00478 IF( NOTRAN ) THEN 00479 IFUNC = IJOB 00480 END IF 00481 SCALE2 = SCALE 00482 CALL ZLACPY( 'F', M, N, C, LDC, WORK, M ) 00483 CALL ZLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 00484 CALL ZLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00485 CALL ZLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 00486 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 00487 CALL ZLACPY( 'F', M, N, WORK, M, C, LDC ) 00488 CALL ZLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) 00489 SCALE = SCALE2 00490 END IF 00491 150 CONTINUE 00492 ELSE 00493 * 00494 * Solve transposed (I, J)-subsystem 00495 * A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J) 00496 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 00497 * for I = 1,2,..., P; J = Q, Q-1,..., 1 00498 * 00499 SCALE = ONE 00500 DO 210 I = 1, P 00501 IS = IWORK( I ) 00502 IE = IWORK( I+1 ) - 1 00503 MB = IE - IS + 1 00504 DO 200 J = Q, P + 2, -1 00505 JS = IWORK( J ) 00506 JE = IWORK( J+1 ) - 1 00507 NB = JE - JS + 1 00508 CALL ZTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA, 00509 $ B( JS, JS ), LDB, C( IS, JS ), LDC, 00510 $ D( IS, IS ), LDD, E( JS, JS ), LDE, 00511 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE, 00512 $ LINFO ) 00513 IF( LINFO.GT.0 ) 00514 $ INFO = LINFO 00515 IF( SCALOC.NE.ONE ) THEN 00516 DO 160 K = 1, JS - 1 00517 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 00518 $ 1 ) 00519 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 00520 $ 1 ) 00521 160 CONTINUE 00522 DO 170 K = JS, JE 00523 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 00524 $ C( 1, K ), 1 ) 00525 CALL ZSCAL( IS-1, DCMPLX( SCALOC, ZERO ), 00526 $ F( 1, K ), 1 ) 00527 170 CONTINUE 00528 DO 180 K = JS, JE 00529 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 00530 $ C( IE+1, K ), 1 ) 00531 CALL ZSCAL( M-IE, DCMPLX( SCALOC, ZERO ), 00532 $ F( IE+1, K ), 1 ) 00533 180 CONTINUE 00534 DO 190 K = JE + 1, N 00535 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ), 00536 $ 1 ) 00537 CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ), 00538 $ 1 ) 00539 190 CONTINUE 00540 SCALE = SCALE*SCALOC 00541 END IF 00542 * 00543 * Substitute R(I,J) and L(I,J) into remaining equation. 00544 * 00545 IF( J.GT.P+2 ) THEN 00546 CALL ZGEMM( 'N', 'C', MB, JS-1, NB, 00547 $ DCMPLX( ONE, ZERO ), C( IS, JS ), LDC, 00548 $ B( 1, JS ), LDB, DCMPLX( ONE, ZERO ), 00549 $ F( IS, 1 ), LDF ) 00550 CALL ZGEMM( 'N', 'C', MB, JS-1, NB, 00551 $ DCMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00552 $ E( 1, JS ), LDE, DCMPLX( ONE, ZERO ), 00553 $ F( IS, 1 ), LDF ) 00554 END IF 00555 IF( I.LT.P ) THEN 00556 CALL ZGEMM( 'C', 'N', M-IE, NB, MB, 00557 $ DCMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA, 00558 $ C( IS, JS ), LDC, DCMPLX( ONE, ZERO ), 00559 $ C( IE+1, JS ), LDC ) 00560 CALL ZGEMM( 'C', 'N', M-IE, NB, MB, 00561 $ DCMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD, 00562 $ F( IS, JS ), LDF, DCMPLX( ONE, ZERO ), 00563 $ C( IE+1, JS ), LDC ) 00564 END IF 00565 200 CONTINUE 00566 210 CONTINUE 00567 END IF 00568 * 00569 WORK( 1 ) = LWMIN 00570 * 00571 RETURN 00572 * 00573 * End of ZTGSYL 00574 * 00575 END
1.7.3 
