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