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