SRC/dtgsyl.f

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