Annotation of rpl/lapack/lapack/dtgsyl.f, revision 1.6

1.1       bertrand    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>