SRC/ctgsyl.f

   1 *> \brief \b CTGSYL
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CTGSYL + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsyl.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsyl.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsyl.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
  22 *                          LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
  23 *                          IWORK, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          TRANS
  27 *       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
  28 *      $                   LWORK, M, N
  29 *       REAL               DIF, SCALE
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            IWORK( * )
  33 *       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
  34 *      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
  35 *      $                   WORK( * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> CTGSYL solves the generalized Sylvester equation:
  45 *>
  46 *>             A * R - L * B = scale * C            (1)
  47 *>             D * R - L * E = scale * F
  48 *>
  49 *> where R and L are unknown m-by-n matrices, (A, D), (B, E) and
  50 *> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
  51 *> respectively, with complex entries. A, B, D and E are upper
  52 *> triangular (i.e., (A,D) and (B,E) in generalized Schur form).
  53 *>
  54 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1
  55 *> is an output scaling factor chosen to avoid overflow.
  56 *>
  57 *> In matrix notation (1) is equivalent to solve Zx = scale*b, where Z
  58 *> is defined as
  59 *>
  60 *>        Z = [ kron(In, A)  -kron(B**H, Im) ]        (2)
  61 *>            [ kron(In, D)  -kron(E**H, Im) ],
  62 *>
  63 *> Here Ix is the identity matrix of size x and X**H is the conjugate
  64 *> transpose of X. Kron(X, Y) is the Kronecker product between the
  65 *> matrices X and Y.
  66 *>
  67 *> If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b
  68 *> is solved for, which is equivalent to solve for R and L in
  69 *>
  70 *>             A**H * R + D**H * L = scale * C           (3)
  71 *>             R * B**H + L * E**H = scale * -F
  72 *>
  73 *> This case (TRANS = 'C') is used to compute an one-norm-based estimate
  74 *> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
  75 *> and (B,E), using CLACON.
  76 *>
  77 *> If IJOB >= 1, CTGSYL computes a Frobenius norm-based estimate of
  78 *> Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
  79 *> reciprocal of the smallest singular value of Z.
  80 *>
  81 *> This is a level-3 BLAS algorithm.
  82 *> \endverbatim
  83 *
  84 *  Arguments:
  85 *  ==========
  86 *
  87 *> \param[in] TRANS
  88 *> \verbatim
  89 *>          TRANS is CHARACTER*1
  90 *>          = 'N': solve the generalized sylvester equation (1).
  91 *>          = 'C': solve the "conjugate transposed" system (3).
  92 *> \endverbatim
  93 *>
  94 *> \param[in] IJOB
  95 *> \verbatim
  96 *>          IJOB is INTEGER
  97 *>          Specifies what kind of functionality to be performed.
  98 *>          =0: solve (1) only.
  99 *>          =1: The functionality of 0 and 3.
 100 *>          =2: The functionality of 0 and 4.
 101 *>          =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
 102 *>              (look ahead strategy is used).
 103 *>          =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
 104 *>              (CGECON on sub-systems is used).
 105 *>          Not referenced if TRANS = 'C'.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] M
 109 *> \verbatim
 110 *>          M is INTEGER
 111 *>          The order of the matrices A and D, and the row dimension of
 112 *>          the matrices C, F, R and L.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] N
 116 *> \verbatim
 117 *>          N is INTEGER
 118 *>          The order of the matrices B and E, and the column dimension
 119 *>          of the matrices C, F, R and L.
 120 *> \endverbatim
 121 *>
 122 *> \param[in] A
 123 *> \verbatim
 124 *>          A is COMPLEX array, dimension (LDA, M)
 125 *>          The upper triangular matrix A.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LDA
 129 *> \verbatim
 130 *>          LDA is INTEGER
 131 *>          The leading dimension of the array A. LDA >= max(1, M).
 132 *> \endverbatim
 133 *>
 134 *> \param[in] B
 135 *> \verbatim
 136 *>          B is COMPLEX array, dimension (LDB, N)
 137 *>          The upper triangular matrix B.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] LDB
 141 *> \verbatim
 142 *>          LDB is INTEGER
 143 *>          The leading dimension of the array B. LDB >= max(1, N).
 144 *> \endverbatim
 145 *>
 146 *> \param[in,out] C
 147 *> \verbatim
 148 *>          C is COMPLEX array, dimension (LDC, N)
 149 *>          On entry, C contains the right-hand-side of the first matrix
 150 *>          equation in (1) or (3).
 151 *>          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
 152 *>          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
 153 *>          the solution achieved during the computation of the
 154 *>          Dif-estimate.
 155 *> \endverbatim
 156 *>
 157 *> \param[in] LDC
 158 *> \verbatim
 159 *>          LDC is INTEGER
 160 *>          The leading dimension of the array C. LDC >= max(1, M).
 161 *> \endverbatim
 162 *>
 163 *> \param[in] D
 164 *> \verbatim
 165 *>          D is COMPLEX array, dimension (LDD, M)
 166 *>          The upper triangular matrix D.
 167 *> \endverbatim
 168 *>
 169 *> \param[in] LDD
 170 *> \verbatim
 171 *>          LDD is INTEGER
 172 *>          The leading dimension of the array D. LDD >= max(1, M).
 173 *> \endverbatim
 174 *>
 175 *> \param[in] E
 176 *> \verbatim
 177 *>          E is COMPLEX array, dimension (LDE, N)
 178 *>          The upper triangular matrix E.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] LDE
 182 *> \verbatim
 183 *>          LDE is INTEGER
 184 *>          The leading dimension of the array E. LDE >= max(1, N).
 185 *> \endverbatim
 186 *>
 187 *> \param[in,out] F
 188 *> \verbatim
 189 *>          F is COMPLEX array, dimension (LDF, N)
 190 *>          On entry, F contains the right-hand-side of the second matrix
 191 *>          equation in (1) or (3).
 192 *>          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
 193 *>          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
 194 *>          the solution achieved during the computation of the
 195 *>          Dif-estimate.
 196 *> \endverbatim
 197 *>
 198 *> \param[in] LDF
 199 *> \verbatim
 200 *>          LDF is INTEGER
 201 *>          The leading dimension of the array F. LDF >= max(1, M).
 202 *> \endverbatim
 203 *>
 204 *> \param[out] DIF
 205 *> \verbatim
 206 *>          DIF is REAL
 207 *>          On exit DIF is the reciprocal of a lower bound of the
 208 *>          reciprocal of the Dif-function, i.e. DIF is an upper bound of
 209 *>          Dif[(A,D), (B,E)] = sigma-min(Z), where Z as in (2).
 210 *>          IF IJOB = 0 or TRANS = 'C', DIF is not referenced.
 211 *> \endverbatim
 212 *>
 213 *> \param[out] SCALE
 214 *> \verbatim
 215 *>          SCALE is REAL
 216 *>          On exit SCALE is the scaling factor in (1) or (3).
 217 *>          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
 218 *>          to a slightly perturbed system but the input matrices A, B,
 219 *>          D and E have not been changed. If SCALE = 0, R and L will
 220 *>          hold the solutions to the homogenious system with C = F = 0.
 221 *> \endverbatim
 222 *>
 223 *> \param[out] WORK
 224 *> \verbatim
 225 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 226 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 227 *> \endverbatim
 228 *>
 229 *> \param[in] LWORK
 230 *> \verbatim
 231 *>          LWORK is INTEGER
 232 *>          The dimension of the array WORK. LWORK > = 1.
 233 *>          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
 234 *>
 235 *>          If LWORK = -1, then a workspace query is assumed; the routine
 236 *>          only calculates the optimal size of the WORK array, returns
 237 *>          this value as the first entry of the WORK array, and no error
 238 *>          message related to LWORK is issued by XERBLA.
 239 *> \endverbatim
 240 *>
 241 *> \param[out] IWORK
 242 *> \verbatim
 243 *>          IWORK is INTEGER array, dimension (M+N+2)
 244 *> \endverbatim
 245 *>
 246 *> \param[out] INFO
 247 *> \verbatim
 248 *>          INFO is INTEGER
 249 *>            =0: successful exit
 250 *>            <0: If INFO = -i, the i-th argument had an illegal value.
 251 *>            >0: (A, D) and (B, E) have common or very close
 252 *>                eigenvalues.
 253 *> \endverbatim
 254 *
 255 *  Authors:
 256 *  ========
 257 *
 258 *> \author Univ. of Tennessee
 259 *> \author Univ. of California Berkeley
 260 *> \author Univ. of Colorado Denver
 261 *> \author NAG Ltd.
 262 *
 263 *> \date November 2011
 264 *
 265 *> \ingroup complexSYcomputational
 266 *
 267 *> \par Contributors:
 268 *  ==================
 269 *>
 270 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 271 *>     Umea University, S-901 87 Umea, Sweden.
 272 *
 273 *> \par References:
 274 *  ================
 275 *>
 276 *>  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
 277 *>      for Solving the Generalized Sylvester Equation and Estimating the
 278 *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
 279 *>      Department of Computing Science, Umea University, S-901 87 Umea,
 280 *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
 281 *>      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
 282 *>      No 1, 1996.
 283 *> \n
 284 *>  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
 285 *>      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
 286 *>      Appl., 15(4):1045-1060, 1994.
 287 *> \n
 288 *>  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
 289 *>      Condition Estimators for Solving the Generalized Sylvester
 290 *>      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
 291 *>      July 1989, pp 745-751.
 292 *>
 293 *  =====================================================================
 294       SUBROUTINE CTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
 295      $                   LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
 296      $                   IWORK, INFO )
 297 *
 298 *  -- LAPACK computational routine (version 3.4.0) --
 299 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 300 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 301 *     November 2011
 302 *
 303 *     .. Scalar Arguments ..
 304       CHARACTER          TRANS
 305       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
 306      $                   LWORK, M, N
 307       REAL               DIF, SCALE
 308 *     ..
 309 *     .. Array Arguments ..
 310       INTEGER            IWORK( * )
 311       COMPLEX            A( LDA, * ), B( LDB, * ), C( LDC, * ),
 312      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
 313      $                   WORK( * )
 314 *     ..
 315 *
 316 *  =====================================================================
 317 *  Replaced various illegal calls to CCOPY by calls to CLASET.
 318 *  Sven Hammarling, 1/5/02.
 319 *
 320 *     .. Parameters ..
 321       REAL               ZERO, ONE
 322       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 323       COMPLEX            CZERO
 324       PARAMETER          ( CZERO = (0.0E+0, 0.0E+0) )
 325 *     ..
 326 *     .. Local Scalars ..
 327       LOGICAL            LQUERY, NOTRAN
 328       INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
 329      $                   LINFO, LWMIN, MB, NB, P, PQ, Q
 330       REAL               DSCALE, DSUM, SCALE2, SCALOC
 331 *     ..
 332 *     .. External Functions ..
 333       LOGICAL            LSAME
 334       INTEGER            ILAENV
 335       EXTERNAL           LSAME, ILAENV
 336 *     ..
 337 *     .. External Subroutines ..
 338       EXTERNAL           CGEMM, CLACPY, CLASET, CSCAL, CTGSY2, XERBLA
 339 *     ..
 340 *     .. Intrinsic Functions ..
 341       INTRINSIC          CMPLX, MAX, REAL, SQRT
 342 *     ..
 343 *     .. Executable Statements ..
 344 *
 345 *     Decode and test input parameters
 346 *
 347       INFO = 0
 348       NOTRAN = LSAME( TRANS, 'N' )
 349       LQUERY = ( LWORK.EQ.-1 )
 350 *
 351       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
 352          INFO = -1
 353       ELSE IF( NOTRAN ) THEN
 354          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
 355             INFO = -2
 356          END IF
 357       END IF
 358       IF( INFO.EQ.0 ) THEN
 359          IF( M.LE.0 ) THEN
 360             INFO = -3
 361          ELSE IF( N.LE.0 ) THEN
 362             INFO = -4
 363          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 364             INFO = -6
 365          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 366             INFO = -8
 367          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
 368             INFO = -10
 369          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
 370             INFO = -12
 371          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
 372             INFO = -14
 373          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
 374             INFO = -16
 375          END IF
 376       END IF
 377 *
 378       IF( INFO.EQ.0 ) THEN
 379          IF( NOTRAN ) THEN
 380             IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
 381                LWMIN = MAX( 1, 2*M*N )
 382             ELSE
 383                LWMIN = 1
 384             END IF
 385          ELSE
 386             LWMIN = 1
 387          END IF
 388          WORK( 1 ) = LWMIN
 389 *
 390          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 391             INFO = -20
 392          END IF
 393       END IF
 394 *
 395       IF( INFO.NE.0 ) THEN
 396          CALL XERBLA( 'CTGSYL', -INFO )
 397          RETURN
 398       ELSE IF( LQUERY ) THEN
 399          RETURN
 400       END IF
 401 *
 402 *     Quick return if possible
 403 *
 404       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 405          SCALE = 1
 406          IF( NOTRAN ) THEN
 407             IF( IJOB.NE.0 ) THEN
 408                DIF = 0
 409             END IF
 410          END IF
 411          RETURN
 412       END IF
 413 *
 414 *     Determine  optimal block sizes MB and NB
 415 *
 416       MB = ILAENV( 2, 'CTGSYL', TRANS, M, N, -1, -1 )
 417       NB = ILAENV( 5, 'CTGSYL', TRANS, M, N, -1, -1 )
 418 *
 419       ISOLVE = 1
 420       IFUNC = 0
 421       IF( NOTRAN ) THEN
 422          IF( IJOB.GE.3 ) THEN
 423             IFUNC = IJOB - 2
 424             CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
 425             CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
 426          ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
 427             ISOLVE = 2
 428          END IF
 429       END IF
 430 *
 431       IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
 432      $     THEN
 433 *
 434 *        Use unblocked Level 2 solver
 435 *
 436          DO 30 IROUND = 1, ISOLVE
 437 *
 438             SCALE = ONE
 439             DSCALE = ZERO
 440             DSUM = ONE
 441             PQ = M*N
 442             CALL CTGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
 443      $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
 444      $                   INFO )
 445             IF( DSCALE.NE.ZERO ) THEN
 446                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
 447                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
 448                ELSE
 449                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
 450                END IF
 451             END IF
 452             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
 453                IF( NOTRAN ) THEN
 454                   IFUNC = IJOB
 455                END IF
 456                SCALE2 = SCALE
 457                CALL CLACPY( 'F', M, N, C, LDC, WORK, M )
 458                CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
 459                CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
 460                CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
 461             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
 462                CALL CLACPY( 'F', M, N, WORK, M, C, LDC )
 463                CALL CLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
 464                SCALE = SCALE2
 465             END IF
 466    30    CONTINUE
 467 *
 468          RETURN
 469 *
 470       END IF
 471 *
 472 *     Determine block structure of A
 473 *
 474       P = 0
 475       I = 1
 476    40 CONTINUE
 477       IF( I.GT.M )
 478      $   GO TO 50
 479       P = P + 1
 480       IWORK( P ) = I
 481       I = I + MB
 482       IF( I.GE.M )
 483      $   GO TO 50
 484       GO TO 40
 485    50 CONTINUE
 486       IWORK( P+1 ) = M + 1
 487       IF( IWORK( P ).EQ.IWORK( P+1 ) )
 488      $   P = P - 1
 489 *
 490 *     Determine block structure of B
 491 *
 492       Q = P + 1
 493       J = 1
 494    60 CONTINUE
 495       IF( J.GT.N )
 496      $   GO TO 70
 497 *
 498       Q = Q + 1
 499       IWORK( Q ) = J
 500       J = J + NB
 501       IF( J.GE.N )
 502      $   GO TO 70
 503       GO TO 60
 504 *
 505    70 CONTINUE
 506       IWORK( Q+1 ) = N + 1
 507       IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
 508      $   Q = Q - 1
 509 *
 510       IF( NOTRAN ) THEN
 511          DO 150 IROUND = 1, ISOLVE
 512 *
 513 *           Solve (I, J) - subsystem
 514 *               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
 515 *               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
 516 *           for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
 517 *
 518             PQ = 0
 519             SCALE = ONE
 520             DSCALE = ZERO
 521             DSUM = ONE
 522             DO 130 J = P + 2, Q
 523                JS = IWORK( J )
 524                JE = IWORK( J+1 ) - 1
 525                NB = JE - JS + 1
 526                DO 120 I = P, 1, -1
 527                   IS = IWORK( I )
 528                   IE = IWORK( I+1 ) - 1
 529                   MB = IE - IS + 1
 530                   CALL CTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
 531      $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
 532      $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
 533      $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
 534      $                         LINFO )
 535                   IF( LINFO.GT.0 )
 536      $               INFO = LINFO
 537                   PQ = PQ + MB*NB
 538                   IF( SCALOC.NE.ONE ) THEN
 539                      DO 80 K = 1, JS - 1
 540                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
 541      $                              1 )
 542                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
 543      $                              1 )
 544    80                CONTINUE
 545                      DO 90 K = JS, JE
 546                         CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ),
 547      $                              C( 1, K ), 1 )
 548                         CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ),
 549      $                              F( 1, K ), 1 )
 550    90                CONTINUE
 551                      DO 100 K = JS, JE
 552                         CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
 553      $                              C( IE+1, K ), 1 )
 554                         CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
 555      $                              F( IE+1, K ), 1 )
 556   100                CONTINUE
 557                      DO 110 K = JE + 1, N
 558                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
 559      $                              1 )
 560                         CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
 561      $                              1 )
 562   110                CONTINUE
 563                      SCALE = SCALE*SCALOC
 564                   END IF
 565 *
 566 *                 Substitute R(I,J) and L(I,J) into remaining equation.
 567 *
 568                   IF( I.GT.1 ) THEN
 569                      CALL CGEMM( 'N', 'N', IS-1, NB, MB,
 570      $                           CMPLX( -ONE, ZERO ), A( 1, IS ), LDA,
 571      $                           C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
 572      $                           C( 1, JS ), LDC )
 573                      CALL CGEMM( 'N', 'N', IS-1, NB, MB,
 574      $                           CMPLX( -ONE, ZERO ), D( 1, IS ), LDD,
 575      $                           C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
 576      $                           F( 1, JS ), LDF )
 577                   END IF
 578                   IF( J.LT.Q ) THEN
 579                      CALL CGEMM( 'N', 'N', MB, N-JE, NB,
 580      $                           CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
 581      $                           B( JS, JE+1 ), LDB, CMPLX( ONE, ZERO ),
 582      $                           C( IS, JE+1 ), LDC )
 583                      CALL CGEMM( 'N', 'N', MB, N-JE, NB,
 584      $                           CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
 585      $                           E( JS, JE+1 ), LDE, CMPLX( ONE, ZERO ),
 586      $                           F( IS, JE+1 ), LDF )
 587                   END IF
 588   120          CONTINUE
 589   130       CONTINUE
 590             IF( DSCALE.NE.ZERO ) THEN
 591                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
 592                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
 593                ELSE
 594                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
 595                END IF
 596             END IF
 597             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
 598                IF( NOTRAN ) THEN
 599                   IFUNC = IJOB
 600                END IF
 601                SCALE2 = SCALE
 602                CALL CLACPY( 'F', M, N, C, LDC, WORK, M )
 603                CALL CLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
 604                CALL CLASET( 'F', M, N, CZERO, CZERO, C, LDC )
 605                CALL CLASET( 'F', M, N, CZERO, CZERO, F, LDF )
 606             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
 607                CALL CLACPY( 'F', M, N, WORK, M, C, LDC )
 608                CALL CLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
 609                SCALE = SCALE2
 610             END IF
 611   150    CONTINUE
 612       ELSE
 613 *
 614 *        Solve transposed (I, J)-subsystem
 615 *            A(I, I)**H * R(I, J) + D(I, I)**H * L(I, J) = C(I, J)
 616 *            R(I, J) * B(J, J)  + L(I, J) * E(J, J) = -F(I, J)
 617 *        for I = 1,2,..., P; J = Q, Q-1,..., 1
 618 *
 619          SCALE = ONE
 620          DO 210 I = 1, P
 621             IS = IWORK( I )
 622             IE = IWORK( I+1 ) - 1
 623             MB = IE - IS + 1
 624             DO 200 J = Q, P + 2, -1
 625                JS = IWORK( J )
 626                JE = IWORK( J+1 ) - 1
 627                NB = JE - JS + 1
 628                CALL CTGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
 629      $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
 630      $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
 631      $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
 632      $                      LINFO )
 633                IF( LINFO.GT.0 )
 634      $            INFO = LINFO
 635                IF( SCALOC.NE.ONE ) THEN
 636                   DO 160 K = 1, JS - 1
 637                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
 638      $                           1 )
 639                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
 640      $                           1 )
 641   160             CONTINUE
 642                   DO 170 K = JS, JE
 643                      CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), C( 1, K ),
 644      $                           1 )
 645                      CALL CSCAL( IS-1, CMPLX( SCALOC, ZERO ), F( 1, K ),
 646      $                           1 )
 647   170             CONTINUE
 648                   DO 180 K = JS, JE
 649                      CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
 650      $                           C( IE+1, K ), 1 )
 651                      CALL CSCAL( M-IE, CMPLX( SCALOC, ZERO ),
 652      $                           F( IE+1, K ), 1 )
 653   180             CONTINUE
 654                   DO 190 K = JE + 1, N
 655                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), C( 1, K ),
 656      $                           1 )
 657                      CALL CSCAL( M, CMPLX( SCALOC, ZERO ), F( 1, K ),
 658      $                           1 )
 659   190             CONTINUE
 660                   SCALE = SCALE*SCALOC
 661                END IF
 662 *
 663 *              Substitute R(I,J) and L(I,J) into remaining equation.
 664 *
 665                IF( J.GT.P+2 ) THEN
 666                   CALL CGEMM( 'N', 'C', MB, JS-1, NB,
 667      $                        CMPLX( ONE, ZERO ), C( IS, JS ), LDC,
 668      $                        B( 1, JS ), LDB, CMPLX( ONE, ZERO ),
 669      $                        F( IS, 1 ), LDF )
 670                   CALL CGEMM( 'N', 'C', MB, JS-1, NB,
 671      $                        CMPLX( ONE, ZERO ), F( IS, JS ), LDF,
 672      $                        E( 1, JS ), LDE, CMPLX( ONE, ZERO ),
 673      $                        F( IS, 1 ), LDF )
 674                END IF
 675                IF( I.LT.P ) THEN
 676                   CALL CGEMM( 'C', 'N', M-IE, NB, MB,
 677      $                        CMPLX( -ONE, ZERO ), A( IS, IE+1 ), LDA,
 678      $                        C( IS, JS ), LDC, CMPLX( ONE, ZERO ),
 679      $                        C( IE+1, JS ), LDC )
 680                   CALL CGEMM( 'C', 'N', M-IE, NB, MB,
 681      $                        CMPLX( -ONE, ZERO ), D( IS, IE+1 ), LDD,
 682      $                        F( IS, JS ), LDF, CMPLX( ONE, ZERO ),
 683      $                        C( IE+1, JS ), LDC )
 684                END IF
 685   200       CONTINUE
 686   210    CONTINUE
 687       END IF
 688 *
 689       WORK( 1 ) = LWMIN
 690 *
 691       RETURN
 692 *
 693 *     End of CTGSYL
 694 *
 695       END