Annotation of rpl/lapack/lapack/dtgexc.f, revision 1.4

1.1       bertrand    1:       SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                      2:      $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
                      3: *
                      4: *  -- LAPACK routine (version 3.2) --
                      5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      7: *     November 2006
                      8: *
                      9: *     .. Scalar Arguments ..
                     10:       LOGICAL            WANTQ, WANTZ
                     11:       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
                     12: *     ..
                     13: *     .. Array Arguments ..
                     14:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
                     15:      $                   WORK( * ), Z( LDZ, * )
                     16: *     ..
                     17: *
                     18: *  Purpose
                     19: *  =======
                     20: *
                     21: *  DTGEXC reorders the generalized real Schur decomposition of a real
                     22: *  matrix pair (A,B) using an orthogonal equivalence transformation
                     23: *
                     24: *                 (A, B) = Q * (A, B) * Z',
                     25: *
                     26: *  so that the diagonal block of (A, B) with row index IFST is moved
                     27: *  to row ILST.
                     28: *
                     29: *  (A, B) must be in generalized real Schur canonical form (as returned
                     30: *  by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
                     31: *  diagonal blocks. B is upper triangular.
                     32: *
                     33: *  Optionally, the matrices Q and Z of generalized Schur vectors are
                     34: *  updated.
                     35: *
                     36: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
                     37: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
                     38: *
                     39: *
                     40: *  Arguments
                     41: *  =========
                     42: *
                     43: *  WANTQ   (input) LOGICAL
                     44: *          .TRUE. : update the left transformation matrix Q;
                     45: *          .FALSE.: do not update Q.
                     46: *
                     47: *  WANTZ   (input) LOGICAL
                     48: *          .TRUE. : update the right transformation matrix Z;
                     49: *          .FALSE.: do not update Z.
                     50: *
                     51: *  N       (input) INTEGER
                     52: *          The order of the matrices A and B. N >= 0.
                     53: *
                     54: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
                     55: *          On entry, the matrix A in generalized real Schur canonical
                     56: *          form.
                     57: *          On exit, the updated matrix A, again in generalized
                     58: *          real Schur canonical form.
                     59: *
                     60: *  LDA     (input)  INTEGER
                     61: *          The leading dimension of the array A. LDA >= max(1,N).
                     62: *
                     63: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
                     64: *          On entry, the matrix B in generalized real Schur canonical
                     65: *          form (A,B).
                     66: *          On exit, the updated matrix B, again in generalized
                     67: *          real Schur canonical form (A,B).
                     68: *
                     69: *  LDB     (input)  INTEGER
                     70: *          The leading dimension of the array B. LDB >= max(1,N).
                     71: *
                     72: *  Q       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
                     73: *          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
                     74: *          On exit, the updated matrix Q.
                     75: *          If WANTQ = .FALSE., Q is not referenced.
                     76: *
                     77: *  LDQ     (input) INTEGER
                     78: *          The leading dimension of the array Q. LDQ >= 1.
                     79: *          If WANTQ = .TRUE., LDQ >= N.
                     80: *
                     81: *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
                     82: *          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
                     83: *          On exit, the updated matrix Z.
                     84: *          If WANTZ = .FALSE., Z is not referenced.
                     85: *
                     86: *  LDZ     (input) INTEGER
                     87: *          The leading dimension of the array Z. LDZ >= 1.
                     88: *          If WANTZ = .TRUE., LDZ >= N.
                     89: *
                     90: *  IFST    (input/output) INTEGER
                     91: *  ILST    (input/output) INTEGER
                     92: *          Specify the reordering of the diagonal blocks of (A, B).
                     93: *          The block with row index IFST is moved to row ILST, by a
                     94: *          sequence of swapping between adjacent blocks.
                     95: *          On exit, if IFST pointed on entry to the second row of
                     96: *          a 2-by-2 block, it is changed to point to the first row;
                     97: *          ILST always points to the first row of the block in its
                     98: *          final position (which may differ from its input value by
                     99: *          +1 or -1). 1 <= IFST, ILST <= N.
                    100: *
                    101: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                    102: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                    103: *
                    104: *  LWORK   (input) INTEGER
                    105: *          The dimension of the array WORK.
                    106: *          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
                    107: *
                    108: *          If LWORK = -1, then a workspace query is assumed; the routine
                    109: *          only calculates the optimal size of the WORK array, returns
                    110: *          this value as the first entry of the WORK array, and no error
                    111: *          message related to LWORK is issued by XERBLA.
                    112: *
                    113: *  INFO    (output) INTEGER
                    114: *           =0:  successful exit.
                    115: *           <0:  if INFO = -i, the i-th argument had an illegal value.
                    116: *           =1:  The transformed matrix pair (A, B) would be too far
                    117: *                from generalized Schur form; the problem is ill-
                    118: *                conditioned. (A, B) may have been partially reordered,
                    119: *                and ILST points to the first row of the current
                    120: *                position of the block being moved.
                    121: *
                    122: *  Further Details
                    123: *  ===============
                    124: *
                    125: *  Based on contributions by
                    126: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
                    127: *     Umea University, S-901 87 Umea, Sweden.
                    128: *
                    129: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
                    130: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
                    131: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
                    132: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
                    133: *
                    134: *  =====================================================================
                    135: *
                    136: *     .. Parameters ..
                    137:       DOUBLE PRECISION   ZERO
                    138:       PARAMETER          ( ZERO = 0.0D+0 )
                    139: *     ..
                    140: *     .. Local Scalars ..
                    141:       LOGICAL            LQUERY
                    142:       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
                    143: *     ..
                    144: *     .. External Subroutines ..
                    145:       EXTERNAL           DTGEX2, XERBLA
                    146: *     ..
                    147: *     .. Intrinsic Functions ..
                    148:       INTRINSIC          MAX
                    149: *     ..
                    150: *     .. Executable Statements ..
                    151: *
                    152: *     Decode and test input arguments.
                    153: *
                    154:       INFO = 0
                    155:       LQUERY = ( LWORK.EQ.-1 )
                    156:       IF( N.LT.0 ) THEN
                    157:          INFO = -3
                    158:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    159:          INFO = -5
                    160:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
                    161:          INFO = -7
                    162:       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
                    163:          INFO = -9
                    164:       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
                    165:          INFO = -11
                    166:       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
                    167:          INFO = -12
                    168:       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
                    169:          INFO = -13
                    170:       END IF
                    171: *
                    172:       IF( INFO.EQ.0 ) THEN
                    173:          IF( N.LE.1 ) THEN
                    174:             LWMIN = 1
                    175:          ELSE
                    176:             LWMIN = 4*N + 16
                    177:          END IF
                    178:          WORK(1) = LWMIN
                    179: *
                    180:          IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
                    181:             INFO = -15
                    182:          END IF
                    183:       END IF
                    184: *
                    185:       IF( INFO.NE.0 ) THEN
                    186:          CALL XERBLA( 'DTGEXC', -INFO )
                    187:          RETURN
                    188:       ELSE IF( LQUERY ) THEN
                    189:          RETURN
                    190:       END IF
                    191: *
                    192: *     Quick return if possible
                    193: *
                    194:       IF( N.LE.1 )
                    195:      $   RETURN
                    196: *
                    197: *     Determine the first row of the specified block and find out
                    198: *     if it is 1-by-1 or 2-by-2.
                    199: *
                    200:       IF( IFST.GT.1 ) THEN
                    201:          IF( A( IFST, IFST-1 ).NE.ZERO )
                    202:      $      IFST = IFST - 1
                    203:       END IF
                    204:       NBF = 1
                    205:       IF( IFST.LT.N ) THEN
                    206:          IF( A( IFST+1, IFST ).NE.ZERO )
                    207:      $      NBF = 2
                    208:       END IF
                    209: *
                    210: *     Determine the first row of the final block
                    211: *     and find out if it is 1-by-1 or 2-by-2.
                    212: *
                    213:       IF( ILST.GT.1 ) THEN
                    214:          IF( A( ILST, ILST-1 ).NE.ZERO )
                    215:      $      ILST = ILST - 1
                    216:       END IF
                    217:       NBL = 1
                    218:       IF( ILST.LT.N ) THEN
                    219:          IF( A( ILST+1, ILST ).NE.ZERO )
                    220:      $      NBL = 2
                    221:       END IF
                    222:       IF( IFST.EQ.ILST )
                    223:      $   RETURN
                    224: *
                    225:       IF( IFST.LT.ILST ) THEN
                    226: *
                    227: *        Update ILST.
                    228: *
                    229:          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
                    230:      $      ILST = ILST - 1
                    231:          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
                    232:      $      ILST = ILST + 1
                    233: *
                    234:          HERE = IFST
                    235: *
                    236:    10    CONTINUE
                    237: *
                    238: *        Swap with next one below.
                    239: *
                    240:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
                    241: *
                    242: *           Current block either 1-by-1 or 2-by-2.
                    243: *
                    244:             NBNEXT = 1
                    245:             IF( HERE+NBF+1.LE.N ) THEN
                    246:                IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
                    247:      $            NBNEXT = 2
                    248:             END IF
                    249:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    250:      $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
                    251:             IF( INFO.NE.0 ) THEN
                    252:                ILST = HERE
                    253:                RETURN
                    254:             END IF
                    255:             HERE = HERE + NBNEXT
                    256: *
                    257: *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
                    258: *
                    259:             IF( NBF.EQ.2 ) THEN
                    260:                IF( A( HERE+1, HERE ).EQ.ZERO )
                    261:      $            NBF = 3
                    262:             END IF
                    263: *
                    264:          ELSE
                    265: *
                    266: *           Current block consists of two 1-by-1 blocks, each of which
                    267: *           must be swapped individually.
                    268: *
                    269:             NBNEXT = 1
                    270:             IF( HERE+3.LE.N ) THEN
                    271:                IF( A( HERE+3, HERE+2 ).NE.ZERO )
                    272:      $            NBNEXT = 2
                    273:             END IF
                    274:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    275:      $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
                    276:             IF( INFO.NE.0 ) THEN
                    277:                ILST = HERE
                    278:                RETURN
                    279:             END IF
                    280:             IF( NBNEXT.EQ.1 ) THEN
                    281: *
                    282: *              Swap two 1-by-1 blocks.
                    283: *
                    284:                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    285:      $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                    286:                IF( INFO.NE.0 ) THEN
                    287:                   ILST = HERE
                    288:                   RETURN
                    289:                END IF
                    290:                HERE = HERE + 1
                    291: *
                    292:             ELSE
                    293: *
                    294: *              Recompute NBNEXT in case of 2-by-2 split.
                    295: *
                    296:                IF( A( HERE+2, HERE+1 ).EQ.ZERO )
                    297:      $            NBNEXT = 1
                    298:                IF( NBNEXT.EQ.2 ) THEN
                    299: *
                    300: *                 2-by-2 block did not split.
                    301: *
                    302:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    303:      $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
                    304:      $                         INFO )
                    305:                   IF( INFO.NE.0 ) THEN
                    306:                      ILST = HERE
                    307:                      RETURN
                    308:                   END IF
                    309:                   HERE = HERE + 2
                    310:                ELSE
                    311: *
                    312: *                 2-by-2 block did split.
                    313: *
                    314:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    315:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                    316:                   IF( INFO.NE.0 ) THEN
                    317:                      ILST = HERE
                    318:                      RETURN
                    319:                   END IF
                    320:                   HERE = HERE + 1
                    321:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    322:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                    323:                   IF( INFO.NE.0 ) THEN
                    324:                      ILST = HERE
                    325:                      RETURN
                    326:                   END IF
                    327:                   HERE = HERE + 1
                    328:                END IF
                    329: *
                    330:             END IF
                    331:          END IF
                    332:          IF( HERE.LT.ILST )
                    333:      $      GO TO 10
                    334:       ELSE
                    335:          HERE = IFST
                    336: *
                    337:    20    CONTINUE
                    338: *
                    339: *        Swap with next one below.
                    340: *
                    341:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
                    342: *
                    343: *           Current block either 1-by-1 or 2-by-2.
                    344: *
                    345:             NBNEXT = 1
                    346:             IF( HERE.GE.3 ) THEN
                    347:                IF( A( HERE-1, HERE-2 ).NE.ZERO )
                    348:      $            NBNEXT = 2
                    349:             END IF
                    350:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    351:      $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
                    352:      $                   INFO )
                    353:             IF( INFO.NE.0 ) THEN
                    354:                ILST = HERE
                    355:                RETURN
                    356:             END IF
                    357:             HERE = HERE - NBNEXT
                    358: *
                    359: *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
                    360: *
                    361:             IF( NBF.EQ.2 ) THEN
                    362:                IF( A( HERE+1, HERE ).EQ.ZERO )
                    363:      $            NBF = 3
                    364:             END IF
                    365: *
                    366:          ELSE
                    367: *
                    368: *           Current block consists of two 1-by-1 blocks, each of which
                    369: *           must be swapped individually.
                    370: *
                    371:             NBNEXT = 1
                    372:             IF( HERE.GE.3 ) THEN
                    373:                IF( A( HERE-1, HERE-2 ).NE.ZERO )
                    374:      $            NBNEXT = 2
                    375:             END IF
                    376:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    377:      $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
                    378:      $                   INFO )
                    379:             IF( INFO.NE.0 ) THEN
                    380:                ILST = HERE
                    381:                RETURN
                    382:             END IF
                    383:             IF( NBNEXT.EQ.1 ) THEN
                    384: *
                    385: *              Swap two 1-by-1 blocks.
                    386: *
                    387:                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
                    388:      $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
                    389:                IF( INFO.NE.0 ) THEN
                    390:                   ILST = HERE
                    391:                   RETURN
                    392:                END IF
                    393:                HERE = HERE - 1
                    394:             ELSE
                    395: *
                    396: *             Recompute NBNEXT in case of 2-by-2 split.
                    397: *
                    398:                IF( A( HERE, HERE-1 ).EQ.ZERO )
                    399:      $            NBNEXT = 1
                    400:                IF( NBNEXT.EQ.2 ) THEN
                    401: *
                    402: *                 2-by-2 block did not split.
                    403: *
                    404:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    405:      $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
                    406:                   IF( INFO.NE.0 ) THEN
                    407:                      ILST = HERE
                    408:                      RETURN
                    409:                   END IF
                    410:                   HERE = HERE - 2
                    411:                ELSE
                    412: *
                    413: *                 2-by-2 block did split.
                    414: *
                    415:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    416:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                    417:                   IF( INFO.NE.0 ) THEN
                    418:                      ILST = HERE
                    419:                      RETURN
                    420:                   END IF
                    421:                   HERE = HERE - 1
                    422:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
                    423:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
                    424:                   IF( INFO.NE.0 ) THEN
                    425:                      ILST = HERE
                    426:                      RETURN
                    427:                   END IF
                    428:                   HERE = HERE - 1
                    429:                END IF
                    430:             END IF
                    431:          END IF
                    432:          IF( HERE.GT.ILST )
                    433:      $      GO TO 20
                    434:       END IF
                    435:       ILST = HERE
                    436:       WORK( 1 ) = LWMIN
                    437:       RETURN
                    438: *
                    439: *     End of DTGEXC
                    440: *
                    441:       END

CVSweb interface <joel.bertrand@systella.fr>