Annotation of rpl/lapack/lapack/ztgsyl.f, revision 1.1.1.1

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

CVSweb interface <joel.bertrand@systella.fr>