LAPACK: SRC/claein.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
00002      $                   EPS3, SMLNUM, INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.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 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       LOGICAL            NOINIT, RIGHTV
00011       INTEGER            INFO, LDB, LDH, N
00012       REAL               EPS3, SMLNUM
00013       COMPLEX            W
00014 *     ..
00015 *     .. Array Arguments ..
00016       REAL               RWORK( * )
00017       COMPLEX            B( LDB, * ), H( LDH, * ), V( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CLAEIN uses inverse iteration to find a right or left eigenvector
00024 *  corresponding to the eigenvalue W of a complex upper Hessenberg
00025 *  matrix H.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  RIGHTV   (input) LOGICAL
00031 *          = .TRUE. : compute right eigenvector;
00032 *          = .FALSE.: compute left eigenvector.
00033 *
00034 *  NOINIT   (input) LOGICAL
00035 *          = .TRUE. : no initial vector supplied in V
00036 *          = .FALSE.: initial vector supplied in V.
00037 *
00038 *  N       (input) INTEGER
00039 *          The order of the matrix H.  N >= 0.
00040 *
00041 *  H       (input) COMPLEX array, dimension (LDH,N)
00042 *          The upper Hessenberg matrix H.
00043 *
00044 *  LDH     (input) INTEGER
00045 *          The leading dimension of the array H.  LDH >= max(1,N).
00046 *
00047 *  W       (input) COMPLEX
00048 *          The eigenvalue of H whose corresponding right or left
00049 *          eigenvector is to be computed.
00050 *
00051 *  V       (input/output) COMPLEX array, dimension (N)
00052 *          On entry, if NOINIT = .FALSE., V must contain a starting
00053 *          vector for inverse iteration; otherwise V need not be set.
00054 *          On exit, V contains the computed eigenvector, normalized so
00055 *          that the component of largest magnitude has magnitude 1; here
00056 *          the magnitude of a complex number (x,y) is taken to be
00057 *          |x| + |y|.
00058 *
00059 *  B       (workspace) COMPLEX array, dimension (LDB,N)
00060 *
00061 *  LDB     (input) INTEGER
00062 *          The leading dimension of the array B.  LDB >= max(1,N).
00063 *
00064 *  RWORK   (workspace) REAL array, dimension (N)
00065 *
00066 *  EPS3    (input) REAL
00067 *          A small machine-dependent value which is used to perturb
00068 *          close eigenvalues, and to replace zero pivots.
00069 *
00070 *  SMLNUM  (input) REAL
00071 *          A machine-dependent value close to the underflow threshold.
00072 *
00073 *  INFO    (output) INTEGER
00074 *          = 0:  successful exit
00075 *          = 1:  inverse iteration did not converge; V is set to the
00076 *                last iterate.
00077 *
00078 *  =====================================================================
00079 *
00080 *     .. Parameters ..
00081       REAL               ONE, TENTH
00082       PARAMETER          ( ONE = 1.0E+0, TENTH = 1.0E-1 )
00083       COMPLEX            ZERO
00084       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ) )
00085 *     ..
00086 *     .. Local Scalars ..
00087       CHARACTER          NORMIN, TRANS
00088       INTEGER            I, IERR, ITS, J
00089       REAL               GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
00090       COMPLEX            CDUM, EI, EJ, TEMP, X
00091 *     ..
00092 *     .. External Functions ..
00093       INTEGER            ICAMAX
00094       REAL               SCASUM, SCNRM2
00095       COMPLEX            CLADIV
00096       EXTERNAL           ICAMAX, SCASUM, SCNRM2, CLADIV
00097 *     ..
00098 *     .. External Subroutines ..
00099       EXTERNAL           CLATRS, CSSCAL
00100 *     ..
00101 *     .. Intrinsic Functions ..
00102       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
00103 *     ..
00104 *     .. Statement Functions ..
00105       REAL               CABS1
00106 *     ..
00107 *     .. Statement Function definitions ..
00108       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00109 *     ..
00110 *     .. Executable Statements ..
00111 *
00112       INFO = 0
00113 *
00114 *     GROWTO is the threshold used in the acceptance test for an
00115 *     eigenvector.
00116 *
00117       ROOTN = SQRT( REAL( N ) )
00118       GROWTO = TENTH / ROOTN
00119       NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
00120 *
00121 *     Form B = H - W*I (except that the subdiagonal elements are not
00122 *     stored).
00123 *
00124       DO 20 J = 1, N
00125          DO 10 I = 1, J - 1
00126             B( I, J ) = H( I, J )
00127    10    CONTINUE
00128          B( J, J ) = H( J, J ) - W
00129    20 CONTINUE
00130 *
00131       IF( NOINIT ) THEN
00132 *
00133 *        Initialize V.
00134 *
00135          DO 30 I = 1, N
00136             V( I ) = EPS3
00137    30    CONTINUE
00138       ELSE
00139 *
00140 *        Scale supplied initial vector.
00141 *
00142          VNORM = SCNRM2( N, V, 1 )
00143          CALL CSSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
00144       END IF
00145 *
00146       IF( RIGHTV ) THEN
00147 *
00148 *        LU decomposition with partial pivoting of B, replacing zero
00149 *        pivots by EPS3.
00150 *
00151          DO 60 I = 1, N - 1
00152             EI = H( I+1, I )
00153             IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
00154 *
00155 *              Interchange rows and eliminate.
00156 *
00157                X = CLADIV( B( I, I ), EI )
00158                B( I, I ) = EI
00159                DO 40 J = I + 1, N
00160                   TEMP = B( I+1, J )
00161                   B( I+1, J ) = B( I, J ) - X*TEMP
00162                   B( I, J ) = TEMP
00163    40          CONTINUE
00164             ELSE
00165 *
00166 *              Eliminate without interchange.
00167 *
00168                IF( B( I, I ).EQ.ZERO )
00169      $            B( I, I ) = EPS3
00170                X = CLADIV( EI, B( I, I ) )
00171                IF( X.NE.ZERO ) THEN
00172                   DO 50 J = I + 1, N
00173                      B( I+1, J ) = B( I+1, J ) - X*B( I, J )
00174    50             CONTINUE
00175                END IF
00176             END IF
00177    60    CONTINUE
00178          IF( B( N, N ).EQ.ZERO )
00179      $      B( N, N ) = EPS3
00180 *
00181          TRANS = 'N'
00182 *
00183       ELSE
00184 *
00185 *        UL decomposition with partial pivoting of B, replacing zero
00186 *        pivots by EPS3.
00187 *
00188          DO 90 J = N, 2, -1
00189             EJ = H( J, J-1 )
00190             IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
00191 *
00192 *              Interchange columns and eliminate.
00193 *
00194                X = CLADIV( B( J, J ), EJ )
00195                B( J, J ) = EJ
00196                DO 70 I = 1, J - 1
00197                   TEMP = B( I, J-1 )
00198                   B( I, J-1 ) = B( I, J ) - X*TEMP
00199                   B( I, J ) = TEMP
00200    70          CONTINUE
00201             ELSE
00202 *
00203 *              Eliminate without interchange.
00204 *
00205                IF( B( J, J ).EQ.ZERO )
00206      $            B( J, J ) = EPS3
00207                X = CLADIV( EJ, B( J, J ) )
00208                IF( X.NE.ZERO ) THEN
00209                   DO 80 I = 1, J - 1
00210                      B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
00211    80             CONTINUE
00212                END IF
00213             END IF
00214    90    CONTINUE
00215          IF( B( 1, 1 ).EQ.ZERO )
00216      $      B( 1, 1 ) = EPS3
00217 *
00218          TRANS = 'C'
00219 *
00220       END IF
00221 *
00222       NORMIN = 'N'
00223       DO 110 ITS = 1, N
00224 *
00225 *        Solve U*x = scale*v for a right eigenvector
00226 *          or U'*x = scale*v for a left eigenvector,
00227 *        overwriting x on v.
00228 *
00229          CALL CLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
00230      $                SCALE, RWORK, IERR )
00231          NORMIN = 'Y'
00232 *
00233 *        Test for sufficient growth in the norm of v.
00234 *
00235          VNORM = SCASUM( N, V, 1 )
00236          IF( VNORM.GE.GROWTO*SCALE )
00237      $      GO TO 120
00238 *
00239 *        Choose new orthogonal starting vector and try again.
00240 *
00241          RTEMP = EPS3 / ( ROOTN+ONE )
00242          V( 1 ) = EPS3
00243          DO 100 I = 2, N
00244             V( I ) = RTEMP
00245   100    CONTINUE
00246          V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
00247   110 CONTINUE
00248 *
00249 *     Failure to find eigenvector in N iterations.
00250 *
00251       INFO = 1
00252 *
00253   120 CONTINUE
00254 *
00255 *     Normalize eigenvector.
00256 *
00257       I = ICAMAX( N, V, 1 )
00258       CALL CSSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
00259 *
00260       RETURN
00261 *
00262 *     End of CLAEIN
00263 *
00264       END