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