 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, 00002 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, 00003 $ IWORK, PQ, INFO ) 00004 * 00005 * -- LAPACK auxiliary 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, M, N, 00013 $ PQ 00014 DOUBLE PRECISION RDSCAL, RDSUM, 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 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * DTGSY2 solves the generalized Sylvester equation: 00026 * 00027 * A * R - L * B = scale * C (1) 00028 * D * R - L * E = scale * F, 00029 * 00030 * using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices, 00031 * (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, 00032 * N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E) 00033 * must be in generalized Schur canonical form, i.e. A, B are upper 00034 * quasi triangular and D, E are upper triangular. The solution (R, L) 00035 * overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor 00036 * chosen to avoid overflow. 00037 * 00038 * In matrix notation solving equation (1) corresponds to solve 00039 * Z*x = scale*b, where Z is defined as 00040 * 00041 * Z = [ kron(In, A) -kron(B**T, Im) ] (2) 00042 * [ kron(In, D) -kron(E**T, Im) ], 00043 * 00044 * Ik is the identity matrix of size k and X**T is the transpose of X. 00045 * kron(X, Y) is the Kronecker product between the matrices X and Y. 00046 * In the process of solving (1), we solve a number of such systems 00047 * where Dim(In), Dim(In) = 1 or 2. 00048 * 00049 * If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y, 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 is used to compute an estimate of Dif[(A, D), (B, E)] = 00056 * sigma_min(Z) using reverse communicaton with DLACON. 00057 * 00058 * DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL 00059 * of an upper bound on the separation between to matrix pairs. Then 00060 * the input (A, D), (B, E) are sub-pencils of the matrix pair in 00061 * DTGSYL. See DTGSYL for details. 00062 * 00063 * Arguments 00064 * ========= 00065 * 00066 * TRANS (input) CHARACTER*1 00067 * = 'N', solve the generalized Sylvester equation (1). 00068 * = 'T': solve the 'transposed' system (3). 00069 * 00070 * IJOB (input) INTEGER 00071 * Specifies what kind of functionality to be performed. 00072 * = 0: solve (1) only. 00073 * = 1: A contribution from this subsystem to a Frobenius 00074 * norm-based estimate of the separation between two matrix 00075 * pairs is computed. (look ahead strategy is used). 00076 * = 2: A contribution from this subsystem to a Frobenius 00077 * norm-based estimate of the separation between two matrix 00078 * pairs is computed. (DGECON on sub-systems is used.) 00079 * Not referenced if TRANS = 'T'. 00080 * 00081 * M (input) INTEGER 00082 * On entry, M specifies the order of A and D, and the row 00083 * dimension of C, F, R and L. 00084 * 00085 * N (input) INTEGER 00086 * On entry, N specifies the order of B and E, and the column 00087 * dimension of C, F, R and L. 00088 * 00089 * A (input) DOUBLE PRECISION array, dimension (LDA, M) 00090 * On entry, A contains an upper quasi triangular matrix. 00091 * 00092 * LDA (input) INTEGER 00093 * The leading dimension of the matrix A. LDA >= max(1, M). 00094 * 00095 * B (input) DOUBLE PRECISION array, dimension (LDB, N) 00096 * On entry, B contains an upper quasi triangular matrix. 00097 * 00098 * LDB (input) INTEGER 00099 * The leading dimension of the matrix B. LDB >= max(1, N). 00100 * 00101 * C (input/output) DOUBLE PRECISION array, dimension (LDC, N) 00102 * On entry, C contains the right-hand-side of the first matrix 00103 * equation in (1). 00104 * On exit, if IJOB = 0, C has been overwritten by the 00105 * solution R. 00106 * 00107 * LDC (input) INTEGER 00108 * The leading dimension of the matrix C. LDC >= max(1, M). 00109 * 00110 * D (input) DOUBLE PRECISION array, dimension (LDD, M) 00111 * On entry, D contains an upper triangular matrix. 00112 * 00113 * LDD (input) INTEGER 00114 * The leading dimension of the matrix D. LDD >= max(1, M). 00115 * 00116 * E (input) DOUBLE PRECISION array, dimension (LDE, N) 00117 * On entry, E contains an upper triangular matrix. 00118 * 00119 * LDE (input) INTEGER 00120 * The leading dimension of the matrix E. LDE >= max(1, N). 00121 * 00122 * F (input/output) DOUBLE PRECISION array, dimension (LDF, N) 00123 * On entry, F contains the right-hand-side of the second matrix 00124 * equation in (1). 00125 * On exit, if IJOB = 0, F has been overwritten by the 00126 * solution L. 00127 * 00128 * LDF (input) INTEGER 00129 * The leading dimension of the matrix F. LDF >= max(1, M). 00130 * 00131 * SCALE (output) DOUBLE PRECISION 00132 * On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions 00133 * R and L (C and F on entry) will hold the solutions to a 00134 * slightly perturbed system but the input matrices A, B, D and 00135 * E have not been changed. If SCALE = 0, R and L will hold the 00136 * solutions to the homogeneous system with C = F = 0. Normally, 00137 * SCALE = 1. 00138 * 00139 * RDSUM (input/output) DOUBLE PRECISION 00140 * On entry, the sum of squares of computed contributions to 00141 * the Dif-estimate under computation by DTGSYL, where the 00142 * scaling factor RDSCAL (see below) has been factored out. 00143 * On exit, the corresponding sum of squares updated with the 00144 * contributions from the current sub-system. 00145 * If TRANS = 'T' RDSUM is not touched. 00146 * NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL. 00147 * 00148 * RDSCAL (input/output) DOUBLE PRECISION 00149 * On entry, scaling factor used to prevent overflow in RDSUM. 00150 * On exit, RDSCAL is updated w.r.t. the current contributions 00151 * in RDSUM. 00152 * If TRANS = 'T', RDSCAL is not touched. 00153 * NOTE: RDSCAL only makes sense when DTGSY2 is called by 00154 * DTGSYL. 00155 * 00156 * IWORK (workspace) INTEGER array, dimension (M+N+2) 00157 * 00158 * PQ (output) INTEGER 00159 * On exit, the number of subsystems (of size 2-by-2, 4-by-4 and 00160 * 8-by-8) solved by this routine. 00161 * 00162 * INFO (output) INTEGER 00163 * On exit, if INFO is set to 00164 * =0: Successful exit 00165 * <0: If INFO = -i, the i-th argument had an illegal value. 00166 * >0: The matrix pairs (A, D) and (B, E) have common or very 00167 * close eigenvalues. 00168 * 00169 * Further Details 00170 * =============== 00171 * 00172 * Based on contributions by 00173 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00174 * Umea University, S-901 87 Umea, Sweden. 00175 * 00176 * ===================================================================== 00177 * Replaced various illegal calls to DCOPY by calls to DLASET. 00178 * Sven Hammarling, 27/5/02. 00179 * 00180 * .. Parameters .. 00181 INTEGER LDZ 00182 PARAMETER ( LDZ = 8 ) 00183 DOUBLE PRECISION ZERO, ONE 00184 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00185 * .. 00186 * .. Local Scalars .. 00187 LOGICAL NOTRAN 00188 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1, 00189 $ K, MB, NB, P, Q, ZDIM 00190 DOUBLE PRECISION ALPHA, SCALOC 00191 * .. 00192 * .. Local Arrays .. 00193 INTEGER IPIV( LDZ ), JPIV( LDZ ) 00194 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ ) 00195 * .. 00196 * .. External Functions .. 00197 LOGICAL LSAME 00198 EXTERNAL LSAME 00199 * .. 00200 * .. External Subroutines .. 00201 EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2, 00202 $ DGETC2, DLASET, DLATDF, DSCAL, XERBLA 00203 * .. 00204 * .. Intrinsic Functions .. 00205 INTRINSIC MAX 00206 * .. 00207 * .. Executable Statements .. 00208 * 00209 * Decode and test input parameters 00210 * 00211 INFO = 0 00212 IERR = 0 00213 NOTRAN = LSAME( TRANS, 'N' ) 00214 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 00215 INFO = -1 00216 ELSE IF( NOTRAN ) THEN 00217 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN 00218 INFO = -2 00219 END IF 00220 END IF 00221 IF( INFO.EQ.0 ) THEN 00222 IF( M.LE.0 ) THEN 00223 INFO = -3 00224 ELSE IF( N.LE.0 ) THEN 00225 INFO = -4 00226 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00227 INFO = -5 00228 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00229 INFO = -8 00230 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00231 INFO = -10 00232 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN 00233 INFO = -12 00234 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN 00235 INFO = -14 00236 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN 00237 INFO = -16 00238 END IF 00239 END IF 00240 IF( INFO.NE.0 ) THEN 00241 CALL XERBLA( 'DTGSY2', -INFO ) 00242 RETURN 00243 END IF 00244 * 00245 * Determine block structure of A 00246 * 00247 PQ = 0 00248 P = 0 00249 I = 1 00250 10 CONTINUE 00251 IF( I.GT.M ) 00252 $ GO TO 20 00253 P = P + 1 00254 IWORK( P ) = I 00255 IF( I.EQ.M ) 00256 $ GO TO 20 00257 IF( A( I+1, I ).NE.ZERO ) THEN 00258 I = I + 2 00259 ELSE 00260 I = I + 1 00261 END IF 00262 GO TO 10 00263 20 CONTINUE 00264 IWORK( P+1 ) = M + 1 00265 * 00266 * Determine block structure of B 00267 * 00268 Q = P + 1 00269 J = 1 00270 30 CONTINUE 00271 IF( J.GT.N ) 00272 $ GO TO 40 00273 Q = Q + 1 00274 IWORK( Q ) = J 00275 IF( J.EQ.N ) 00276 $ GO TO 40 00277 IF( B( J+1, J ).NE.ZERO ) THEN 00278 J = J + 2 00279 ELSE 00280 J = J + 1 00281 END IF 00282 GO TO 30 00283 40 CONTINUE 00284 IWORK( Q+1 ) = N + 1 00285 PQ = P*( Q-P-1 ) 00286 * 00287 IF( NOTRAN ) THEN 00288 * 00289 * Solve (I, J) - subsystem 00290 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) 00291 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) 00292 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q 00293 * 00294 SCALE = ONE 00295 SCALOC = ONE 00296 DO 120 J = P + 2, Q 00297 JS = IWORK( J ) 00298 JSP1 = JS + 1 00299 JE = IWORK( J+1 ) - 1 00300 NB = JE - JS + 1 00301 DO 110 I = P, 1, -1 00302 * 00303 IS = IWORK( I ) 00304 ISP1 = IS + 1 00305 IE = IWORK( I+1 ) - 1 00306 MB = IE - IS + 1 00307 ZDIM = MB*NB*2 00308 * 00309 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 00310 * 00311 * Build a 2-by-2 system Z * x = RHS 00312 * 00313 Z( 1, 1 ) = A( IS, IS ) 00314 Z( 2, 1 ) = D( IS, IS ) 00315 Z( 1, 2 ) = -B( JS, JS ) 00316 Z( 2, 2 ) = -E( JS, JS ) 00317 * 00318 * Set up right hand side(s) 00319 * 00320 RHS( 1 ) = C( IS, JS ) 00321 RHS( 2 ) = F( IS, JS ) 00322 * 00323 * Solve Z * x = RHS 00324 * 00325 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00326 IF( IERR.GT.0 ) 00327 $ INFO = IERR 00328 * 00329 IF( IJOB.EQ.0 ) THEN 00330 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 00331 $ SCALOC ) 00332 IF( SCALOC.NE.ONE ) THEN 00333 DO 50 K = 1, N 00334 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00335 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00336 50 CONTINUE 00337 SCALE = SCALE*SCALOC 00338 END IF 00339 ELSE 00340 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 00341 $ RDSCAL, IPIV, JPIV ) 00342 END IF 00343 * 00344 * Unpack solution vector(s) 00345 * 00346 C( IS, JS ) = RHS( 1 ) 00347 F( IS, JS ) = RHS( 2 ) 00348 * 00349 * Substitute R(I, J) and L(I, J) into remaining 00350 * equation. 00351 * 00352 IF( I.GT.1 ) THEN 00353 ALPHA = -RHS( 1 ) 00354 CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ), 00355 $ 1 ) 00356 CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ), 00357 $ 1 ) 00358 END IF 00359 IF( J.LT.Q ) THEN 00360 CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB, 00361 $ C( IS, JE+1 ), LDC ) 00362 CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE, 00363 $ F( IS, JE+1 ), LDF ) 00364 END IF 00365 * 00366 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 00367 * 00368 * Build a 4-by-4 system Z * x = RHS 00369 * 00370 Z( 1, 1 ) = A( IS, IS ) 00371 Z( 2, 1 ) = ZERO 00372 Z( 3, 1 ) = D( IS, IS ) 00373 Z( 4, 1 ) = ZERO 00374 * 00375 Z( 1, 2 ) = ZERO 00376 Z( 2, 2 ) = A( IS, IS ) 00377 Z( 3, 2 ) = ZERO 00378 Z( 4, 2 ) = D( IS, IS ) 00379 * 00380 Z( 1, 3 ) = -B( JS, JS ) 00381 Z( 2, 3 ) = -B( JS, JSP1 ) 00382 Z( 3, 3 ) = -E( JS, JS ) 00383 Z( 4, 3 ) = -E( JS, JSP1 ) 00384 * 00385 Z( 1, 4 ) = -B( JSP1, JS ) 00386 Z( 2, 4 ) = -B( JSP1, JSP1 ) 00387 Z( 3, 4 ) = ZERO 00388 Z( 4, 4 ) = -E( JSP1, JSP1 ) 00389 * 00390 * Set up right hand side(s) 00391 * 00392 RHS( 1 ) = C( IS, JS ) 00393 RHS( 2 ) = C( IS, JSP1 ) 00394 RHS( 3 ) = F( IS, JS ) 00395 RHS( 4 ) = F( IS, JSP1 ) 00396 * 00397 * Solve Z * x = RHS 00398 * 00399 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00400 IF( IERR.GT.0 ) 00401 $ INFO = IERR 00402 * 00403 IF( IJOB.EQ.0 ) THEN 00404 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 00405 $ SCALOC ) 00406 IF( SCALOC.NE.ONE ) THEN 00407 DO 60 K = 1, N 00408 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00409 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00410 60 CONTINUE 00411 SCALE = SCALE*SCALOC 00412 END IF 00413 ELSE 00414 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 00415 $ RDSCAL, IPIV, JPIV ) 00416 END IF 00417 * 00418 * Unpack solution vector(s) 00419 * 00420 C( IS, JS ) = RHS( 1 ) 00421 C( IS, JSP1 ) = RHS( 2 ) 00422 F( IS, JS ) = RHS( 3 ) 00423 F( IS, JSP1 ) = RHS( 4 ) 00424 * 00425 * Substitute R(I, J) and L(I, J) into remaining 00426 * equation. 00427 * 00428 IF( I.GT.1 ) THEN 00429 CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ), 00430 $ 1, C( 1, JS ), LDC ) 00431 CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ), 00432 $ 1, F( 1, JS ), LDF ) 00433 END IF 00434 IF( J.LT.Q ) THEN 00435 CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB, 00436 $ C( IS, JE+1 ), LDC ) 00437 CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE, 00438 $ F( IS, JE+1 ), LDF ) 00439 CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB, 00440 $ C( IS, JE+1 ), LDC ) 00441 CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE, 00442 $ F( IS, JE+1 ), LDF ) 00443 END IF 00444 * 00445 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 00446 * 00447 * Build a 4-by-4 system Z * x = RHS 00448 * 00449 Z( 1, 1 ) = A( IS, IS ) 00450 Z( 2, 1 ) = A( ISP1, IS ) 00451 Z( 3, 1 ) = D( IS, IS ) 00452 Z( 4, 1 ) = ZERO 00453 * 00454 Z( 1, 2 ) = A( IS, ISP1 ) 00455 Z( 2, 2 ) = A( ISP1, ISP1 ) 00456 Z( 3, 2 ) = D( IS, ISP1 ) 00457 Z( 4, 2 ) = D( ISP1, ISP1 ) 00458 * 00459 Z( 1, 3 ) = -B( JS, JS ) 00460 Z( 2, 3 ) = ZERO 00461 Z( 3, 3 ) = -E( JS, JS ) 00462 Z( 4, 3 ) = ZERO 00463 * 00464 Z( 1, 4 ) = ZERO 00465 Z( 2, 4 ) = -B( JS, JS ) 00466 Z( 3, 4 ) = ZERO 00467 Z( 4, 4 ) = -E( JS, JS ) 00468 * 00469 * Set up right hand side(s) 00470 * 00471 RHS( 1 ) = C( IS, JS ) 00472 RHS( 2 ) = C( ISP1, JS ) 00473 RHS( 3 ) = F( IS, JS ) 00474 RHS( 4 ) = F( ISP1, JS ) 00475 * 00476 * Solve Z * x = RHS 00477 * 00478 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00479 IF( IERR.GT.0 ) 00480 $ INFO = IERR 00481 IF( IJOB.EQ.0 ) THEN 00482 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 00483 $ SCALOC ) 00484 IF( SCALOC.NE.ONE ) THEN 00485 DO 70 K = 1, N 00486 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00487 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00488 70 CONTINUE 00489 SCALE = SCALE*SCALOC 00490 END IF 00491 ELSE 00492 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 00493 $ RDSCAL, IPIV, JPIV ) 00494 END IF 00495 * 00496 * Unpack solution vector(s) 00497 * 00498 C( IS, JS ) = RHS( 1 ) 00499 C( ISP1, JS ) = RHS( 2 ) 00500 F( IS, JS ) = RHS( 3 ) 00501 F( ISP1, JS ) = RHS( 4 ) 00502 * 00503 * Substitute R(I, J) and L(I, J) into remaining 00504 * equation. 00505 * 00506 IF( I.GT.1 ) THEN 00507 CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA, 00508 $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 ) 00509 CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD, 00510 $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 ) 00511 END IF 00512 IF( J.LT.Q ) THEN 00513 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1, 00514 $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC ) 00515 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1, 00516 $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF ) 00517 END IF 00518 * 00519 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 00520 * 00521 * Build an 8-by-8 system Z * x = RHS 00522 * 00523 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 00524 * 00525 Z( 1, 1 ) = A( IS, IS ) 00526 Z( 2, 1 ) = A( ISP1, IS ) 00527 Z( 5, 1 ) = D( IS, IS ) 00528 * 00529 Z( 1, 2 ) = A( IS, ISP1 ) 00530 Z( 2, 2 ) = A( ISP1, ISP1 ) 00531 Z( 5, 2 ) = D( IS, ISP1 ) 00532 Z( 6, 2 ) = D( ISP1, ISP1 ) 00533 * 00534 Z( 3, 3 ) = A( IS, IS ) 00535 Z( 4, 3 ) = A( ISP1, IS ) 00536 Z( 7, 3 ) = D( IS, IS ) 00537 * 00538 Z( 3, 4 ) = A( IS, ISP1 ) 00539 Z( 4, 4 ) = A( ISP1, ISP1 ) 00540 Z( 7, 4 ) = D( IS, ISP1 ) 00541 Z( 8, 4 ) = D( ISP1, ISP1 ) 00542 * 00543 Z( 1, 5 ) = -B( JS, JS ) 00544 Z( 3, 5 ) = -B( JS, JSP1 ) 00545 Z( 5, 5 ) = -E( JS, JS ) 00546 Z( 7, 5 ) = -E( JS, JSP1 ) 00547 * 00548 Z( 2, 6 ) = -B( JS, JS ) 00549 Z( 4, 6 ) = -B( JS, JSP1 ) 00550 Z( 6, 6 ) = -E( JS, JS ) 00551 Z( 8, 6 ) = -E( JS, JSP1 ) 00552 * 00553 Z( 1, 7 ) = -B( JSP1, JS ) 00554 Z( 3, 7 ) = -B( JSP1, JSP1 ) 00555 Z( 7, 7 ) = -E( JSP1, JSP1 ) 00556 * 00557 Z( 2, 8 ) = -B( JSP1, JS ) 00558 Z( 4, 8 ) = -B( JSP1, JSP1 ) 00559 Z( 8, 8 ) = -E( JSP1, JSP1 ) 00560 * 00561 * Set up right hand side(s) 00562 * 00563 K = 1 00564 II = MB*NB + 1 00565 DO 80 JJ = 0, NB - 1 00566 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 00567 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 00568 K = K + MB 00569 II = II + MB 00570 80 CONTINUE 00571 * 00572 * Solve Z * x = RHS 00573 * 00574 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00575 IF( IERR.GT.0 ) 00576 $ INFO = IERR 00577 IF( IJOB.EQ.0 ) THEN 00578 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, 00579 $ SCALOC ) 00580 IF( SCALOC.NE.ONE ) THEN 00581 DO 90 K = 1, N 00582 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00583 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00584 90 CONTINUE 00585 SCALE = SCALE*SCALOC 00586 END IF 00587 ELSE 00588 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM, 00589 $ RDSCAL, IPIV, JPIV ) 00590 END IF 00591 * 00592 * Unpack solution vector(s) 00593 * 00594 K = 1 00595 II = MB*NB + 1 00596 DO 100 JJ = 0, NB - 1 00597 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 00598 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 00599 K = K + MB 00600 II = II + MB 00601 100 CONTINUE 00602 * 00603 * Substitute R(I, J) and L(I, J) into remaining 00604 * equation. 00605 * 00606 IF( I.GT.1 ) THEN 00607 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 00608 $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE, 00609 $ C( 1, JS ), LDC ) 00610 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE, 00611 $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE, 00612 $ F( 1, JS ), LDF ) 00613 END IF 00614 IF( J.LT.Q ) THEN 00615 K = MB*NB + 1 00616 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 00617 $ MB, B( JS, JE+1 ), LDB, ONE, 00618 $ C( IS, JE+1 ), LDC ) 00619 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ), 00620 $ MB, E( JS, JE+1 ), LDE, ONE, 00621 $ F( IS, JE+1 ), LDF ) 00622 END IF 00623 * 00624 END IF 00625 * 00626 110 CONTINUE 00627 120 CONTINUE 00628 ELSE 00629 * 00630 * Solve (I, J) - subsystem 00631 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J) 00632 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) 00633 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1 00634 * 00635 SCALE = ONE 00636 SCALOC = ONE 00637 DO 200 I = 1, P 00638 * 00639 IS = IWORK( I ) 00640 ISP1 = IS + 1 00641 IE = IWORK ( I+1 ) - 1 00642 MB = IE - IS + 1 00643 DO 190 J = Q, P + 2, -1 00644 * 00645 JS = IWORK( J ) 00646 JSP1 = JS + 1 00647 JE = IWORK( J+1 ) - 1 00648 NB = JE - JS + 1 00649 ZDIM = MB*NB*2 00650 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN 00651 * 00652 * Build a 2-by-2 system Z**T * x = RHS 00653 * 00654 Z( 1, 1 ) = A( IS, IS ) 00655 Z( 2, 1 ) = -B( JS, JS ) 00656 Z( 1, 2 ) = D( IS, IS ) 00657 Z( 2, 2 ) = -E( JS, JS ) 00658 * 00659 * Set up right hand side(s) 00660 * 00661 RHS( 1 ) = C( IS, JS ) 00662 RHS( 2 ) = F( IS, JS ) 00663 * 00664 * Solve Z**T * x = RHS 00665 * 00666 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00667 IF( IERR.GT.0 ) 00668 $ INFO = IERR 00669 * 00670 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00671 IF( SCALOC.NE.ONE ) THEN 00672 DO 130 K = 1, N 00673 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00674 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00675 130 CONTINUE 00676 SCALE = SCALE*SCALOC 00677 END IF 00678 * 00679 * Unpack solution vector(s) 00680 * 00681 C( IS, JS ) = RHS( 1 ) 00682 F( IS, JS ) = RHS( 2 ) 00683 * 00684 * Substitute R(I, J) and L(I, J) into remaining 00685 * equation. 00686 * 00687 IF( J.GT.P+2 ) THEN 00688 ALPHA = RHS( 1 ) 00689 CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ), 00690 $ LDF ) 00691 ALPHA = RHS( 2 ) 00692 CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ), 00693 $ LDF ) 00694 END IF 00695 IF( I.LT.P ) THEN 00696 ALPHA = -RHS( 1 ) 00697 CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA, 00698 $ C( IE+1, JS ), 1 ) 00699 ALPHA = -RHS( 2 ) 00700 CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD, 00701 $ C( IE+1, JS ), 1 ) 00702 END IF 00703 * 00704 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN 00705 * 00706 * Build a 4-by-4 system Z**T * x = RHS 00707 * 00708 Z( 1, 1 ) = A( IS, IS ) 00709 Z( 2, 1 ) = ZERO 00710 Z( 3, 1 ) = -B( JS, JS ) 00711 Z( 4, 1 ) = -B( JSP1, JS ) 00712 * 00713 Z( 1, 2 ) = ZERO 00714 Z( 2, 2 ) = A( IS, IS ) 00715 Z( 3, 2 ) = -B( JS, JSP1 ) 00716 Z( 4, 2 ) = -B( JSP1, JSP1 ) 00717 * 00718 Z( 1, 3 ) = D( IS, IS ) 00719 Z( 2, 3 ) = ZERO 00720 Z( 3, 3 ) = -E( JS, JS ) 00721 Z( 4, 3 ) = ZERO 00722 * 00723 Z( 1, 4 ) = ZERO 00724 Z( 2, 4 ) = D( IS, IS ) 00725 Z( 3, 4 ) = -E( JS, JSP1 ) 00726 Z( 4, 4 ) = -E( JSP1, JSP1 ) 00727 * 00728 * Set up right hand side(s) 00729 * 00730 RHS( 1 ) = C( IS, JS ) 00731 RHS( 2 ) = C( IS, JSP1 ) 00732 RHS( 3 ) = F( IS, JS ) 00733 RHS( 4 ) = F( IS, JSP1 ) 00734 * 00735 * Solve Z**T * x = RHS 00736 * 00737 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00738 IF( IERR.GT.0 ) 00739 $ INFO = IERR 00740 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00741 IF( SCALOC.NE.ONE ) THEN 00742 DO 140 K = 1, N 00743 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00744 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00745 140 CONTINUE 00746 SCALE = SCALE*SCALOC 00747 END IF 00748 * 00749 * Unpack solution vector(s) 00750 * 00751 C( IS, JS ) = RHS( 1 ) 00752 C( IS, JSP1 ) = RHS( 2 ) 00753 F( IS, JS ) = RHS( 3 ) 00754 F( IS, JSP1 ) = RHS( 4 ) 00755 * 00756 * Substitute R(I, J) and L(I, J) into remaining 00757 * equation. 00758 * 00759 IF( J.GT.P+2 ) THEN 00760 CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1, 00761 $ F( IS, 1 ), LDF ) 00762 CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1, 00763 $ F( IS, 1 ), LDF ) 00764 CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1, 00765 $ F( IS, 1 ), LDF ) 00766 CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1, 00767 $ F( IS, 1 ), LDF ) 00768 END IF 00769 IF( I.LT.P ) THEN 00770 CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA, 00771 $ RHS( 1 ), 1, C( IE+1, JS ), LDC ) 00772 CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD, 00773 $ RHS( 3 ), 1, C( IE+1, JS ), LDC ) 00774 END IF 00775 * 00776 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN 00777 * 00778 * Build a 4-by-4 system Z**T * x = RHS 00779 * 00780 Z( 1, 1 ) = A( IS, IS ) 00781 Z( 2, 1 ) = A( IS, ISP1 ) 00782 Z( 3, 1 ) = -B( JS, JS ) 00783 Z( 4, 1 ) = ZERO 00784 * 00785 Z( 1, 2 ) = A( ISP1, IS ) 00786 Z( 2, 2 ) = A( ISP1, ISP1 ) 00787 Z( 3, 2 ) = ZERO 00788 Z( 4, 2 ) = -B( JS, JS ) 00789 * 00790 Z( 1, 3 ) = D( IS, IS ) 00791 Z( 2, 3 ) = D( IS, ISP1 ) 00792 Z( 3, 3 ) = -E( JS, JS ) 00793 Z( 4, 3 ) = ZERO 00794 * 00795 Z( 1, 4 ) = ZERO 00796 Z( 2, 4 ) = D( ISP1, ISP1 ) 00797 Z( 3, 4 ) = ZERO 00798 Z( 4, 4 ) = -E( JS, JS ) 00799 * 00800 * Set up right hand side(s) 00801 * 00802 RHS( 1 ) = C( IS, JS ) 00803 RHS( 2 ) = C( ISP1, JS ) 00804 RHS( 3 ) = F( IS, JS ) 00805 RHS( 4 ) = F( ISP1, JS ) 00806 * 00807 * Solve Z**T * x = RHS 00808 * 00809 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00810 IF( IERR.GT.0 ) 00811 $ INFO = IERR 00812 * 00813 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00814 IF( SCALOC.NE.ONE ) THEN 00815 DO 150 K = 1, N 00816 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00817 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00818 150 CONTINUE 00819 SCALE = SCALE*SCALOC 00820 END IF 00821 * 00822 * Unpack solution vector(s) 00823 * 00824 C( IS, JS ) = RHS( 1 ) 00825 C( ISP1, JS ) = RHS( 2 ) 00826 F( IS, JS ) = RHS( 3 ) 00827 F( ISP1, JS ) = RHS( 4 ) 00828 * 00829 * Substitute R(I, J) and L(I, J) into remaining 00830 * equation. 00831 * 00832 IF( J.GT.P+2 ) THEN 00833 CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ), 00834 $ 1, F( IS, 1 ), LDF ) 00835 CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ), 00836 $ 1, F( IS, 1 ), LDF ) 00837 END IF 00838 IF( I.LT.P ) THEN 00839 CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ), 00840 $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ), 00841 $ 1 ) 00842 CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ), 00843 $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ), 00844 $ 1 ) 00845 END IF 00846 * 00847 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN 00848 * 00849 * Build an 8-by-8 system Z**T * x = RHS 00850 * 00851 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ ) 00852 * 00853 Z( 1, 1 ) = A( IS, IS ) 00854 Z( 2, 1 ) = A( IS, ISP1 ) 00855 Z( 5, 1 ) = -B( JS, JS ) 00856 Z( 7, 1 ) = -B( JSP1, JS ) 00857 * 00858 Z( 1, 2 ) = A( ISP1, IS ) 00859 Z( 2, 2 ) = A( ISP1, ISP1 ) 00860 Z( 6, 2 ) = -B( JS, JS ) 00861 Z( 8, 2 ) = -B( JSP1, JS ) 00862 * 00863 Z( 3, 3 ) = A( IS, IS ) 00864 Z( 4, 3 ) = A( IS, ISP1 ) 00865 Z( 5, 3 ) = -B( JS, JSP1 ) 00866 Z( 7, 3 ) = -B( JSP1, JSP1 ) 00867 * 00868 Z( 3, 4 ) = A( ISP1, IS ) 00869 Z( 4, 4 ) = A( ISP1, ISP1 ) 00870 Z( 6, 4 ) = -B( JS, JSP1 ) 00871 Z( 8, 4 ) = -B( JSP1, JSP1 ) 00872 * 00873 Z( 1, 5 ) = D( IS, IS ) 00874 Z( 2, 5 ) = D( IS, ISP1 ) 00875 Z( 5, 5 ) = -E( JS, JS ) 00876 * 00877 Z( 2, 6 ) = D( ISP1, ISP1 ) 00878 Z( 6, 6 ) = -E( JS, JS ) 00879 * 00880 Z( 3, 7 ) = D( IS, IS ) 00881 Z( 4, 7 ) = D( IS, ISP1 ) 00882 Z( 5, 7 ) = -E( JS, JSP1 ) 00883 Z( 7, 7 ) = -E( JSP1, JSP1 ) 00884 * 00885 Z( 4, 8 ) = D( ISP1, ISP1 ) 00886 Z( 6, 8 ) = -E( JS, JSP1 ) 00887 Z( 8, 8 ) = -E( JSP1, JSP1 ) 00888 * 00889 * Set up right hand side(s) 00890 * 00891 K = 1 00892 II = MB*NB + 1 00893 DO 160 JJ = 0, NB - 1 00894 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 ) 00895 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 ) 00896 K = K + MB 00897 II = II + MB 00898 160 CONTINUE 00899 * 00900 * 00901 * Solve Z**T * x = RHS 00902 * 00903 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR ) 00904 IF( IERR.GT.0 ) 00905 $ INFO = IERR 00906 * 00907 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC ) 00908 IF( SCALOC.NE.ONE ) THEN 00909 DO 170 K = 1, N 00910 CALL DSCAL( M, SCALOC, C( 1, K ), 1 ) 00911 CALL DSCAL( M, SCALOC, F( 1, K ), 1 ) 00912 170 CONTINUE 00913 SCALE = SCALE*SCALOC 00914 END IF 00915 * 00916 * Unpack solution vector(s) 00917 * 00918 K = 1 00919 II = MB*NB + 1 00920 DO 180 JJ = 0, NB - 1 00921 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 ) 00922 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 ) 00923 K = K + MB 00924 II = II + MB 00925 180 CONTINUE 00926 * 00927 * Substitute R(I, J) and L(I, J) into remaining 00928 * equation. 00929 * 00930 IF( J.GT.P+2 ) THEN 00931 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE, 00932 $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE, 00933 $ F( IS, 1 ), LDF ) 00934 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE, 00935 $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE, 00936 $ F( IS, 1 ), LDF ) 00937 END IF 00938 IF( I.LT.P ) THEN 00939 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 00940 $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC, 00941 $ ONE, C( IE+1, JS ), LDC ) 00942 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE, 00943 $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF, 00944 $ ONE, C( IE+1, JS ), LDC ) 00945 END IF 00946 * 00947 END IF 00948 * 00949 190 CONTINUE 00950 200 CONTINUE 00951 * 00952 END IF 00953 RETURN 00954 * 00955 * End of DTGSY2 00956 * 00957 END
1.7.3 
