Annotation of rpl/lapack/lapack/dtrsna.f, revision 1.19

1.9       bertrand    1: *> \brief \b DTRSNA
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.15      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.9       bertrand    7: *
                      8: *> \htmlonly
1.15      bertrand    9: *> Download DTRSNA + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
1.9       bertrand   15: *> [TXT]</a>
1.15      bertrand   16: *> \endhtmlonly
1.9       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
                     22: *                          LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
                     23: *                          INFO )
1.15      bertrand   24: *
1.9       bertrand   25: *       .. Scalar Arguments ..
                     26: *       CHARACTER          HOWMNY, JOB
                     27: *       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
                     28: *       ..
                     29: *       .. Array Arguments ..
                     30: *       LOGICAL            SELECT( * )
                     31: *       INTEGER            IWORK( * )
                     32: *       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
                     33: *      $                   VR( LDVR, * ), WORK( LDWORK, * )
                     34: *       ..
1.15      bertrand   35: *
1.9       bertrand   36: *
                     37: *> \par Purpose:
                     38: *  =============
                     39: *>
                     40: *> \verbatim
                     41: *>
                     42: *> DTRSNA estimates reciprocal condition numbers for specified
                     43: *> eigenvalues and/or right eigenvectors of a real upper
                     44: *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
                     45: *> orthogonal).
                     46: *>
                     47: *> T must be in Schur canonical form (as returned by DHSEQR), that is,
                     48: *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
                     49: *> 2-by-2 diagonal block has its diagonal elements equal and its
                     50: *> off-diagonal elements of opposite sign.
                     51: *> \endverbatim
                     52: *
                     53: *  Arguments:
                     54: *  ==========
                     55: *
                     56: *> \param[in] JOB
                     57: *> \verbatim
                     58: *>          JOB is CHARACTER*1
                     59: *>          Specifies whether condition numbers are required for
                     60: *>          eigenvalues (S) or eigenvectors (SEP):
                     61: *>          = 'E': for eigenvalues only (S);
                     62: *>          = 'V': for eigenvectors only (SEP);
                     63: *>          = 'B': for both eigenvalues and eigenvectors (S and SEP).
                     64: *> \endverbatim
                     65: *>
                     66: *> \param[in] HOWMNY
                     67: *> \verbatim
                     68: *>          HOWMNY is CHARACTER*1
                     69: *>          = 'A': compute condition numbers for all eigenpairs;
                     70: *>          = 'S': compute condition numbers for selected eigenpairs
                     71: *>                 specified by the array SELECT.
                     72: *> \endverbatim
                     73: *>
                     74: *> \param[in] SELECT
                     75: *> \verbatim
                     76: *>          SELECT is LOGICAL array, dimension (N)
                     77: *>          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
                     78: *>          condition numbers are required. To select condition numbers
                     79: *>          for the eigenpair corresponding to a real eigenvalue w(j),
                     80: *>          SELECT(j) must be set to .TRUE.. To select condition numbers
                     81: *>          corresponding to a complex conjugate pair of eigenvalues w(j)
                     82: *>          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
                     83: *>          set to .TRUE..
                     84: *>          If HOWMNY = 'A', SELECT is not referenced.
                     85: *> \endverbatim
                     86: *>
                     87: *> \param[in] N
                     88: *> \verbatim
                     89: *>          N is INTEGER
                     90: *>          The order of the matrix T. N >= 0.
                     91: *> \endverbatim
                     92: *>
                     93: *> \param[in] T
                     94: *> \verbatim
                     95: *>          T is DOUBLE PRECISION array, dimension (LDT,N)
                     96: *>          The upper quasi-triangular matrix T, in Schur canonical form.
                     97: *> \endverbatim
                     98: *>
                     99: *> \param[in] LDT
                    100: *> \verbatim
                    101: *>          LDT is INTEGER
                    102: *>          The leading dimension of the array T. LDT >= max(1,N).
                    103: *> \endverbatim
                    104: *>
                    105: *> \param[in] VL
                    106: *> \verbatim
                    107: *>          VL is DOUBLE PRECISION array, dimension (LDVL,M)
                    108: *>          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
                    109: *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
                    110: *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
                    111: *>          must be stored in consecutive columns of VL, as returned by
                    112: *>          DHSEIN or DTREVC.
                    113: *>          If JOB = 'V', VL is not referenced.
                    114: *> \endverbatim
                    115: *>
                    116: *> \param[in] LDVL
                    117: *> \verbatim
                    118: *>          LDVL is INTEGER
                    119: *>          The leading dimension of the array VL.
                    120: *>          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
                    121: *> \endverbatim
                    122: *>
                    123: *> \param[in] VR
                    124: *> \verbatim
                    125: *>          VR is DOUBLE PRECISION array, dimension (LDVR,M)
                    126: *>          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
                    127: *>          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
                    128: *>          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
                    129: *>          must be stored in consecutive columns of VR, as returned by
                    130: *>          DHSEIN or DTREVC.
                    131: *>          If JOB = 'V', VR is not referenced.
                    132: *> \endverbatim
                    133: *>
                    134: *> \param[in] LDVR
                    135: *> \verbatim
                    136: *>          LDVR is INTEGER
                    137: *>          The leading dimension of the array VR.
                    138: *>          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
                    139: *> \endverbatim
                    140: *>
                    141: *> \param[out] S
                    142: *> \verbatim
                    143: *>          S is DOUBLE PRECISION array, dimension (MM)
                    144: *>          If JOB = 'E' or 'B', the reciprocal condition numbers of the
                    145: *>          selected eigenvalues, stored in consecutive elements of the
                    146: *>          array. For a complex conjugate pair of eigenvalues two
                    147: *>          consecutive elements of S are set to the same value. Thus
                    148: *>          S(j), SEP(j), and the j-th columns of VL and VR all
                    149: *>          correspond to the same eigenpair (but not in general the
                    150: *>          j-th eigenpair, unless all eigenpairs are selected).
                    151: *>          If JOB = 'V', S is not referenced.
                    152: *> \endverbatim
                    153: *>
                    154: *> \param[out] SEP
                    155: *> \verbatim
                    156: *>          SEP is DOUBLE PRECISION array, dimension (MM)
                    157: *>          If JOB = 'V' or 'B', the estimated reciprocal condition
                    158: *>          numbers of the selected eigenvectors, stored in consecutive
                    159: *>          elements of the array. For a complex eigenvector two
                    160: *>          consecutive elements of SEP are set to the same value. If
                    161: *>          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
                    162: *>          is set to 0; this can only occur when the true value would be
                    163: *>          very small anyway.
                    164: *>          If JOB = 'E', SEP is not referenced.
                    165: *> \endverbatim
                    166: *>
                    167: *> \param[in] MM
                    168: *> \verbatim
                    169: *>          MM is INTEGER
                    170: *>          The number of elements in the arrays S (if JOB = 'E' or 'B')
                    171: *>           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
                    172: *> \endverbatim
                    173: *>
                    174: *> \param[out] M
                    175: *> \verbatim
                    176: *>          M is INTEGER
                    177: *>          The number of elements of the arrays S and/or SEP actually
                    178: *>          used to store the estimated condition numbers.
                    179: *>          If HOWMNY = 'A', M is set to N.
                    180: *> \endverbatim
                    181: *>
                    182: *> \param[out] WORK
                    183: *> \verbatim
                    184: *>          WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)
                    185: *>          If JOB = 'E', WORK is not referenced.
                    186: *> \endverbatim
                    187: *>
                    188: *> \param[in] LDWORK
                    189: *> \verbatim
                    190: *>          LDWORK is INTEGER
                    191: *>          The leading dimension of the array WORK.
                    192: *>          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
                    193: *> \endverbatim
                    194: *>
                    195: *> \param[out] IWORK
                    196: *> \verbatim
                    197: *>          IWORK is INTEGER array, dimension (2*(N-1))
                    198: *>          If JOB = 'E', IWORK is not referenced.
                    199: *> \endverbatim
                    200: *>
                    201: *> \param[out] INFO
                    202: *> \verbatim
                    203: *>          INFO is INTEGER
                    204: *>          = 0: successful exit
                    205: *>          < 0: if INFO = -i, the i-th argument had an illegal value
                    206: *> \endverbatim
                    207: *
                    208: *  Authors:
                    209: *  ========
                    210: *
