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

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: *
1.15      bertrand  216: *> \date December 2016
1.9       bertrand  217: *
                    218: *> \ingroup doubleOTHERcomputational
                    219: *
                    220: *> \par Further Details:
                    221: *  =====================
                    222: *>
                    223: *> \verbatim
                    224: *>
                    225: *>  The reciprocal of the condition number of an eigenvalue lambda is
                    226: *>  defined as
                    227: *>
                    228: *>          S(lambda) = |v**T*u| / (norm(u)*norm(v))
                    229: *>
                    230: *>  where u and v are the right and left eigenvectors of T corresponding
                    231: *>  to lambda; v**T denotes the transpose of v, and norm(u)
                    232: *>  denotes the Euclidean norm. These reciprocal condition numbers always
                    233: *>  lie between zero (very badly conditioned) and one (very well
                    234: *>  conditioned). If n = 1, S(lambda) is defined to be 1.
                    235: *>
                    236: *>  An approximate error bound for a computed eigenvalue W(i) is given by
                    237: *>
                    238: *>                      EPS * norm(T) / S(i)
                    239: *>
                    240: *>  where EPS is the machine precision.
                    241: *>
                    242: *>  The reciprocal of the condition number of the right eigenvector u
                    243: *>  corresponding to lambda is defined as follows. Suppose
                    244: *>
                    245: *>              T = ( lambda  c  )
                    246: *>                  (   0    T22 )
                    247: *>
                    248: *>  Then the reciprocal condition number is
                    249: *>
                    250: *>          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
                    251: *>
                    252: *>  where sigma-min denotes the smallest singular value. We approximate
                    253: *>  the smallest singular value by the reciprocal of an estimate of the
                    254: *>  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
                    255: *>  defined to be abs(T(1,1)).
                    256: *>
                    257: *>  An approximate error bound for a computed right eigenvector VR(i)
                    258: *>  is given by
                    259: *>
                    260: *>                      EPS * norm(T) / SEP(i)
                    261: *> \endverbatim
                    262: *>
                    263: *  =====================================================================
1.1       bertrand  264:       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
                    265:      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
                    266:      $                   INFO )
                    267: *
