Annotation of rpl/lapack/lapack/dlaein.f, revision 1.3

1.1       bertrand    1:       SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
                      2:      $                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
                      3: *
                      4: *  -- LAPACK auxiliary 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            NOINIT, RIGHTV
                     11:       INTEGER            INFO, LDB, LDH, N
                     12:       DOUBLE PRECISION   BIGNUM, EPS3, SMLNUM, WI, WR
                     13: *     ..
                     14: *     .. Array Arguments ..
                     15:       DOUBLE PRECISION   B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
                     16:      $                   WORK( * )
                     17: *     ..
                     18: *
                     19: *  Purpose
                     20: *  =======
                     21: *
                     22: *  DLAEIN uses inverse iteration to find a right or left eigenvector
                     23: *  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
                     24: *  matrix H.
                     25: *
                     26: *  Arguments
                     27: *  =========
                     28: *
                     29: *  RIGHTV   (input) LOGICAL
                     30: *          = .TRUE. : compute right eigenvector;
                     31: *          = .FALSE.: compute left eigenvector.
                     32: *
                     33: *  NOINIT   (input) LOGICAL
                     34: *          = .TRUE. : no initial vector supplied in (VR,VI).
                     35: *          = .FALSE.: initial vector supplied in (VR,VI).
                     36: *
                     37: *  N       (input) INTEGER
                     38: *          The order of the matrix H.  N >= 0.
                     39: *
                     40: *  H       (input) DOUBLE PRECISION array, dimension (LDH,N)
                     41: *          The upper Hessenberg matrix H.
                     42: *
                     43: *  LDH     (input) INTEGER
                     44: *          The leading dimension of the array H.  LDH >= max(1,N).
                     45: *
                     46: *  WR      (input) DOUBLE PRECISION
                     47: *  WI      (input) DOUBLE PRECISION
                     48: *          The real and imaginary parts of the eigenvalue of H whose
                     49: *          corresponding right or left eigenvector is to be computed.
                     50: *
                     51: *  VR      (input/output) DOUBLE PRECISION array, dimension (N)
                     52: *  VI      (input/output) DOUBLE PRECISION array, dimension (N)
                     53: *          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
                     54: *          a real starting vector for inverse iteration using the real
                     55: *          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
                     56: *          must contain the real and imaginary parts of a complex
                     57: *          starting vector for inverse iteration using the complex
                     58: *          eigenvalue (WR,WI); otherwise VR and VI need not be set.
                     59: *          On exit, if WI = 0.0 (real eigenvalue), VR contains the
                     60: *          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
                     61: *          VR and VI contain the real and imaginary parts of the
                     62: *          computed complex eigenvector. The eigenvector is normalized
                     63: *          so that the component of largest magnitude has magnitude 1;
                     64: *          here the magnitude of a complex number (x,y) is taken to be
                     65: *          |x| + |y|.
                     66: *          VI is not referenced if WI = 0.0.
                     67: *
                     68: *  B       (workspace) DOUBLE PRECISION array, dimension (LDB,N)
                     69: *
                     70: *  LDB     (input) INTEGER
                     71: *          The leading dimension of the array B.  LDB >= N+1.
                     72: *
                     73: *  WORK   (workspace) DOUBLE PRECISION array, dimension (N)
                     74: *
                     75: *  EPS3    (input) DOUBLE PRECISION
                     76: *          A small machine-dependent value which is used to perturb
                     77: *          close eigenvalues, and to replace zero pivots.
                     78: *
                     79: *  SMLNUM  (input) DOUBLE PRECISION
                     80: *          A machine-dependent value close to the underflow threshold.
                     81: *
                     82: *  BIGNUM  (input) DOUBLE PRECISION
                     83: *          A machine-dependent value close to the overflow threshold.
                     84: *
                     85: *  INFO    (output) INTEGER
                     86: *          = 0:  successful exit
                     87: *          = 1:  inverse iteration did not converge; VR is set to the
                     88: *                last iterate, and so is VI if WI.ne.0.0.
                     89: *
                     90: *  =====================================================================
                     91: *
                     92: *     .. Parameters ..
                     93:       DOUBLE PRECISION   ZERO, ONE, TENTH
                     94:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
                     95: *     ..
                     96: *     .. Local Scalars ..
                     97:       CHARACTER          NORMIN, TRANS
                     98:       INTEGER            I, I1, I2, I3, IERR, ITS, J
                     99:       DOUBLE PRECISION   ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
                    100:      $                   REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
                    101:      $                   W1, X, XI, XR, Y
                    102: *     ..
                    103: *     .. External Functions ..
                    104:       INTEGER            IDAMAX
                    105:       DOUBLE PRECISION   DASUM, DLAPY2, DNRM2
                    106:       EXTERNAL           IDAMAX, DASUM, DLAPY2, DNRM2
                    107: *     ..
                    108: *     .. External Subroutines ..
                    109:       EXTERNAL           DLADIV, DLATRS, DSCAL
                    110: *     ..
                    111: *     .. Intrinsic Functions ..
                    112:       INTRINSIC          ABS, DBLE, MAX, SQRT
                    113: *     ..
                    114: *     .. Executable Statements ..
                    115: *
                    116:       INFO = 0
                    117: *
                    118: *     GROWTO is the threshold used in the acceptance test for an
                    119: *     eigenvector.
                    120: *
                    121:       ROOTN = SQRT( DBLE( N ) )
                    122:       GROWTO = TENTH / ROOTN
                    123:       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
                    124: *
                    125: *     Form B = H - (WR,WI)*I (except that the subdiagonal elements and
                    126: *     the imaginary parts of the diagonal elements are not stored).
                    127: *
                    128:       DO 20 J = 1, N
                    129:          DO 10 I = 1, J - 1
                    130:             B( I, J ) = H( I, J )
                    131:    10    CONTINUE
                    132:          B( J, J ) = H( J, J ) - WR
                    133:    20 CONTINUE
                    134: *
                    135:       IF( WI.EQ.ZERO ) THEN
                    136: *
                    137: *        Real eigenvalue.
                    138: *
                    139:          IF( NOINIT ) THEN
                    140: *
                    141: *           Set initial vector.
                    142: *
                    143:             DO 30 I = 1, N
                    144:                VR( I ) = EPS3
                    145:    30       CONTINUE
                    146:          ELSE
                    147: *
                    148: *           Scale supplied initial vector.
                    149: *
                    150:             VNORM = DNRM2( N, VR, 1 )
                    151:             CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
                    152:      $                  1 )
                    153:          END IF
                    154: *
                    155:          IF( RIGHTV ) THEN
                    156: *
                    157: *           LU decomposition with partial pivoting of B, replacing zero
                    158: *           pivots by EPS3.
                    159: *
                    160:             DO 60 I = 1, N - 1
                    161:                EI = H( I+1, I )
                    162:                IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
                    163: *
                    164: *                 Interchange rows and eliminate.
                    165: *
                    166:                   X = B( I, I ) / EI
                    167:                   B( I, I ) = EI
                    168:                   DO 40 J = I + 1, N
                    169:                      TEMP = B( I+1, J )
                    170:                      B( I+1, J ) = B( I, J ) - X*TEMP
                    171:                      B( I, J ) = TEMP
                    172:    40             CONTINUE
                    173:                ELSE
                    174: *
                    175: *                 Eliminate without interchange.
                    176: *
                    177:                   IF( B( I, I ).EQ.ZERO )
                    178:      $               B( I, I ) = EPS3
                    179:                   X = EI / B( I, I )
                    180:                   IF( X.NE.ZERO ) THEN
                    181:                      DO 50 J = I + 1, N
                    182:                         B( I+1, J ) = B( I+1, J ) - X*B( I, J )
                    183:    50                CONTINUE
                    184:                   END IF
                    185:                END IF
                    186:    60       CONTINUE
                    187:             IF( B( N, N ).EQ.ZERO )
                    188:      $         B( N, N ) = EPS3
                    189: *
                    190:             TRANS = 'N'
                    191: *
                    192:          ELSE
                    193: *
                    194: *           UL decomposition with partial pivoting of B, replacing zero
                    195: *           pivots by EPS3.
                    196: *
                    197:             DO 90 J = N, 2, -1
                    198:                EJ = H( J, J-1 )
                    199:                IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
                    200: *
                    201: *                 Interchange columns and eliminate.
                    202: *
                    203:                   X = B( J, J ) / EJ
                    204:                   B( J, J ) = EJ
                    205:                   DO 70 I = 1, J - 1
                    206:                      TEMP = B( I, J-1 )
                    207:                      B( I, J-1 ) = B( I, J ) - X*TEMP
                    208:                      B( I, J ) = TEMP
                    209:    70             CONTINUE
                    210:                ELSE
                    211: *
                    212: *                 Eliminate without interchange.
                    213: *
                    214:                   IF( B( J, J ).EQ.ZERO )
                    215:      $               B( J, J ) = EPS3
                    216:                   X = EJ / B( J, J )
                    217:                   IF( X.NE.ZERO ) THEN
                    218:                      DO 80 I = 1, J - 1
                    219:                         B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
                    220:    80                CONTINUE
                    221:                   END IF
                    222:                END IF
                    223:    90       CONTINUE
                    224:             IF( B( 1, 1 ).EQ.ZERO )
                    225:      $         B( 1, 1 ) = EPS3
                    226: *
                    227:             TRANS = 'T'
                    228: *
                    229:          END IF
                    230: *
                    231:          NORMIN = 'N'
                    232:          DO 110 ITS = 1, N
                    233: *
                    234: *           Solve U*x = scale*v for a right eigenvector
                    235: *             or U'*x = scale*v for a left eigenvector,
                    236: *           overwriting x on v.
                    237: *
                    238:             CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
                    239:      $                   VR, SCALE, WORK, IERR )
                    240:             NORMIN = 'Y'
                    241: *
                    242: *           Test for sufficient growth in the norm of v.
                    243: *
                    244:             VNORM = DASUM( N, VR, 1 )
                    245:             IF( VNORM.GE.GROWTO*SCALE )
                    246:      $         GO TO 120
                    247: *
                    248: *           Choose new orthogonal starting vector and try again.
                    249: *
                    250:             TEMP = EPS3 / ( ROOTN+ONE )
                    251:             VR( 1 ) = EPS3
                    252:             DO 100 I = 2, N
                    253:                VR( I ) = TEMP
                    254:   100       CONTINUE
                    255:             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
                    256:   110    CONTINUE
                    257: *
                    258: *        Failure to find eigenvector in N iterations.
                    259: *
                    260:          INFO = 1
                    261: *
                    262:   120    CONTINUE
                    263: *
                    264: *        Normalize eigenvector.
                    265: *
                    266:          I = IDAMAX( N, VR, 1 )
                    267:          CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
                    268:       ELSE
                    269: *
                    270: *        Complex eigenvalue.
                    271: *
                    272:          IF( NOINIT ) THEN
                    273: *
                    274: *           Set initial vector.
                    275: *
                    276:             DO 130 I = 1, N
                    277:                VR( I ) = EPS3
                    278:                VI( I ) = ZERO
                    279:   130       CONTINUE
                    280:          ELSE
                    281: *
                    282: *           Scale supplied initial vector.
                    283: *
                    284:             NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
                    285:             REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
                    286:             CALL DSCAL( N, REC, VR, 1 )
                    287:             CALL DSCAL( N, REC, VI, 1 )
                    288:          END IF
                    289: *
                    290:          IF( RIGHTV ) THEN
                    291: *
                    292: *           LU decomposition with partial pivoting of B, replacing zero
                    293: *           pivots by EPS3.
                    294: *
                    295: *           The imaginary part of the (i,j)-th element of U is stored in
                    296: *           B(j+1,i).
                    297: *
                    298:             B( 2, 1 ) = -WI
                    299:             DO 140 I = 2, N
                    300:                B( I+1, 1 ) = ZERO
                    301:   140       CONTINUE
                    302: *
                    303:             DO 170 I = 1, N - 1
                    304:                ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
                    305:                EI = H( I+1, I )
                    306:                IF( ABSBII.LT.ABS( EI ) ) THEN
                    307: *
                    308: *                 Interchange rows and eliminate.
                    309: *
                    310:                   XR = B( I, I ) / EI
                    311:                   XI = B( I+1, I ) / EI
                    312:                   B( I, I ) = EI
                    313:                   B( I+1, I ) = ZERO
                    314:                   DO 150 J = I + 1, N
                    315:                      TEMP = B( I+1, J )
                    316:                      B( I+1, J ) = B( I, J ) - XR*TEMP
                    317:                      B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
                    318:                      B( I, J ) = TEMP
                    319:                      B( J+1, I ) = ZERO
                    320:   150             CONTINUE
                    321:                   B( I+2, I ) = -WI
                    322:                   B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
                    323:                   B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
                    324:                ELSE
                    325: *
                    326: *                 Eliminate without interchanging rows.
                    327: *
                    328:                   IF( ABSBII.EQ.ZERO ) THEN
                    329:                      B( I, I ) = EPS3
                    330:                      B( I+1, I ) = ZERO
                    331:                      ABSBII = EPS3
                    332:                   END IF
                    333:                   EI = ( EI / ABSBII ) / ABSBII
                    334:                   XR = B( I, I )*EI
                    335:                   XI = -B( I+1, I )*EI
                    336:                   DO 160 J = I + 1, N
                    337:                      B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
                    338:      $                             XI*B( J+1, I )
                    339:                      B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
                    340:   160             CONTINUE
                    341:                   B( I+2, I+1 ) = B( I+2, I+1 ) - WI
                    342:                END IF
                    343: *
                    344: *              Compute 1-norm of offdiagonal elements of i-th row.
                    345: *
                    346:                WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
                    347:      $                     DASUM( N-I, B( I+2, I ), 1 )
                    348:   170       CONTINUE
                    349:             IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
                    350:      $         B( N, N ) = EPS3
                    351:             WORK( N ) = ZERO
                    352: *
                    353:             I1 = N
                    354:             I2 = 1
                    355:             I3 = -1
                    356:          ELSE
                    357: *
                    358: *           UL decomposition with partial pivoting of conjg(B),
                    359: *           replacing zero pivots by EPS3.
                    360: *
                    361: *           The imaginary part of the (i,j)-th element of U is stored in
                    362: *           B(j+1,i).
                    363: *
                    364:             B( N+1, N ) = WI
                    365:             DO 180 J = 1, N - 1
                    366:                B( N+1, J ) = ZERO
                    367:   180       CONTINUE
                    368: *
                    369:             DO 210 J = N, 2, -1
                    370:                EJ = H( J, J-1 )
                    371:                ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
                    372:                IF( ABSBJJ.LT.ABS( EJ ) ) THEN
                    373: *
                    374: *                 Interchange columns and eliminate
                    375: *
                    376:                   XR = B( J, J ) / EJ
                    377:                   XI = B( J+1, J ) / EJ
                    378:                   B( J, J ) = EJ
                    379:                   B( J+1, J ) = ZERO
                    380:                   DO 190 I = 1, J - 1
                    381:                      TEMP = B( I, J-1 )
                    382:                      B( I, J-1 ) = B( I, J ) - XR*TEMP
                    383:                      B( J, I ) = B( J+1, I ) - XI*TEMP
                    384:                      B( I, J ) = TEMP
                    385:                      B( J+1, I ) = ZERO
                    386:   190             CONTINUE
                    387:                   B( J+1, J-1 ) = WI
                    388:                   B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
                    389:                   B( J, J-1 ) = B( J, J-1 ) - XR*WI
                    390:                ELSE
                    391: *
                    392: *                 Eliminate without interchange.
                    393: *
                    394:                   IF( ABSBJJ.EQ.ZERO ) THEN
                    395:                      B( J, J ) = EPS3
                    396:                      B( J+1, J ) = ZERO
                    397:                      ABSBJJ = EPS3
                    398:                   END IF
                    399:                   EJ = ( EJ / ABSBJJ ) / ABSBJJ
                    400:                   XR = B( J, J )*EJ
                    401:                   XI = -B( J+1, J )*EJ
                    402:                   DO 200 I = 1, J - 1
                    403:                      B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
                    404:      $                             XI*B( J+1, I )
                    405:                      B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
                    406:   200             CONTINUE
                    407:                   B( J, J-1 ) = B( J, J-1 ) + WI
                    408:                END IF
                    409: *
                    410: *              Compute 1-norm of offdiagonal elements of j-th column.
                    411: *
                    412:                WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
                    413:      $                     DASUM( J-1, B( J+1, 1 ), LDB )
                    414:   210       CONTINUE
                    415:             IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
                    416:      $         B( 1, 1 ) = EPS3
                    417:             WORK( 1 ) = ZERO
                    418: *
                    419:             I1 = 1
                    420:             I2 = N
                    421:             I3 = 1
                    422:          END IF
                    423: *
                    424:          DO 270 ITS = 1, N
                    425:             SCALE = ONE
                    426:             VMAX = ONE
                    427:             VCRIT = BIGNUM
                    428: *
                    429: *           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
                    430: *             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
                    431: *           overwriting (xr,xi) on (vr,vi).
                    432: *
                    433:             DO 250 I = I1, I2, I3
                    434: *
                    435:                IF( WORK( I ).GT.VCRIT ) THEN
                    436:                   REC = ONE / VMAX
                    437:                   CALL DSCAL( N, REC, VR, 1 )
                    438:                   CALL DSCAL( N, REC, VI, 1 )
                    439:                   SCALE = SCALE*REC
                    440:                   VMAX = ONE
                    441:                   VCRIT = BIGNUM
                    442:                END IF
                    443: *
                    444:                XR = VR( I )
                    445:                XI = VI( I )
                    446:                IF( RIGHTV ) THEN
                    447:                   DO 220 J = I + 1, N
                    448:                      XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
                    449:                      XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
                    450:   220             CONTINUE
                    451:                ELSE
                    452:                   DO 230 J = 1, I - 1
                    453:                      XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
                    454:                      XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
                    455:   230             CONTINUE
                    456:                END IF
                    457: *
                    458:                W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
                    459:                IF( W.GT.SMLNUM ) THEN
                    460:                   IF( W.LT.ONE ) THEN
                    461:                      W1 = ABS( XR ) + ABS( XI )
                    462:                      IF( W1.GT.W*BIGNUM ) THEN
                    463:                         REC = ONE / W1
                    464:                         CALL DSCAL( N, REC, VR, 1 )
                    465:                         CALL DSCAL( N, REC, VI, 1 )
                    466:                         XR = VR( I )
                    467:                         XI = VI( I )
                    468:                         SCALE = SCALE*REC
                    469:                         VMAX = VMAX*REC
                    470:                      END IF
                    471:                   END IF
                    472: *
                    473: *                 Divide by diagonal element of B.
                    474: *
                    475:                   CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
                    476:      $                         VI( I ) )
                    477:                   VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
                    478:                   VCRIT = BIGNUM / VMAX
                    479:                ELSE
                    480:                   DO 240 J = 1, N
                    481:                      VR( J ) = ZERO
                    482:                      VI( J ) = ZERO
                    483:   240             CONTINUE
                    484:                   VR( I ) = ONE
                    485:                   VI( I ) = ONE
                    486:                   SCALE = ZERO
                    487:                   VMAX = ONE
                    488:                   VCRIT = BIGNUM
                    489:                END IF
                    490:   250       CONTINUE
                    491: *
                    492: *           Test for sufficient growth in the norm of (VR,VI).
                    493: *
                    494:             VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
                    495:             IF( VNORM.GE.GROWTO*SCALE )
                    496:      $         GO TO 280
                    497: *
                    498: *           Choose a new orthogonal starting vector and try again.
                    499: *
                    500:             Y = EPS3 / ( ROOTN+ONE )
                    501:             VR( 1 ) = EPS3
                    502:             VI( 1 ) = ZERO
                    503: *
                    504:             DO 260 I = 2, N
                    505:                VR( I ) = Y
                    506:                VI( I ) = ZERO
                    507:   260       CONTINUE
                    508:             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
                    509:   270    CONTINUE
                    510: *
                    511: *        Failure to find eigenvector in N iterations
                    512: *
                    513:          INFO = 1
                    514: *
                    515:   280    CONTINUE
                    516: *
                    517: *        Normalize eigenvector.
                    518: *
                    519:          VNORM = ZERO
                    520:          DO 290 I = 1, N
                    521:             VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
                    522:   290    CONTINUE
                    523:          CALL DSCAL( N, ONE / VNORM, VR, 1 )
                    524:          CALL DSCAL( N, ONE / VNORM, VI, 1 )
                    525: *
                    526:       END IF
                    527: *
                    528:       RETURN
                    529: *
                    530: *     End of DLAEIN
                    531: *
                    532:       END

CVSweb interface <joel.bertrand@systella.fr>