1.15      bertrand  211: *> \author Univ. of Tennessee
                    212: *> \author Univ. of California Berkeley
                    213: *> \author Univ. of Colorado Denver
                    214: *> \author NAG Ltd.
1.9       bertrand  215: *
                    216: *> \ingroup doubleOTHERcomputational
                    217: *
                    218: *> \par Further Details:
                    219: *  =====================
                    220: *>
                    221: *> \verbatim
                    222: *>
                    223: *>  The reciprocal of the condition number of an eigenvalue lambda is
                    224: *>  defined as
                    225: *>
                    226: *>          S(lambda) = |v**T*u| / (norm(u)*norm(v))
                    227: *>
                    228: *>  where u and v are the right and left eigenvectors of T corresponding
                    229: *>  to lambda; v**T denotes the transpose of v, and norm(u)
                    230: *>  denotes the Euclidean norm. These reciprocal condition numbers always
                    231: *>  lie between zero (very badly conditioned) and one (very well
                    232: *>  conditioned). If n = 1, S(lambda) is defined to be 1.
                    233: *>
                    234: *>  An approximate error bound for a computed eigenvalue W(i) is given by
                    235: *>
                    236: *>                      EPS * norm(T) / S(i)
                    237: *>
                    238: *>  where EPS is the machine precision.
                    239: *>
                    240: *>  The reciprocal of the condition number of the right eigenvector u
                    241: *>  corresponding to lambda is defined as follows. Suppose
                    242: *>
                    243: *>              T = ( lambda  c  )
                    244: *>                  (   0    T22 )
                    245: *>
                    246: *>  Then the reciprocal condition number is
                    247: *>
                    248: *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
                    249: *>
                    250: *>  where sigma-min denotes the smallest singular value. We approximate
                    251: *>  the smallest singular value by the reciprocal of an estimate of the
                    252: *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
                    253: *>  defined to be abs(T(1,1)).
                    254: *>
                    255: *>  An approximate error bound for a computed right eigenvector VR(i)
                    256: *>  is given by
                    257: *>
                    258: *>                      EPS * norm(T) / SEP(i)
                    259: *> \endverbatim
                    260: *>
                    261: *  =====================================================================
