File:  [local] / rpl / lapack / lapack / dtgsyl.f
Revision 1.1.1.1 (vendor branch): download - view: text, annotated - select for diffs - revision graph
Tue Jan 26 15:22:45 2010 UTC (14 years, 3 months ago) by bertrand
Branches: JKB
CVS tags: start, rpl-4_0_14, rpl-4_0_13, rpl-4_0_12, rpl-4_0_11, rpl-4_0_10


Commit initial.

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

CVSweb interface <joel.bertrand@systella.fr>