 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CTGSYL( 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.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER TRANS 00012 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, 00013 $ LWORK, M, N 00014 REAL DIF, SCALE 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IWORK( * ) 00018 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ), 00019 $ D( LDD, * ), E( LDE, * ), F( LDF, * ), 00020 $ WORK( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * CTGSYL 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**H, Im) ] (2) 00043 * [ kron(In, D) -kron(E**H, Im) ], 00044 * 00045 * Here Ix is the identity matrix of size x and X**H 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**H *y = scale*b 00050 * is solved for, which is equivalent to solve for R and L in 00051 * 00052 * A**H * R + D**H * L = scale * C (3) 00053 * R * B**H + L * E**H = 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 CLACON. 00058 * 00059 * If IJOB >= 1, CTGSYL 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 * (CGECON 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 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 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 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 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 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 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) REAL 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) REAL 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 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 REAL ZERO, ONE 00200 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00201 COMPLEX CZERO 00202 PARAMETER ( CZERO = (0.0E+0, 0.0E+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 REAL DSCALE, DSUM, SCALE2, SCALOC 00209 * .. 00210 * .. External Functions .. 00211 LOGICAL LSAME 00212 INTEGER ILAENV 00213 EXTERNAL LSAME, ILAENV 00214 * .. 00215 * .. External Subroutines .. 00216 EXTERNAL CGEMM, CLACPY, CLASET, CSCAL, CTGSY2, XERBLA 00217 * .. 00218 * .. Intrinsic Functions .. 00219 INTRINSIC CMPLX, MAX, REAL, 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( 'CTGSYL', -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, 'CTGSYL', TRANS, M, N, -1, -1 ) 00295 NB = ILAENV( 5, 'CTGSYL', 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 CLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00303 CALL CLASET( '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 CTGSY2( 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( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 00326 ELSE 00327 DIF = SQRT( REAL( 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 CLACPY( 'F', M, N, C, LDC, WORK, M ) 00336 CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 00337 CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00338 CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 00339 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 00340 CALL CLACPY( 'F', M, N, WORK, M, C, LDC ) 00341 CALL CLACPY( '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 CTGSY2( 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 CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00419 $ 1 ) 00420 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00421 $ 1 ) 00422 80 CONTINUE 00423 DO 90 K = JS, JE 00424 CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), 00425 $ C( 1, K ), 1 ) 00426 CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), 00427 $ F( 1, K ), 1 ) 00428 90 CONTINUE 00429 DO 100 K = JS, JE 00430 CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ), 00431 $ C( IE+1, K ), 1 ) 00432 CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ), 00433 $ F( IE+1, K ), 1 ) 00434 100 CONTINUE 00435 DO 110 K = JE + 1, N 00436 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00437 $ 1 ) 00438 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00439 $ 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 CGEMM( 'N', 'N', IS-1, NB, MB, 00448 $ CMPLX( -ONE, ZERO ), A( 1, IS ), LDA, 00449 $ C( IS, JS ), LDC, CMPLX( ONE, ZERO ), 00450 $ C( 1, JS ), LDC ) 00451 CALL CGEMM( 'N', 'N', IS-1, NB, MB, 00452 $ CMPLX( -ONE, ZERO ), D( 1, IS ), LDD, 00453 $ C( IS, JS ), LDC, CMPLX( ONE, ZERO ), 00454 $ F( 1, JS ), LDF ) 00455 END IF 00456 IF( J.LT.Q ) THEN 00457 CALL CGEMM( 'N', 'N', MB, N-JE, NB, 00458 $ CMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00459 $ B( JS, JE+1 ), LDB, CMPLX( ONE, ZERO ), 00460 $ C( IS, JE+1 ), LDC ) 00461 CALL CGEMM( 'N', 'N', MB, N-JE, NB, 00462 $ CMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00463 $ E( JS, JE+1 ), LDE, CMPLX( ONE, ZERO ), 00464 $ F( IS, JE+1 ), LDF ) 00465 END IF 00466 120 CONTINUE 00467 130 CONTINUE 00468 IF( DSCALE.NE.ZERO ) THEN 00469 IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN 00470 DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) ) 00471 ELSE 00472 DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) ) 00473 END IF 00474 END IF 00475 IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN 00476 IF( NOTRAN ) THEN 00477 IFUNC = IJOB 00478 END IF 00479 SCALE2 = SCALE 00480 CALL CLACPY( 'F', M, N, C, LDC, WORK, M ) 00481 CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M ) 00482 CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC ) 00483 CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF ) 00484 ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN 00485 CALL CLACPY( 'F', M, N, WORK, M, C, LDC ) 00486 CALL CLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF ) 00487 SCALE = SCALE2 00488 END IF 00489 150 CONTINUE 00490 ELSE 00491 * 00492 * Solve transposed (I, J)-subsystem 00493 * A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J) 00494 * R(I, J) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 00495 * for I = 1,2,..., P; J = Q, Q-1,..., 1 00496 * 00497 SCALE = ONE 00498 DO 210 I = 1, P 00499 IS = IWORK( I ) 00500 IE = IWORK( I+1 ) - 1 00501 MB = IE - IS + 1 00502 DO 200 J = Q, P + 2, -1 00503 JS = IWORK( J ) 00504 JE = IWORK( J+1 ) - 1 00505 NB = JE - JS + 1 00506 CALL CTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA, 00507 $ B( JS, JS ), LDB, C( IS, JS ), LDC, 00508 $ D( IS, IS ), LDD, E( JS, JS ), LDE, 00509 $ F( IS, JS ), LDF, SCALOC, DSUM, DSCALE, 00510 $ LINFO ) 00511 IF( LINFO.GT.0 ) 00512 $ INFO = LINFO 00513 IF( SCALOC.NE.ONE ) THEN 00514 DO 160 K = 1, JS - 1 00515 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00516 $ 1 ) 00517 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00518 $ 1 ) 00519 160 CONTINUE 00520 DO 170 K = JS, JE 00521 CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), C( 1, K ), 00522 $ 1 ) 00523 CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), F( 1, K ), 00524 $ 1 ) 00525 170 CONTINUE 00526 DO 180 K = JS, JE 00527 CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ), 00528 $ C( IE+1, K ), 1 ) 00529 CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ), 00530 $ F( IE+1, K ), 1 ) 00531 180 CONTINUE 00532 DO 190 K = JE + 1, N 00533 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ), 00534 $ 1 ) 00535 CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ), 00536 $ 1 ) 00537 190 CONTINUE 00538 SCALE = SCALE*SCALOC 00539 END IF 00540 * 00541 * Substitute R(I,J) and L(I,J) into remaining equation. 00542 * 00543 IF( J.GT.P+2 ) THEN 00544 CALL CGEMM( 'N', 'C', MB, JS-1, NB, 00545 $ CMPLX( ONE, ZERO ), C( IS, JS ), LDC, 00546 $ B( 1, JS ), LDB, CMPLX( ONE, ZERO ), 00547 $ F( IS, 1 ), LDF ) 00548 CALL CGEMM( 'N', 'C', MB, JS-1, NB, 00549 $ CMPLX( ONE, ZERO ), F( IS, JS ), LDF, 00550 $ E( 1, JS ), LDE, CMPLX( ONE, ZERO ), 00551 $ F( IS, 1 ), LDF ) 00552 END IF 00553 IF( I.LT.P ) THEN 00554 CALL CGEMM( 'C', 'N', M-IE, NB, MB, 00555 $ CMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA, 00556 $ C( IS, JS ), LDC, CMPLX( ONE, ZERO ), 00557 $ C( IE+1, JS ), LDC ) 00558 CALL CGEMM( 'C', 'N', M-IE, NB, MB, 00559 $ CMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD, 00560 $ F( IS, JS ), LDF, CMPLX( ONE, ZERO ), 00561 $ C( IE+1, JS ), LDC ) 00562 END IF 00563 200 CONTINUE 00564 210 CONTINUE 00565 END IF 00566 * 00567 WORK( 1 ) = LWMIN 00568 * 00569 RETURN 00570 * 00571 * End of CTGSYL 00572 * 00573 END
1.7.3 