1.1       bertrand  262:       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
                    263:      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
                    264:      $                   INFO )
                    265: *
1.19    ! bertrand  266: *  -- LAPACK computational routine --
1.1       bertrand  267: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    268: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    269: *
                    270: *     .. Scalar Arguments ..
                    271:       CHARACTER          HOWMNY, JOB
                    272:       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
                    273: *     ..
                    274: *     .. Array Arguments ..
                    275:       LOGICAL            SELECT( * )
                    276:       INTEGER            IWORK( * )
                    277:       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
                    278:      $                   VR( LDVR, * ), WORK( LDWORK, * )
                    279: *     ..
                    280: *
                    281: *  =====================================================================
                    282: *
                    283: *     .. Parameters ..
                    284:       DOUBLE PRECISION   ZERO, ONE, TWO
                    285:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
                    286: *     ..
                    287: *     .. Local Scalars ..
                    288:       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
                    289:       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
                    290:       DOUBLE PRECISION   BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
                    291:      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
                    292: *     ..
                    293: *     .. Local Arrays ..
                    294:       INTEGER            ISAVE( 3 )
                    295:       DOUBLE PRECISION   DUMMY( 1 )
                    296: *     ..
                    297: *     .. External Functions ..
                    298:       LOGICAL            LSAME
                    299:       DOUBLE PRECISION   DDOT, DLAMCH, DLAPY2, DNRM2
                    300:       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
                    301: *     ..
                    302: *     .. External Subroutines ..
