|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 |
SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
$ EPS3, SMLNUM, INFO ) * * -- LAPACK auxiliary routine (version 3.3.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- April 2011 -- * * .. Scalar Arguments .. LOGICAL NOINIT, RIGHTV INTEGER INFO, LDB, LDH, N DOUBLE PRECISION EPS3, SMLNUM COMPLEX*16 W * .. * .. Array Arguments .. DOUBLE PRECISION RWORK( * ) COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * ) * .. * * Purpose * ======= * * ZLAEIN uses inverse iteration to find a right or left eigenvector * corresponding to the eigenvalue W of a complex upper Hessenberg * matrix H. * * Arguments * ========= * * RIGHTV (input) LOGICAL * = .TRUE. : compute right eigenvector; * = .FALSE.: compute left eigenvector. * * NOINIT (input) LOGICAL * = .TRUE. : no initial vector supplied in V * = .FALSE.: initial vector supplied in V. * * N (input) INTEGER * The order of the matrix H. N >= 0. * * H (input) COMPLEX*16 array, dimension (LDH,N) * The upper Hessenberg matrix H. * * LDH (input) INTEGER * The leading dimension of the array H. LDH >= max(1,N). * * W (input) COMPLEX*16 * The eigenvalue of H whose corresponding right or left * eigenvector is to be computed. * * V (input/output) COMPLEX*16 array, dimension (N) * On entry, if NOINIT = .FALSE., V must contain a starting * vector for inverse iteration; otherwise V need not be set. * On exit, V contains the computed eigenvector, normalized so * that the component of largest magnitude has magnitude 1; here * the magnitude of a complex number (x,y) is taken to be * |x| + |y|. * * B (workspace) COMPLEX*16 array, dimension (LDB,N) * * LDB (input) INTEGER * The leading dimension of the array B. LDB >= max(1,N). * * RWORK (workspace) DOUBLE PRECISION array, dimension (N) * * EPS3 (input) DOUBLE PRECISION * A small machine-dependent value which is used to perturb * close eigenvalues, and to replace zero pivots. * * SMLNUM (input) DOUBLE PRECISION * A machine-dependent value close to the underflow threshold. * * INFO (output) INTEGER * = 0: successful exit * = 1: inverse iteration did not converge; V is set to the * last iterate. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, TENTH PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Local Scalars .. CHARACTER NORMIN, TRANS INTEGER I, IERR, ITS, J DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM COMPLEX*16 CDUM, EI, EJ, TEMP, X * .. * .. External Functions .. INTEGER IZAMAX DOUBLE PRECISION DZASUM, DZNRM2 COMPLEX*16 ZLADIV EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZLATRS * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT * .. * .. Statement Functions .. DOUBLE PRECISION CABS1 * .. * .. Statement Function definitions .. CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) * .. * .. Executable Statements .. * INFO = 0 * * GROWTO is the threshold used in the acceptance test for an * eigenvector. * ROOTN = SQRT( DBLE( N ) ) GROWTO = TENTH / ROOTN NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM * * Form B = H - W*I (except that the subdiagonal elements are not * stored). * DO 20 J = 1, N DO 10 I = 1, J - 1 B( I, J ) = H( I, J ) 10 CONTINUE B( J, J ) = H( J, J ) - W 20 CONTINUE * IF( NOINIT ) THEN * * Initialize V. * DO 30 I = 1, N V( I ) = EPS3 30 CONTINUE ELSE * * Scale supplied initial vector. * VNORM = DZNRM2( N, V, 1 ) CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 ) END IF * IF( RIGHTV ) THEN * * LU decomposition with partial pivoting of B, replacing zero * pivots by EPS3. * DO 60 I = 1, N - 1 EI = H( I+1, I ) IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN * * Interchange rows and eliminate. * X = ZLADIV( B( I, I ), EI ) B( I, I ) = EI DO 40 J = I + 1, N TEMP = B( I+1, J ) B( I+1, J ) = B( I, J ) - X*TEMP B( I, J ) = TEMP 40 CONTINUE ELSE * * Eliminate without interchange. * IF( B( I, I ).EQ.ZERO ) $ B( I, I ) = EPS3 X = ZLADIV( EI, B( I, I ) ) IF( X.NE.ZERO ) THEN DO 50 J = I + 1, N B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 50 CONTINUE END IF END IF 60 CONTINUE IF( B( N, N ).EQ.ZERO ) $ B( N, N ) = EPS3 * TRANS = 'N' * ELSE * * UL decomposition with partial pivoting of B, replacing zero * pivots by EPS3. * DO 90 J = N, 2, -1 EJ = H( J, J-1 ) IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN * * Interchange columns and eliminate. * X = ZLADIV( B( J, J ), EJ ) B( J, J ) = EJ DO 70 I = 1, J - 1 TEMP = B( I, J-1 ) B( I, J-1 ) = B( I, J ) - X*TEMP B( I, J ) = TEMP 70 CONTINUE ELSE * * Eliminate without interchange. * IF( B( J, J ).EQ.ZERO ) $ B( J, J ) = EPS3 X = ZLADIV( EJ, B( J, J ) ) IF( X.NE.ZERO ) THEN DO 80 I = 1, J - 1 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 80 CONTINUE END IF END IF 90 CONTINUE IF( B( 1, 1 ).EQ.ZERO ) $ B( 1, 1 ) = EPS3 * TRANS = 'C' * END IF * NORMIN = 'N' DO 110 ITS = 1, N * * Solve U*x = scale*v for a right eigenvector * or U**H *x = scale*v for a left eigenvector, * overwriting x on v. * CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V, $ SCALE, RWORK, IERR ) NORMIN = 'Y' * * Test for sufficient growth in the norm of v. * VNORM = DZASUM( N, V, 1 ) IF( VNORM.GE.GROWTO*SCALE ) $ GO TO 120 * * Choose new orthogonal starting vector and try again. * RTEMP = EPS3 / ( ROOTN+ONE ) V( 1 ) = EPS3 DO 100 I = 2, N V( I ) = RTEMP 100 CONTINUE V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN 110 CONTINUE * * Failure to find eigenvector in N iterations. * INFO = 1 * 120 CONTINUE * * Normalize eigenvector. * I = IZAMAX( N, V, 1 ) CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 ) * RETURN * * End of ZLAEIN * END |