1.15      bertrand  268: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  269: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    270: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15      bertrand  271: *     December 2016
1.1       bertrand  272: *
                    273: *     .. Scalar Arguments ..
                    274:       CHARACTER          HOWMNY, JOB
                    275:       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
                    276: *     ..
                    277: *     .. Array Arguments ..
                    278:       LOGICAL            SELECT( * )
                    279:       INTEGER            IWORK( * )
                    280:       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
                    281:      $                   VR( LDVR, * ), WORK( LDWORK, * )
                    282: *     ..
                    283: *
                    284: *  =====================================================================
                    285: *
                    286: *     .. Parameters ..
                    287:       DOUBLE PRECISION   ZERO, ONE, TWO
                    288:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
                    289: *     ..
                    290: *     .. Local Scalars ..
                    291:       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
                    292:       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
                    293:       DOUBLE PRECISION   BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
                    294:      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
                    295: *     ..
                    296: *     .. Local Arrays ..
                    297:       INTEGER            ISAVE( 3 )
                    298:       DOUBLE PRECISION   DUMMY( 1 )
                    299: *     ..
                    300: *     .. External Functions ..
                    301:       LOGICAL            LSAME
                    302:       DOUBLE PRECISION   DDOT, DLAMCH, DLAPY2, DNRM2
                    303:       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
                    304: *     ..
                    305: *     .. External Subroutines ..
                    306:       EXTERNAL           DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
                    307: *     ..
                    308: *     .. Intrinsic Functions ..
                    309:       INTRINSIC          ABS, MAX, SQRT
                    310: *     ..
                    311: *     .. Executable Statements ..
                    312: *
                    313: *     Decode and test the input parameters
                    314: *
                    315:       WANTBH = LSAME( JOB, 'B' )
                    316:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
                    317:       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
                    318: *
                    319:       SOMCON = LSAME( HOWMNY, 'S' )
                    320: *
                    321:       INFO = 0
                    322:       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
                    323:          INFO = -1
                    324:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
                    325:          INFO = -2
                    326:       ELSE IF( N.LT.0 ) THEN
                    327:          INFO = -4
                    328:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
                    329:          INFO = -6
                    330:       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
                    331:          INFO = -8
                    332:       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
                    333:          INFO = -10
                    334:       ELSE
                    335: *
                    336: *        Set M to the number of eigenpairs for which condition numbers
                    337: *        are required, and test MM.
                    338: *
                    339:          IF( SOMCON ) THEN
                    340:             M = 0
                    341:             PAIR = .FALSE.
                    342:             DO 10 K = 1, N
                    343:                IF( PAIR ) THEN
                    344:                   PAIR = .FALSE.
                    345:                ELSE
                    346:                   IF( K.LT.N ) THEN
                    347:                      IF( T( K+1, K ).EQ.ZERO ) THEN
                    348:                         IF( SELECT( K ) )
                    349:      $                     M = M + 1
                    350:                      ELSE
                    351:                         PAIR = .TRUE.
                    352:                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
                    353:      $                     M = M + 2
                    354:                      END IF
                    355:                   ELSE
                    356:                      IF( SELECT( N ) )
                    357:      $                  M = M + 1
                    358:                   END IF
                    359:                END IF
                    360:    10       CONTINUE
                    361:          ELSE
                    362:             M = N
                    363:          END IF
                    364: *
                    365:          IF( MM.LT.M ) THEN
                    366:             INFO = -13
                    367:          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
                    368:             INFO = -16
                    369:          END IF
                    370:       END IF
                    371:       IF( INFO.NE.0 ) THEN
                    372:          CALL XERBLA( 'DTRSNA', -INFO )
                    373:          RETURN
                    374:       END IF
                    375: *
                    376: *     Quick return if possible
                    377: *
                    378:       IF( N.EQ.0 )
                    379:      $   RETURN
                    380: *
                    381:       IF( N.EQ.1 ) THEN
                    382:          IF( SOMCON ) THEN
                    383:             IF( .NOT.SELECT( 1 ) )
                    384:      $         RETURN
                    385:          END IF
                    386:          IF( WANTS )
                    387:      $      S( 1 ) = ONE
                    388:          IF( WANTSP )
                    389:      $      SEP( 1 ) = ABS( T( 1, 1 ) )
                    390:          RETURN
                    391:       END IF
                    392: *
                    393: *     Get machine constants
                    394: *
                    395:       EPS = DLAMCH( 'P' )
                    396:       SMLNUM = DLAMCH( 'S' ) / EPS
                    397:       BIGNUM = ONE / SMLNUM
                    398:       CALL DLABAD( SMLNUM, BIGNUM )
                    399: *
                    400:       KS = 0
                    401:       PAIR = .FALSE.
                    402:       DO 60 K = 1, N
                    403: *
                    404: *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
                    405: *
                    406:          IF( PAIR ) THEN
                    407:             PAIR = .FALSE.
                    408:             GO TO 60
                    409:          ELSE
                    410:             IF( K.LT.N )
                    411:      $         PAIR = T( K+1, K ).NE.ZERO
                    412:          END IF
                    413: *
                    414: *        Determine whether condition numbers are required for the k-th
                    415: *        eigenpair.
                    416: *
                    417:          IF( SOMCON ) THEN
                    418:             IF( PAIR ) THEN
                    419:                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
                    420:      $            GO TO 60
                    421:             ELSE
                    422:                IF( .NOT.SELECT( K ) )
                    423:      $            GO TO 60
                    424:             END IF
                    425:          END IF
                    426: *
                    427:          KS = KS + 1
                    428: *
                    429:          IF( WANTS ) THEN
                    430: *
                    431: *           Compute the reciprocal condition number of the k-th
                    432: *           eigenvalue.
                    433: *
                    434:             IF( .NOT.PAIR ) THEN
                    435: *
                    436: *              Real eigenvalue.
                    437: *
                    438:                PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
                    439:                RNRM = DNRM2( N, VR( 1, KS ), 1 )
                    440:                LNRM = DNRM2( N, VL( 1, KS ), 1 )
                    441:                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
                    442:             ELSE
                    443: *
                    444: *              Complex eigenvalue.
                    445: *
                    446:                PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
                    447:                PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
                    448:      $                 1 )
                    449:                PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
                    450:                PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
                    451:      $                 1 )
                    452:                RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
                    453:      $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
                    454:                LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
                    455:      $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
                    456:                COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
                    457:                S( KS ) = COND
                    458:                S( KS+1 ) = COND
                    459:             END IF
                    460:          END IF
                    461: *
                    462:          IF( WANTSP ) THEN
                    463: *
                    464: *           Estimate the reciprocal condition number of the k-th
                    465: *           eigenvector.
                    466: *
                    467: *           Copy the matrix T to the array WORK and swap the diagonal
                    468: *           block beginning at T(k,k) to the (1,1) position.
                    469: *
                    470:             CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
                    471:             IFST = K
                    472:             ILST = 1
                    473:             CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
                    474:      $                   WORK( 1, N+1 ), IERR )
                    475: *
                    476:             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
                    477: *
                    478: *              Could not swap because blocks not well separated
                    479: *
                    480:                SCALE = ONE
                    481:                EST = BIGNUM
                    482:             ELSE
                    483: *
                    484: *              Reordering successful
                    485: *
                    486:                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
                    487: *
                    488: *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
                    489: *
                    490:                   DO 20 I = 2, N
                    491:                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
                    492:    20             CONTINUE
                    493:                   N2 = 1
                    494:                   NN = N - 1
                    495:                ELSE
                    496: *
                    497: *                 Triangularize the 2 by 2 block by unitary
                    498: *                 transformation U = [  cs   i*ss ]
                    499: *                                    [ i*ss   cs  ].
                    500: *                 such that the (1,1) position of WORK is complex
                    501: *                 eigenvalue lambda with positive imaginary part. (2,2)
                    502: *                 position of WORK is the complex eigenvalue lambda
                    503: *                 with negative imaginary  part.
                    504: *
                    505:                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
                    506:      $                 SQRT( ABS( WORK( 2, 1 ) ) )
                    507:                   DELTA = DLAPY2( MU, WORK( 2, 1 ) )
                    508:                   CS = MU / DELTA
                    509:                   SN = -WORK( 2, 1 ) / DELTA
                    510: *
                    511: *                 Form
                    512: *