1.17      bertrand  303:       EXTERNAL           DLABAD, DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
1.1       bertrand  304: *     ..
                    305: *     .. Intrinsic Functions ..
                    306:       INTRINSIC          ABS, MAX, SQRT
                    307: *     ..
                    308: *     .. Executable Statements ..
                    309: *
                    310: *     Decode and test the input parameters
                    311: *
                    312:       WANTBH = LSAME( JOB, 'B' )
                    313:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
                    314:       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
                    315: *
                    316:       SOMCON = LSAME( HOWMNY, 'S' )
                    317: *
                    318:       INFO = 0
                    319:       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
                    320:          INFO = -1
                    321:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
                    322:          INFO = -2
                    323:       ELSE IF( N.LT.0 ) THEN
                    324:          INFO = -4
                    325:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
                    326:          INFO = -6
                    327:       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
                    328:          INFO = -8
                    329:       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
                    330:          INFO = -10
                    331:       ELSE
                    332: *
                    333: *        Set M to the number of eigenpairs for which condition numbers
                    334: *        are required, and test MM.
                    335: *
                    336:          IF( SOMCON ) THEN
                    337:             M = 0
                    338:             PAIR = .FALSE.
                    339:             DO 10 K = 1, N
                    340:                IF( PAIR ) THEN
                    341:                   PAIR = .FALSE.
                    342:                ELSE
                    343:                   IF( K.LT.N ) THEN
                    344:                      IF( T( K+1, K ).EQ.ZERO ) THEN
                    345:                         IF( SELECT( K ) )
                    346:      $                     M = M + 1
                    347:                      ELSE
                    348:                         PAIR = .TRUE.
                    349:                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
                    350:      $                     M = M + 2
                    351:                      END IF
                    352:                   ELSE
                    353:                      IF( SELECT( N ) )
                    354:      $                  M = M + 1
                    355:                   END IF
                    356:                END IF
                    357:    10       CONTINUE
                    358:          ELSE
                    359:             M = N
                    360:          END IF
                    361: *
                    362:          IF( MM.LT.M ) THEN
                    363:             INFO = -13
                    364:          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
                    365:             INFO = -16
                    366:          END IF
                    367:       END IF
                    368:       IF( INFO.NE.0 ) THEN
                    369:          CALL XERBLA( 'DTRSNA', -INFO )
                    370:          RETURN
                    371:       END IF
                    372: *
                    373: *     Quick return if possible
                    374: *
                    375:       IF( N.EQ.0 )
                    376:      $   RETURN
                    377: *
                    378:       IF( N.EQ.1 ) THEN
                    379:          IF( SOMCON ) THEN
                    380:             IF( .NOT.SELECT( 1 ) )
                    381:      $         RETURN
                    382:          END IF
                    383:          IF( WANTS )
                    384:      $      S( 1 ) = ONE
                    385:          IF( WANTSP )
                    386:      $      SEP( 1 ) = ABS( T( 1, 1 ) )
                    387:          RETURN
                    388:       END IF
                    389: *
                    390: *     Get machine constants
                    391: *
                    392:       EPS = DLAMCH( 'P' )
                    393:       SMLNUM = DLAMCH( 'S' ) / EPS
                    394:       BIGNUM = ONE / SMLNUM
                    395:       CALL DLABAD( SMLNUM, BIGNUM )
                    396: *
                    397:       KS = 0
                    398:       PAIR = .FALSE.
                    399:       DO 60 K = 1, N
                    400: *
                    401: *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
                    402: *
                    403:          IF( PAIR ) THEN
                    404:             PAIR = .FALSE.
                    405:             GO TO 60
                    406:          ELSE
                    407:             IF( K.LT.N )
                    408:      $         PAIR = T( K+1, K ).NE.ZERO
                    409:          END IF
                    410: *
                    411: *        Determine whether condition numbers are required for the k-th
                    412: *        eigenpair.
                    413: *
                    414:          IF( SOMCON ) THEN
                    415:             IF( PAIR ) THEN
                    416:                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
                    417:      $            GO TO 60
                    418:             ELSE
                    419:                IF( .NOT.SELECT( K ) )
                    420:      $            GO TO 60
                    421:             END IF
                    422:          END IF
                    423: *
                    424:          KS = KS + 1
                    425: *
                    426:          IF( WANTS ) THEN
                    427: *
                    428: *           Compute the reciprocal condition number of the k-th
                    429: *           eigenvalue.
                    430: *
                    431:             IF( .NOT.PAIR ) THEN
                    432: *
                    433: *              Real eigenvalue.
                    434: *
                    435:                PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
                    436:                RNRM = DNRM2( N, VR( 1, KS ), 1 )
                    437:                LNRM = DNRM2( N, VL( 1, KS ), 1 )
                    438:                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
                    439:             ELSE
                    440: *
                    441: *              Complex eigenvalue.
                    442: *
                    443:                PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
                    444:                PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
                    445:      $                 1 )
                    446:                PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
                    447:                PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
                    448:      $                 1 )
                    449:                RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
                    450:      $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
                    451:                LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
                    452:      $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
                    453:                COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
                    454:                S( KS ) = COND
                    455:                S( KS+1 ) = COND
                    456:             END IF
                    457:          END IF
                    458: *
                    459:          IF( WANTSP ) THEN
                    460: *
                    461: *           Estimate the reciprocal condition number of the k-th
                    462: *           eigenvector.
                    463: *
                    464: *           Copy the matrix T to the array WORK and swap the diagonal
                    465: *           block beginning at T(k,k) to the (1,1) position.
                    466: *
                    467:             CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
                    468:             IFST = K
                    469:             ILST = 1
                    470:             CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
                    471:      $                   WORK( 1, N+1 ), IERR )
                    472: *
                    473:             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
                    474: *
                    475: *              Could not swap because blocks not well separated
                    476: *
                    477:                SCALE = ONE
                    478:                EST = BIGNUM
                    479:             ELSE
                    480: *
                    481: *              Reordering successful
                    482: *
                    483:                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
                    484: *
                    485: *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
                    486: *
                    487:                   DO 20 I = 2, N
                    488:                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
                    489:    20             CONTINUE
                    490:                   N2 = 1
                    491:                   NN = N - 1
                    492:                ELSE
                    493: *
                    494: *                 Triangularize the 2 by 2 block by unitary
                    495: *                 transformation U = [  cs   i*ss ]
                    496: *                                    [ i*ss   cs  ].
                    497: *                 such that the (1,1) position of WORK is complex
                    498: *                 eigenvalue lambda with positive imaginary part. (2,2)
                    499: *                 position of WORK is the complex eigenvalue lambda
                    500: *                 with negative imaginary  part.
                    501: *
                    502:                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
                    503:      $                 SQRT( ABS( WORK( 2, 1 ) ) )
                    504:                   DELTA = DLAPY2( MU, WORK( 2, 1 ) )
                    505:                   CS = MU / DELTA
                    506:                   SN = -WORK( 2, 1 ) / DELTA
                    507: *
                    508: *                 Form
                    509: *
