![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLAEIN 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLAEIN + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, 00022 * EPS3, SMLNUM, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL NOINIT, RIGHTV 00026 * INTEGER INFO, LDB, LDH, N 00027 * DOUBLE PRECISION EPS3, SMLNUM 00028 * COMPLEX*16 W 00029 * .. 00030 * .. Array Arguments .. 00031 * DOUBLE PRECISION RWORK( * ) 00032 * COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * ) 00033 * .. 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> ZLAEIN uses inverse iteration to find a right or left eigenvector 00042 *> corresponding to the eigenvalue W of a complex upper Hessenberg 00043 *> matrix H. 00044 *> \endverbatim 00045 * 00046 * Arguments: 00047 * ========== 00048 * 00049 *> \param[in] RIGHTV 00050 *> \verbatim 00051 *> RIGHTV is LOGICAL 00052 *> = .TRUE. : compute right eigenvector; 00053 *> = .FALSE.: compute left eigenvector. 00054 *> \endverbatim 00055 *> 00056 *> \param[in] NOINIT 00057 *> \verbatim 00058 *> NOINIT is LOGICAL 00059 *> = .TRUE. : no initial vector supplied in V 00060 *> = .FALSE.: initial vector supplied in V. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] N 00064 *> \verbatim 00065 *> N is INTEGER 00066 *> The order of the matrix H. N >= 0. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] H 00070 *> \verbatim 00071 *> H is COMPLEX*16 array, dimension (LDH,N) 00072 *> The upper Hessenberg matrix H. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] LDH 00076 *> \verbatim 00077 *> LDH is INTEGER 00078 *> The leading dimension of the array H. LDH >= max(1,N). 00079 *> \endverbatim 00080 *> 00081 *> \param[in] W 00082 *> \verbatim 00083 *> W is COMPLEX*16 00084 *> The eigenvalue of H whose corresponding right or left 00085 *> eigenvector is to be computed. 00086 *> \endverbatim 00087 *> 00088 *> \param[in,out] V 00089 *> \verbatim 00090 *> V is COMPLEX*16 array, dimension (N) 00091 *> On entry, if NOINIT = .FALSE., V must contain a starting 00092 *> vector for inverse iteration; otherwise V need not be set. 00093 *> On exit, V contains the computed eigenvector, normalized so 00094 *> that the component of largest magnitude has magnitude 1; here 00095 *> the magnitude of a complex number (x,y) is taken to be 00096 *> |x| + |y|. 00097 *> \endverbatim 00098 *> 00099 *> \param[out] B 00100 *> \verbatim 00101 *> B is COMPLEX*16 array, dimension (LDB,N) 00102 *> \endverbatim 00103 *> 00104 *> \param[in] LDB 00105 *> \verbatim 00106 *> LDB is INTEGER 00107 *> The leading dimension of the array B. LDB >= max(1,N). 00108 *> \endverbatim 00109 *> 00110 *> \param[out] RWORK 00111 *> \verbatim 00112 *> RWORK is DOUBLE PRECISION array, dimension (N) 00113 *> \endverbatim 00114 *> 00115 *> \param[in] EPS3 00116 *> \verbatim 00117 *> EPS3 is DOUBLE PRECISION 00118 *> A small machine-dependent value which is used to perturb 00119 *> close eigenvalues, and to replace zero pivots. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] SMLNUM 00123 *> \verbatim 00124 *> SMLNUM is DOUBLE PRECISION 00125 *> A machine-dependent value close to the underflow threshold. 00126 *> \endverbatim 00127 *> 00128 *> \param[out] INFO 00129 *> \verbatim 00130 *> INFO is INTEGER 00131 *> = 0: successful exit 00132 *> = 1: inverse iteration did not converge; V is set to the 00133 *> last iterate. 00134 *> \endverbatim 00135 * 00136 * Authors: 00137 * ======== 00138 * 00139 *> \author Univ. of Tennessee 00140 *> \author Univ. of California Berkeley 00141 *> \author Univ. of Colorado Denver 00142 *> \author NAG Ltd. 00143 * 00144 *> \date November 2011 00145 * 00146 *> \ingroup complex16OTHERauxiliary 00147 * 00148 * ===================================================================== 00149 SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, 00150 $ EPS3, SMLNUM, INFO ) 00151 * 00152 * -- LAPACK auxiliary routine (version 3.4.0) -- 00153 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00155 * November 2011 00156 * 00157 * .. Scalar Arguments .. 00158 LOGICAL NOINIT, RIGHTV 00159 INTEGER INFO, LDB, LDH, N 00160 DOUBLE PRECISION EPS3, SMLNUM 00161 COMPLEX*16 W 00162 * .. 00163 * .. Array Arguments .. 00164 DOUBLE PRECISION RWORK( * ) 00165 COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * ) 00166 * .. 00167 * 00168 * ===================================================================== 00169 * 00170 * .. Parameters .. 00171 DOUBLE PRECISION ONE, TENTH 00172 PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 ) 00173 COMPLEX*16 ZERO 00174 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 00175 * .. 00176 * .. Local Scalars .. 00177 CHARACTER NORMIN, TRANS 00178 INTEGER I, IERR, ITS, J 00179 DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM 00180 COMPLEX*16 CDUM, EI, EJ, TEMP, X 00181 * .. 00182 * .. External Functions .. 00183 INTEGER IZAMAX 00184 DOUBLE PRECISION DZASUM, DZNRM2 00185 COMPLEX*16 ZLADIV 00186 EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV 00187 * .. 00188 * .. External Subroutines .. 00189 EXTERNAL ZDSCAL, ZLATRS 00190 * .. 00191 * .. Intrinsic Functions .. 00192 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 00193 * .. 00194 * .. Statement Functions .. 00195 DOUBLE PRECISION CABS1 00196 * .. 00197 * .. Statement Function definitions .. 00198 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 00199 * .. 00200 * .. Executable Statements .. 00201 * 00202 INFO = 0 00203 * 00204 * GROWTO is the threshold used in the acceptance test for an 00205 * eigenvector. 00206 * 00207 ROOTN = SQRT( DBLE( N ) ) 00208 GROWTO = TENTH / ROOTN 00209 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM 00210 * 00211 * Form B = H - W*I (except that the subdiagonal elements are not 00212 * stored). 00213 * 00214 DO 20 J = 1, N 00215 DO 10 I = 1, J - 1 00216 B( I, J ) = H( I, J ) 00217 10 CONTINUE 00218 B( J, J ) = H( J, J ) - W 00219 20 CONTINUE 00220 * 00221 IF( NOINIT ) THEN 00222 * 00223 * Initialize V. 00224 * 00225 DO 30 I = 1, N 00226 V( I ) = EPS3 00227 30 CONTINUE 00228 ELSE 00229 * 00230 * Scale supplied initial vector. 00231 * 00232 VNORM = DZNRM2( N, V, 1 ) 00233 CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 ) 00234 END IF 00235 * 00236 IF( RIGHTV ) THEN 00237 * 00238 * LU decomposition with partial pivoting of B, replacing zero 00239 * pivots by EPS3. 00240 * 00241 DO 60 I = 1, N - 1 00242 EI = H( I+1, I ) 00243 IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN 00244 * 00245 * Interchange rows and eliminate. 00246 * 00247 X = ZLADIV( B( I, I ), EI ) 00248 B( I, I ) = EI 00249 DO 40 J = I + 1, N 00250 TEMP = B( I+1, J ) 00251 B( I+1, J ) = B( I, J ) - X*TEMP 00252 B( I, J ) = TEMP 00253 40 CONTINUE 00254 ELSE 00255 * 00256 * Eliminate without interchange. 00257 * 00258 IF( B( I, I ).EQ.ZERO ) 00259 $ B( I, I ) = EPS3 00260 X = ZLADIV( EI, B( I, I ) ) 00261 IF( X.NE.ZERO ) THEN 00262 DO 50 J = I + 1, N 00263 B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 00264 50 CONTINUE 00265 END IF 00266 END IF 00267 60 CONTINUE 00268 IF( B( N, N ).EQ.ZERO ) 00269 $ B( N, N ) = EPS3 00270 * 00271 TRANS = 'N' 00272 * 00273 ELSE 00274 * 00275 * UL decomposition with partial pivoting of B, replacing zero 00276 * pivots by EPS3. 00277 * 00278 DO 90 J = N, 2, -1 00279 EJ = H( J, J-1 ) 00280 IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN 00281 * 00282 * Interchange columns and eliminate. 00283 * 00284 X = ZLADIV( B( J, J ), EJ ) 00285 B( J, J ) = EJ 00286 DO 70 I = 1, J - 1 00287 TEMP = B( I, J-1 ) 00288 B( I, J-1 ) = B( I, J ) - X*TEMP 00289 B( I, J ) = TEMP 00290 70 CONTINUE 00291 ELSE 00292 * 00293 * Eliminate without interchange. 00294 * 00295 IF( B( J, J ).EQ.ZERO ) 00296 $ B( J, J ) = EPS3 00297 X = ZLADIV( EJ, B( J, J ) ) 00298 IF( X.NE.ZERO ) THEN 00299 DO 80 I = 1, J - 1 00300 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 00301 80 CONTINUE 00302 END IF 00303 END IF 00304 90 CONTINUE 00305 IF( B( 1, 1 ).EQ.ZERO ) 00306 $ B( 1, 1 ) = EPS3 00307 * 00308 TRANS = 'C' 00309 * 00310 END IF 00311 * 00312 NORMIN = 'N' 00313 DO 110 ITS = 1, N 00314 * 00315 * Solve U*x = scale*v for a right eigenvector 00316 * or U**H *x = scale*v for a left eigenvector, 00317 * overwriting x on v. 00318 * 00319 CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V, 00320 $ SCALE, RWORK, IERR ) 00321 NORMIN = 'Y' 00322 * 00323 * Test for sufficient growth in the norm of v. 00324 * 00325 VNORM = DZASUM( N, V, 1 ) 00326 IF( VNORM.GE.GROWTO*SCALE ) 00327 $ GO TO 120 00328 * 00329 * Choose new orthogonal starting vector and try again. 00330 * 00331 RTEMP = EPS3 / ( ROOTN+ONE ) 00332 V( 1 ) = EPS3 00333 DO 100 I = 2, N 00334 V( I ) = RTEMP 00335 100 CONTINUE 00336 V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN 00337 110 CONTINUE 00338 * 00339 * Failure to find eigenvector in N iterations. 00340 * 00341 INFO = 1 00342 * 00343 120 CONTINUE 00344 * 00345 * Normalize eigenvector. 00346 * 00347 I = IZAMAX( N, V, 1 ) 00348 CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 ) 00349 * 00350 RETURN 00351 * 00352 * End of ZLAEIN 00353 * 00354 END