1.8       bertrand  513: *                 C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
                    514: *                                          [   mu                     ]
                    515: *                                          [         ..               ]
                    516: *                                          [             ..           ]
                    517: *                                          [                  mu      ]
                    518: *                 where C**T is transpose of matrix C,
1.1       bertrand  519: *                 and RWORK is stored starting in the N+1-st column of
                    520: *                 WORK.
                    521: *
                    522:                   DO 30 J = 3, N
                    523:                      WORK( 2, J ) = CS*WORK( 2, J )
                    524:                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
                    525:    30             CONTINUE
                    526:                   WORK( 2, 2 ) = ZERO
                    527: *
                    528:                   WORK( 1, N+1 ) = TWO*MU
                    529:                   DO 40 I = 2, N - 1
                    530:                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
                    531:    40             CONTINUE
                    532:                   N2 = 2
                    533:                   NN = 2*( N-1 )
                    534:                END IF
                    535: *
1.8       bertrand  536: *              Estimate norm(inv(C**T))
1.1       bertrand  537: *
                    538:                EST = ZERO
                    539:                KASE = 0
                    540:    50          CONTINUE
                    541:                CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
                    542:      $                      EST, KASE, ISAVE )
                    543:                IF( KASE.NE.0 ) THEN
                    544:                   IF( KASE.EQ.1 ) THEN
                    545:                      IF( N2.EQ.1 ) THEN
                    546: *
1.8       bertrand  547: *                       Real eigenvalue: solve C**T*x = scale*c.
1.1       bertrand  548: *
                    549:                         CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
                    550:      $                               LDWORK, DUMMY, DUMM, SCALE,
                    551:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    552:      $                               IERR )
                    553:                      ELSE
                    554: *
                    555: *                       Complex eigenvalue: solve
1.8       bertrand  556: *                       C**T*(p+iq) = scale*(c+id) in real arithmetic.
1.1       bertrand  557: *
                    558:                         CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
                    559:      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
                    560:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    561:      $                               IERR )
                    562:                      END IF
                    563:                   ELSE
                    564:                      IF( N2.EQ.1 ) THEN
                    565: *
                    566: *                       Real eigenvalue: solve C*x = scale*c.
                    567: *
                    568:                         CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
                    569:      $                               LDWORK, DUMMY, DUMM, SCALE,
                    570:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    571:      $                               IERR )
                    572:                      ELSE
                    573: *
                    574: *                       Complex eigenvalue: solve
                    575: *                       C*(p+iq) = scale*(c+id) in real arithmetic.
                    576: *
                    577:                         CALL DLAQTR( .FALSE., .FALSE., N-1,
                    578:      $                               WORK( 2, 2 ), LDWORK,
                    579:      $                               WORK( 1, N+1 ), MU, SCALE,
                    580:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
                    581:      $                               IERR )
                    582: *
                    583:                      END IF
                    584:                   END IF
                    585: *
                    586:                   GO TO 50
                    587:                END IF
                    588:             END IF
                    589: *
                    590:             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
                    591:             IF( PAIR )
                    592:      $         SEP( KS+1 ) = SEP( KS )
                    593:          END IF
                    594: *
                    595:          IF( PAIR )
                    596:      $      KS = KS + 1
                    597: *
                    598:    60 CONTINUE
                    599:       RETURN
                    600: *
                    601: *     End of DTRSNA
                    602: *
                    603:       END

CVSweb interface <joel.bertrand@systella.fr>