1.8       bertrand  510: *                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
                    511: *                                          [   mu                     ]
                    512: *                                          [         ..               ]
                    513: *                                          [             ..           ]
                    514: *                                          [                  mu      ]
                    515: *                 where C**T is transpose of matrix C,
1.1       bertrand  516: *                 and RWORK is stored starting in the N+1-st column of
                    517: *                 WORK.
                    518: *
                    519:                   DO 30 J = 3, N
                    520:                      WORK( 2, J ) = CS*WORK( 2, J )
                    521:                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
                    522:    30             CONTINUE
                    523:                   WORK( 2, 2 ) = ZERO
                    524: *
                    525:                   WORK( 1, N+1 ) = TWO*MU
                    526:                   DO 40 I = 2, N - 1
                    527:                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
                    528:    40             CONTINUE
                    529:                   N2 = 2
                    530:                   NN = 2*( N-1 )
                    531:                END IF
                    532: *
1.8       bertrand  533: *              Estimate norm(inv(C**T))
1.1       bertrand  534: *
                    535:                EST = ZERO
                    536:                KASE = 0
                    537:    50          CONTINUE
                    538:                CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
                    539:      $                      EST, KASE, ISAVE )
                    540:                IF( KASE.NE.0 ) THEN
                    541:                   IF( KASE.EQ.1 ) THEN
                    542:                      IF( N2.EQ.1 ) THEN
                    543: *
1.8       bertrand  544: *                       Real eigenvalue: solve C**T*x = scale*c.
1.1       bertrand  545: *
                    546:                         CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
                    547:      $                               LDWORK, DUMMY, DUMM, SCALE,
                    548:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    549:      $                               IERR )
                    550:                      ELSE
                    551: *
                    552: *                       Complex eigenvalue: solve
1.8       bertrand  553: *                       C**T*(p+iq) = scale*(c+id) in real arithmetic.
1.1       bertrand  554: *
                    555:                         CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
                    556:      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
                    557:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    558:      $                               IERR )
                    559:                      END IF
                    560:                   ELSE
                    561:                      IF( N2.EQ.1 ) THEN
                    562: *
                    563: *                       Real eigenvalue: solve C*x = scale*c.
                    564: *
                    565:                         CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
                    566:      $                               LDWORK, DUMMY, DUMM, SCALE,
                    567:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    568:      $                               IERR )
                    569:                      ELSE
                    570: *
                    571: *                       Complex eigenvalue: solve
                    572: *                       C*(p+iq) = scale*(c+id) in real arithmetic.
                    573: *
                    574:                         CALL DLAQTR( .FALSE., .FALSE., N-1,
                    575:      $                               WORK( 2, 2 ), LDWORK,
                    576:      $                               WORK( 1, N+1 ), MU, SCALE,
                    577:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    578:      $                               IERR )
                    579: *
                    580:                      END IF
                    581:                   END IF
                    582: *
                    583:                   GO TO 50
                    584:                END IF
                    585:             END IF
                    586: *
                    587:             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
                    588:             IF( PAIR )
                    589:      $         SEP( KS+1 ) = SEP( KS )
                    590:          END IF
                    591: *
                    592:          IF( PAIR )
                    593:      $      KS = KS + 1
                    594: *
                    595:    60 CONTINUE
                    596:       RETURN
                    597: *
                    598: *     End of DTRSNA
                    599: *
                    600:       END

CVSweb interface <joel.bertrand@systella.fr>