 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGET52( 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 REAL RESULT( 2 ), RWORK( * ) 00014 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00015 $ BETA( * ), E( LDE, * ), WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CGET52 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 * CGET52 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 * Arguments 00049 * ========= 00050 * 00051 * LEFT (input) LOGICAL 00052 * =.TRUE.: The eigenvectors in the columns of E are assumed 00053 * to be *left* eigenvectors. 00054 * =.FALSE.: The eigenvectors in the columns of E are assumed 00055 * to be *right* eigenvectors. 00056 * 00057 * N (input) INTEGER 00058 * The size of the matrices. If it is zero, CGET52 does 00059 * nothing. It must be at least zero. 00060 * 00061 * A (input) COMPLEX array, dimension (LDA, N) 00062 * The matrix A. 00063 * 00064 * LDA (input) INTEGER 00065 * The leading dimension of A. It must be at least 1 00066 * and at least N. 00067 * 00068 * B (input) COMPLEX array, dimension (LDB, N) 00069 * The matrix B. 00070 * 00071 * LDB (input) INTEGER 00072 * The leading dimension of B. It must be at least 1 00073 * and at least N. 00074 * 00075 * E (input) COMPLEX array, dimension (LDE, N) 00076 * The matrix of eigenvectors. It must be O( 1 ). 00077 * 00078 * LDE (input) INTEGER 00079 * The leading dimension of E. It must be at least 1 and at 00080 * least N. 00081 * 00082 * ALPHA (input) COMPLEX array, dimension (N) 00083 * The values a(i) as described above, which, along with b(i), 00084 * define the generalized eigenvalues. 00085 * 00086 * BETA (input) COMPLEX array, dimension (N) 00087 * The values b(i) as described above, which, along with a(i), 00088 * define the generalized eigenvalues. 00089 * 00090 * WORK (workspace) COMPLEX array, dimension (N**2) 00091 * 00092 * RWORK (workspace) REAL array, dimension (N) 00093 * 00094 * RESULT (output) REAL array, dimension (2) 00095 * The values computed by the test described above. If A E or 00096 * B E is likely to overflow, then RESULT(1:2) is set to 00097 * 10 / ulp. 00098 * 00099 * ===================================================================== 00100 * 00101 * .. Parameters .. 00102 REAL ZERO, ONE 00103 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00104 COMPLEX CZERO, CONE 00105 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00106 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00107 * .. 00108 * .. Local Scalars .. 00109 CHARACTER NORMAB, TRANS 00110 INTEGER J, JVEC 00111 REAL ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM, 00112 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1, 00113 $ ULP 00114 COMPLEX ACOEFF, ALPHAI, BCOEFF, BETAI, X 00115 * .. 00116 * .. External Functions .. 00117 REAL CLANGE, SLAMCH 00118 EXTERNAL CLANGE, SLAMCH 00119 * .. 00120 * .. External Subroutines .. 00121 EXTERNAL CGEMV 00122 * .. 00123 * .. Intrinsic Functions .. 00124 INTRINSIC ABS, AIMAG, CONJG, MAX, REAL 00125 * .. 00126 * .. Statement Functions .. 00127 REAL ABS1 00128 * .. 00129 * .. Statement Function definitions .. 00130 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00131 * .. 00132 * .. Executable Statements .. 00133 * 00134 RESULT( 1 ) = ZERO 00135 RESULT( 2 ) = ZERO 00136 IF( N.LE.0 ) 00137 $ RETURN 00138 * 00139 SAFMIN = SLAMCH( 'Safe minimum' ) 00140 SAFMAX = ONE / SAFMIN 00141 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00142 * 00143 IF( LEFT ) THEN 00144 TRANS = 'C' 00145 NORMAB = 'I' 00146 ELSE 00147 TRANS = 'N' 00148 NORMAB = 'O' 00149 END IF 00150 * 00151 * Norm of A, B, and E: 00152 * 00153 ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN ) 00154 BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN ) 00155 ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP ) 00156 ALFMAX = SAFMAX / MAX( ONE, BNORM ) 00157 BETMAX = SAFMAX / MAX( ONE, ANORM ) 00158 * 00159 * Compute error matrix. 00160 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) 00161 * 00162 DO 10 JVEC = 1, N 00163 ALPHAI = ALPHA( JVEC ) 00164 BETAI = BETA( JVEC ) 00165 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) ) 00166 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR. 00167 $ ABMAX.LT.ONE ) THEN 00168 SCALE = ONE / MAX( ABMAX, SAFMIN ) 00169 ALPHAI = SCALE*ALPHAI 00170 BETAI = SCALE*BETAI 00171 END IF 00172 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM, 00173 $ SAFMIN ) 00174 ACOEFF = SCALE*BETAI 00175 BCOEFF = SCALE*ALPHAI 00176 IF( LEFT ) THEN 00177 ACOEFF = CONJG( ACOEFF ) 00178 BCOEFF = CONJG( BCOEFF ) 00179 END IF 00180 CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1, 00181 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 ) 00182 CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1, 00183 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00184 10 CONTINUE 00185 * 00186 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM 00187 * 00188 * Compute RESULT(1) 00189 * 00190 RESULT( 1 ) = ERRNRM / ULP 00191 * 00192 * Normalization of E: 00193 * 00194 ENRMER = ZERO 00195 DO 30 JVEC = 1, N 00196 TEMP1 = ZERO 00197 DO 20 J = 1, N 00198 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) ) 00199 20 CONTINUE 00200 ENRMER = MAX( ENRMER, TEMP1-ONE ) 00201 30 CONTINUE 00202 * 00203 * Compute RESULT(2) : the normalization error in E. 00204 * 00205 RESULT( 2 ) = ENRMER / ( REAL( N )*ULP ) 00206 * 00207 RETURN 00208 * 00209 * End of CGET52 00210 * 00211 END
1.7.3 
