|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, 00002 $ WORK, RWORK, RESULT ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 LOGICAL LEFT 00010 INTEGER LDA, LDB, LDE, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION RESULT( 2 ), RWORK( * ) 00014 COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), 00015 $ BETA( * ), E( LDE, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZGET52 does an eigenvector check for the generalized eigenvalue 00022 * problem. 00023 * 00024 * The basic test for right eigenvectors is: 00025 * 00026 * | b(i) A E(i) - a(i) B E(i) | 00027 * RESULT(1) = max ------------------------------- 00028 * i n ulp max( |b(i) A|, |a(i) B| ) 00029 * 00030 * using the 1-norm. Here, a(i)/b(i) = w is the i-th generalized 00031 * eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th 00032 * generalized eigenvalue of m A - B. 00033 * 00034 * H H _ _ 00035 * For left eigenvectors, A , B , a, and b are used. 00036 * 00037 * ZGET52 also tests the normalization of E. Each eigenvector is 00038 * supposed to be normalized so that the maximum "absolute value" 00039 * of its elements is 1, where in this case, "absolute value" 00040 * of a complex value x is |Re(x)| + |Im(x)| ; let us call this 00041 * maximum "absolute value" norm of a vector v M(v). 00042 * If a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate 00043 * vector. The normalization test is: 00044 * 00045 * RESULT(2) = max | M(v(i)) - 1 | / ( n ulp ) 00046 * eigenvectors v(i) 00047 * 00048 * 00049 * Arguments 00050 * ========= 00051 * 00052 * LEFT (input) LOGICAL 00053 * =.TRUE.: The eigenvectors in the columns of E are assumed 00054 * to be *left* eigenvectors. 00055 * =.FALSE.: The eigenvectors in the columns of E are assumed 00056 * to be *right* eigenvectors. 00057 * 00058 * N (input) INTEGER 00059 * The size of the matrices. If it is zero, ZGET52 does 00060 * nothing. It must be at least zero. 00061 * 00062 * A (input) COMPLEX*16 array, dimension (LDA, N) 00063 * The matrix A. 00064 * 00065 * LDA (input) INTEGER 00066 * The leading dimension of A. It must be at least 1 00067 * and at least N. 00068 * 00069 * B (input) COMPLEX*16 array, dimension (LDB, N) 00070 * The matrix B. 00071 * 00072 * LDB (input) INTEGER 00073 * The leading dimension of B. It must be at least 1 00074 * and at least N. 00075 * 00076 * E (input) COMPLEX*16 array, dimension (LDE, N) 00077 * The matrix of eigenvectors. It must be O( 1 ). 00078 * 00079 * LDE (input) INTEGER 00080 * The leading dimension of E. It must be at least 1 and at 00081 * least N. 00082 * 00083 * ALPHA (input) COMPLEX*16 array, dimension (N) 00084 * The values a(i) as described above, which, along with b(i), 00085 * define the generalized eigenvalues. 00086 * 00087 * BETA (input) COMPLEX*16 array, dimension (N) 00088 * The values b(i) as described above, which, along with a(i), 00089 * define the generalized eigenvalues. 00090 * 00091 * WORK (workspace) COMPLEX*16 array, dimension (N**2) 00092 * 00093 * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 00094 * 00095 * RESULT (output) DOUBLE PRECISION array, dimension (2) 00096 * The values computed by the test described above. If A E or 00097 * B E is likely to overflow, then RESULT(1:2) is set to 00098 * 10 / ulp. 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Parameters .. 00103 DOUBLE PRECISION ZERO, ONE 00104 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00105 COMPLEX*16 CZERO, CONE 00106 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00107 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00108 * .. 00109 * .. Local Scalars .. 00110 CHARACTER NORMAB, TRANS 00111 INTEGER J, JVEC 00112 DOUBLE PRECISION ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM, 00113 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1, 00114 $ ULP 00115 COMPLEX*16 ACOEFF, ALPHAI, BCOEFF, BETAI, X 00116 * .. 00117 * .. External Functions .. 00118 DOUBLE PRECISION DLAMCH, ZLANGE 00119 EXTERNAL DLAMCH, ZLANGE 00120 * .. 00121 * .. External Subroutines .. 00122 EXTERNAL ZGEMV 00123 * .. 00124 * .. Intrinsic Functions .. 00125 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX 00126 * .. 00127 * .. Statement Functions .. 00128 DOUBLE PRECISION ABS1 00129 * .. 00130 * .. Statement Function definitions .. 00131 ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) ) 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 RESULT( 1 ) = ZERO 00136 RESULT( 2 ) = ZERO 00137 IF( N.LE.0 ) 00138 $ RETURN 00139 * 00140 SAFMIN = DLAMCH( 'Safe minimum' ) 00141 SAFMAX = ONE / SAFMIN 00142 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00143 * 00144 IF( LEFT ) THEN 00145 TRANS = 'C' 00146 NORMAB = 'I' 00147 ELSE 00148 TRANS = 'N' 00149 NORMAB = 'O' 00150 END IF 00151 * 00152 * Norm of A, B, and E: 00153 * 00154 ANORM = MAX( ZLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN ) 00155 BNORM = MAX( ZLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN ) 00156 ENORM = MAX( ZLANGE( 'O', N, N, E, LDE, RWORK ), ULP ) 00157 ALFMAX = SAFMAX / MAX( ONE, BNORM ) 00158 BETMAX = SAFMAX / MAX( ONE, ANORM ) 00159 * 00160 * Compute error matrix. 00161 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) 00162 * 00163 DO 10 JVEC = 1, N 00164 ALPHAI = ALPHA( JVEC ) 00165 BETAI = BETA( JVEC ) 00166 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) ) 00167 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR. 00168 $ ABMAX.LT.ONE ) THEN 00169 SCALE = ONE / MAX( ABMAX, SAFMIN ) 00170 ALPHAI = SCALE*ALPHAI 00171 BETAI = SCALE*BETAI 00172 END IF 00173 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM, 00174 $ SAFMIN ) 00175 ACOEFF = SCALE*BETAI 00176 BCOEFF = SCALE*ALPHAI 00177 IF( LEFT ) THEN 00178 ACOEFF = DCONJG( ACOEFF ) 00179 BCOEFF = DCONJG( BCOEFF ) 00180 END IF 00181 CALL ZGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1, 00182 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 ) 00183 CALL ZGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1, 00184 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00185 10 CONTINUE 00186 * 00187 ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM 00188 * 00189 * Compute RESULT(1) 00190 * 00191 RESULT( 1 ) = ERRNRM / ULP 00192 * 00193 * Normalization of E: 00194 * 00195 ENRMER = ZERO 00196 DO 30 JVEC = 1, N 00197 TEMP1 = ZERO 00198 DO 20 J = 1, N 00199 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) ) 00200 20 CONTINUE 00201 ENRMER = MAX( ENRMER, TEMP1-ONE ) 00202 30 CONTINUE 00203 * 00204 * Compute RESULT(2) : the normalization error in E. 00205 * 00206 RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP ) 00207 * 00208 RETURN 00209 * 00210 * End of ZGET52 00211 * 00212 END
1.7.3 
