LAPACK: SRC/dlaein.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
00002      $                   LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            NOINIT, RIGHTV
00011       INTEGER            INFO, LDB, LDH, N
00012       DOUBLE PRECISION   BIGNUM, EPS3, SMLNUM, WI, WR
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
00016      $                   WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DLAEIN uses inverse iteration to find a right or left eigenvector
00023 *  corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
00024 *  matrix H.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  RIGHTV  (input) LOGICAL
00030 *          = .TRUE. : compute right eigenvector;
00031 *          = .FALSE.: compute left eigenvector.
00032 *
00033 *  NOINIT  (input) LOGICAL
00034 *          = .TRUE. : no initial vector supplied in (VR,VI).
00035 *          = .FALSE.: initial vector supplied in (VR,VI).
00036 *
00037 *  N       (input) INTEGER
00038 *          The order of the matrix H.  N >= 0.
00039 *
00040 *  H       (input) DOUBLE PRECISION array, dimension (LDH,N)
00041 *          The upper Hessenberg matrix H.
00042 *
00043 *  LDH     (input) INTEGER
00044 *          The leading dimension of the array H.  LDH >= max(1,N).
00045 *
00046 *  WR      (input) DOUBLE PRECISION
00047 *  WI      (input) DOUBLE PRECISION
00048 *          The real and imaginary parts of the eigenvalue of H whose
00049 *          corresponding right or left eigenvector is to be computed.
00050 *
00051 *  VR      (input/output) DOUBLE PRECISION array, dimension (N)
00052 *  VI      (input/output) DOUBLE PRECISION array, dimension (N)
00053 *          On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
00054 *          a real starting vector for inverse iteration using the real
00055 *          eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
00056 *          must contain the real and imaginary parts of a complex
00057 *          starting vector for inverse iteration using the complex
00058 *          eigenvalue (WR,WI); otherwise VR and VI need not be set.
00059 *          On exit, if WI = 0.0 (real eigenvalue), VR contains the
00060 *          computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
00061 *          VR and VI contain the real and imaginary parts of the
00062 *          computed complex eigenvector. The eigenvector is normalized
00063 *          so that the component of largest magnitude has magnitude 1;
00064 *          here the magnitude of a complex number (x,y) is taken to be
00065 *          |x| + |y|.
00066 *          VI is not referenced if WI = 0.0.
00067 *
00068 *  B       (workspace) DOUBLE PRECISION array, dimension (LDB,N)
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of the array B.  LDB >= N+1.
00072 *
00073 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00074 *
00075 *  EPS3    (input) DOUBLE PRECISION
00076 *          A small machine-dependent value which is used to perturb
00077 *          close eigenvalues, and to replace zero pivots.
00078 *
00079 *  SMLNUM  (input) DOUBLE PRECISION
00080 *          A machine-dependent value close to the underflow threshold.
00081 *
00082 *  BIGNUM  (input) DOUBLE PRECISION
00083 *          A machine-dependent value close to the overflow threshold.
00084 *
00085 *  INFO    (output) INTEGER
00086 *          = 0:  successful exit
00087 *          = 1:  inverse iteration did not converge; VR is set to the
00088 *                last iterate, and so is VI if WI.ne.0.0.
00089 *
00090 *  =====================================================================
00091 *
00092 *     .. Parameters ..
00093       DOUBLE PRECISION   ZERO, ONE, TENTH
00094       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       CHARACTER          NORMIN, TRANS
00098       INTEGER            I, I1, I2, I3, IERR, ITS, J
00099       DOUBLE PRECISION   ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
00100      $                   REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
00101      $                   W1, X, XI, XR, Y
00102 *     ..
00103 *     .. External Functions ..
00104       INTEGER            IDAMAX
00105       DOUBLE PRECISION   DASUM, DLAPY2, DNRM2
00106       EXTERNAL           IDAMAX, DASUM, DLAPY2, DNRM2
00107 *     ..
00108 *     .. External Subroutines ..
00109       EXTERNAL           DLADIV, DLATRS, DSCAL
00110 *     ..
00111 *     .. Intrinsic Functions ..
00112       INTRINSIC          ABS, DBLE, MAX, SQRT
00113 *     ..
00114 *     .. Executable Statements ..
00115 *
00116       INFO = 0
00117 *
00118 *     GROWTO is the threshold used in the acceptance test for an
00119 *     eigenvector.
00120 *
00121       ROOTN = SQRT( DBLE( N ) )
00122       GROWTO = TENTH / ROOTN
00123       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
00124 *
00125 *     Form B = H - (WR,WI)*I (except that the subdiagonal elements and
00126 *     the imaginary parts of the diagonal elements are not stored).
00127 *
00128       DO 20 J = 1, N
00129          DO 10 I = 1, J - 1
00130             B( I, J ) = H( I, J )
00131    10    CONTINUE
00132          B( J, J ) = H( J, J ) - WR
00133    20 CONTINUE
00134 *
00135       IF( WI.EQ.ZERO ) THEN
00136 *
00137 *        Real eigenvalue.
00138 *
00139          IF( NOINIT ) THEN
00140 *
00141 *           Set initial vector.
00142 *
00143             DO 30 I = 1, N
00144                VR( I ) = EPS3
00145    30       CONTINUE
00146          ELSE
00147 *
00148 *           Scale supplied initial vector.
00149 *
00150             VNORM = DNRM2( N, VR, 1 )
00151             CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
00152      $                  1 )
00153          END IF
00154 *
00155          IF( RIGHTV ) THEN
00156 *
00157 *           LU decomposition with partial pivoting of B, replacing zero
00158 *           pivots by EPS3.
00159 *
00160             DO 60 I = 1, N - 1
00161                EI = H( I+1, I )
00162                IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
00163 *
00164 *                 Interchange rows and eliminate.
00165 *
00166                   X = B( I, I ) / EI
00167                   B( I, I ) = EI
00168                   DO 40 J = I + 1, N
00169                      TEMP = B( I+1, J )
00170                      B( I+1, J ) = B( I, J ) - X*TEMP
00171                      B( I, J ) = TEMP
00172    40             CONTINUE
00173                ELSE
00174 *
00175 *                 Eliminate without interchange.
00176 *
00177                   IF( B( I, I ).EQ.ZERO )
00178      $               B( I, I ) = EPS3
00179                   X = EI / B( I, I )
00180                   IF( X.NE.ZERO ) THEN
00181                      DO 50 J = I + 1, N
00182                         B( I+1, J ) = B( I+1, J ) - X*B( I, J )
00183    50                CONTINUE
00184                   END IF
00185                END IF
00186    60       CONTINUE
00187             IF( B( N, N ).EQ.ZERO )
00188      $         B( N, N ) = EPS3
00189 *
00190             TRANS = 'N'
00191 *
00192          ELSE
00193 *
00194 *           UL decomposition with partial pivoting of B, replacing zero
00195 *           pivots by EPS3.
00196 *
00197             DO 90 J = N, 2, -1
00198                EJ = H( J, J-1 )
00199                IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
00200 *
00201 *                 Interchange columns and eliminate.
00202 *
00203                   X = B( J, J ) / EJ
00204                   B( J, J ) = EJ
00205                   DO 70 I = 1, J - 1
00206                      TEMP = B( I, J-1 )
00207                      B( I, J-1 ) = B( I, J ) - X*TEMP
00208                      B( I, J ) = TEMP
00209    70             CONTINUE
00210                ELSE
00211 *
00212 *                 Eliminate without interchange.
00213 *
00214                   IF( B( J, J ).EQ.ZERO )
00215      $               B( J, J ) = EPS3
00216                   X = EJ / B( J, J )
00217                   IF( X.NE.ZERO ) THEN
00218                      DO 80 I = 1, J - 1
00219                         B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
00220    80                CONTINUE
00221                   END IF
00222                END IF
00223    90       CONTINUE
00224             IF( B( 1, 1 ).EQ.ZERO )
00225      $         B( 1, 1 ) = EPS3
00226 *
00227             TRANS = 'T'
00228 *
00229          END IF
00230 *
00231          NORMIN = 'N'
00232          DO 110 ITS = 1, N
00233 *
00234 *           Solve U*x = scale*v for a right eigenvector
00235 *             or U'*x = scale*v for a left eigenvector,
00236 *           overwriting x on v.
00237 *
00238             CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
00239      $                   VR, SCALE, WORK, IERR )
00240             NORMIN = 'Y'
00241 *
00242 *           Test for sufficient growth in the norm of v.
00243 *
00244             VNORM = DASUM( N, VR, 1 )
00245             IF( VNORM.GE.GROWTO*SCALE )
00246      $         GO TO 120
00247 *
00248 *           Choose new orthogonal starting vector and try again.
00249 *
00250             TEMP = EPS3 / ( ROOTN+ONE )
00251             VR( 1 ) = EPS3
00252             DO 100 I = 2, N
00253                VR( I ) = TEMP
00254   100       CONTINUE
00255             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
00256   110    CONTINUE
00257 *
00258 *        Failure to find eigenvector in N iterations.
00259 *
00260          INFO = 1
00261 *
00262   120    CONTINUE
00263 *
00264 *        Normalize eigenvector.
00265 *
00266          I = IDAMAX( N, VR, 1 )
00267          CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
00268       ELSE
00269 *
00270 *        Complex eigenvalue.
00271 *
00272          IF( NOINIT ) THEN
00273 *
00274 *           Set initial vector.
00275 *
00276             DO 130 I = 1, N
00277                VR( I ) = EPS3
00278                VI( I ) = ZERO
00279   130       CONTINUE
00280          ELSE
00281 *
00282 *           Scale supplied initial vector.
00283 *
00284             NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
00285             REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
00286             CALL DSCAL( N, REC, VR, 1 )
00287             CALL DSCAL( N, REC, VI, 1 )
00288          END IF
00289 *
00290          IF( RIGHTV ) THEN
00291 *
00292 *           LU decomposition with partial pivoting of B, replacing zero
00293 *           pivots by EPS3.
00294 *
00295 *           The imaginary part of the (i,j)-th element of U is stored in
00296 *           B(j+1,i).
00297 *
00298             B( 2, 1 ) = -WI
00299             DO 140 I = 2, N
00300                B( I+1, 1 ) = ZERO
00301   140       CONTINUE
00302 *
00303             DO 170 I = 1, N - 1
00304                ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
00305                EI = H( I+1, I )
00306                IF( ABSBII.LT.ABS( EI ) ) THEN
00307 *
00308 *                 Interchange rows and eliminate.
00309 *
00310                   XR = B( I, I ) / EI
00311                   XI = B( I+1, I ) / EI
00312                   B( I, I ) = EI
00313                   B( I+1, I ) = ZERO
00314                   DO 150 J = I + 1, N
00315                      TEMP = B( I+1, J )
00316                      B( I+1, J ) = B( I, J ) - XR*TEMP
00317                      B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
00318                      B( I, J ) = TEMP
00319                      B( J+1, I ) = ZERO
00320   150             CONTINUE
00321                   B( I+2, I ) = -WI
00322                   B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
00323                   B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
00324                ELSE
00325 *
00326 *                 Eliminate without interchanging rows.
00327 *
00328                   IF( ABSBII.EQ.ZERO ) THEN
00329                      B( I, I ) = EPS3
00330                      B( I+1, I ) = ZERO
00331                      ABSBII = EPS3
00332                   END IF
00333                   EI = ( EI / ABSBII ) / ABSBII
00334                   XR = B( I, I )*EI
00335                   XI = -B( I+1, I )*EI
00336                   DO 160 J = I + 1, N
00337                      B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
00338      $                             XI*B( J+1, I )
00339                      B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
00340   160             CONTINUE
00341                   B( I+2, I+1 ) = B( I+2, I+1 ) - WI
00342                END IF
00343 *
00344 *              Compute 1-norm of offdiagonal elements of i-th row.
00345 *
00346                WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
00347      $                     DASUM( N-I, B( I+2, I ), 1 )
00348   170       CONTINUE
00349             IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
00350      $         B( N, N ) = EPS3
00351             WORK( N ) = ZERO
00352 *
00353             I1 = N
00354             I2 = 1
00355             I3 = -1
00356          ELSE
00357 *
00358 *           UL decomposition with partial pivoting of conjg(B),
00359 *           replacing zero pivots by EPS3.
00360 *
00361 *           The imaginary part of the (i,j)-th element of U is stored in
00362 *           B(j+1,i).
00363 *
00364             B( N+1, N ) = WI
00365             DO 180 J = 1, N - 1
00366                B( N+1, J ) = ZERO
00367   180       CONTINUE
00368 *
00369             DO 210 J = N, 2, -1
00370                EJ = H( J, J-1 )
00371                ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
00372                IF( ABSBJJ.LT.ABS( EJ ) ) THEN
00373 *
00374 *                 Interchange columns and eliminate
00375 *
00376                   XR = B( J, J ) / EJ
00377                   XI = B( J+1, J ) / EJ
00378                   B( J, J ) = EJ
00379                   B( J+1, J ) = ZERO
00380                   DO 190 I = 1, J - 1
00381                      TEMP = B( I, J-1 )
00382                      B( I, J-1 ) = B( I, J ) - XR*TEMP
00383                      B( J, I ) = B( J+1, I ) - XI*TEMP
00384                      B( I, J ) = TEMP
00385                      B( J+1, I ) = ZERO
00386   190             CONTINUE
00387                   B( J+1, J-1 ) = WI
00388                   B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
00389                   B( J, J-1 ) = B( J, J-1 ) - XR*WI
00390                ELSE
00391 *
00392 *                 Eliminate without interchange.
00393 *
00394                   IF( ABSBJJ.EQ.ZERO ) THEN
00395                      B( J, J ) = EPS3
00396                      B( J+1, J ) = ZERO
00397                      ABSBJJ = EPS3
00398                   END IF
00399                   EJ = ( EJ / ABSBJJ ) / ABSBJJ
00400                   XR = B( J, J )*EJ
00401                   XI = -B( J+1, J )*EJ
00402                   DO 200 I = 1, J - 1
00403                      B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
00404      $                             XI*B( J+1, I )
00405                      B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
00406   200             CONTINUE
00407                   B( J, J-1 ) = B( J, J-1 ) + WI
00408                END IF
00409 *
00410 *              Compute 1-norm of offdiagonal elements of j-th column.
00411 *
00412                WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
00413      $                     DASUM( J-1, B( J+1, 1 ), LDB )
00414   210       CONTINUE
00415             IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
00416      $         B( 1, 1 ) = EPS3
00417             WORK( 1 ) = ZERO
00418 *
00419             I1 = 1
00420             I2 = N
00421             I3 = 1
00422          END IF
00423 *
00424          DO 270 ITS = 1, N
00425             SCALE = ONE
00426             VMAX = ONE
00427             VCRIT = BIGNUM
00428 *
00429 *           Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
00430 *             or U'*(xr,xi) = scale*(vr,vi) for a left eigenvector,
00431 *           overwriting (xr,xi) on (vr,vi).
00432 *
00433             DO 250 I = I1, I2, I3
00434 *
00435                IF( WORK( I ).GT.VCRIT ) THEN
00436                   REC = ONE / VMAX
00437                   CALL DSCAL( N, REC, VR, 1 )
00438                   CALL DSCAL( N, REC, VI, 1 )
00439                   SCALE = SCALE*REC
00440                   VMAX = ONE
00441                   VCRIT = BIGNUM
00442                END IF
00443 *
00444                XR = VR( I )
00445                XI = VI( I )
00446                IF( RIGHTV ) THEN
00447                   DO 220 J = I + 1, N
00448                      XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
00449                      XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
00450   220             CONTINUE
00451                ELSE
00452                   DO 230 J = 1, I - 1
00453                      XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
00454                      XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
00455   230             CONTINUE
00456                END IF
00457 *
00458                W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
00459                IF( W.GT.SMLNUM ) THEN
00460                   IF( W.LT.ONE ) THEN
00461                      W1 = ABS( XR ) + ABS( XI )
00462                      IF( W1.GT.W*BIGNUM ) THEN
00463                         REC = ONE / W1
00464                         CALL DSCAL( N, REC, VR, 1 )
00465                         CALL DSCAL( N, REC, VI, 1 )
00466                         XR = VR( I )
00467                         XI = VI( I )
00468                         SCALE = SCALE*REC
00469                         VMAX = VMAX*REC
00470                      END IF
00471                   END IF
00472 *
00473 *                 Divide by diagonal element of B.
00474 *
00475                   CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
00476      $                         VI( I ) )
00477                   VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
00478                   VCRIT = BIGNUM / VMAX
00479                ELSE
00480                   DO 240 J = 1, N
00481                      VR( J ) = ZERO
00482                      VI( J ) = ZERO
00483   240             CONTINUE
00484                   VR( I ) = ONE
00485                   VI( I ) = ONE
00486                   SCALE = ZERO
00487                   VMAX = ONE
00488                   VCRIT = BIGNUM
00489                END IF
00490   250       CONTINUE
00491 *
00492 *           Test for sufficient growth in the norm of (VR,VI).
00493 *
00494             VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
00495             IF( VNORM.GE.GROWTO*SCALE )
00496      $         GO TO 280
00497 *
00498 *           Choose a new orthogonal starting vector and try again.
00499 *
00500             Y = EPS3 / ( ROOTN+ONE )
00501             VR( 1 ) = EPS3
00502             VI( 1 ) = ZERO
00503 *
00504             DO 260 I = 2, N
00505                VR( I ) = Y
00506                VI( I ) = ZERO
00507   260       CONTINUE
00508             VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
00509   270    CONTINUE
00510 *
00511 *        Failure to find eigenvector in N iterations
00512 *
00513          INFO = 1
00514 *
00515   280    CONTINUE
00516 *
00517 *        Normalize eigenvector.
00518 *
00519          VNORM = ZERO
00520          DO 290 I = 1, N
00521             VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
00522   290    CONTINUE
00523          CALL DSCAL( N, ONE / VNORM, VR, 1 )
00524          CALL DSCAL( N, ONE / VNORM, VI, 1 )
00525 *
00526       END IF
00527 *
00528       RETURN
00529 *
00530 *     End of DLAEIN
00531 *
00532       END