|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZLATM6( TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, 00002 $ BETA, WX, WY, S, DIF ) 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 INTEGER LDA, LDX, LDY, N, TYPE 00010 COMPLEX*16 ALPHA, BETA, WX, WY 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION DIF( * ), S( * ) 00014 COMPLEX*16 A( LDA, * ), B( LDA, * ), X( LDX, * ), 00015 $ Y( LDY, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZLATM6 generates test matrices for the generalized eigenvalue 00022 * problem, their corresponding right and left eigenvector matrices, 00023 * and also reciprocal condition numbers for all eigenvalues and 00024 * the reciprocal condition numbers of eigenvectors corresponding to 00025 * the 1th and 5th eigenvalues. 00026 * 00027 * Test Matrices 00028 * ============= 00029 * 00030 * Two kinds of test matrix pairs 00031 * (A, B) = inverse(YH) * (Da, Db) * inverse(X) 00032 * are used in the tests: 00033 * 00034 * Type 1: 00035 * Da = 1+a 0 0 0 0 Db = 1 0 0 0 0 00036 * 0 2+a 0 0 0 0 1 0 0 0 00037 * 0 0 3+a 0 0 0 0 1 0 0 00038 * 0 0 0 4+a 0 0 0 0 1 0 00039 * 0 0 0 0 5+a , 0 0 0 0 1 00040 * and Type 2: 00041 * Da = 1+i 0 0 0 0 Db = 1 0 0 0 0 00042 * 0 1-i 0 0 0 0 1 0 0 0 00043 * 0 0 1 0 0 0 0 1 0 0 00044 * 0 0 0 (1+a)+(1+b)i 0 0 0 0 1 0 00045 * 0 0 0 0 (1+a)-(1+b)i, 0 0 0 0 1 . 00046 * 00047 * In both cases the same inverse(YH) and inverse(X) are used to compute 00048 * (A, B), giving the exact eigenvectors to (A,B) as (YH, X): 00049 * 00050 * YH: = 1 0 -y y -y X = 1 0 -x -x x 00051 * 0 1 -y y -y 0 1 x -x -x 00052 * 0 0 1 0 0 0 0 1 0 0 00053 * 0 0 0 1 0 0 0 0 1 0 00054 * 0 0 0 0 1, 0 0 0 0 1 , where 00055 * 00056 * a, b, x and y will have all values independently of each other. 00057 * 00058 * Arguments 00059 * ========= 00060 * 00061 * TYPE (input) INTEGER 00062 * Specifies the problem type (see futher details). 00063 * 00064 * N (input) INTEGER 00065 * Size of the matrices A and B. 00066 * 00067 * A (output) COMPLEX*16 array, dimension (LDA, N). 00068 * On exit A N-by-N is initialized according to TYPE. 00069 * 00070 * LDA (input) INTEGER 00071 * The leading dimension of A and of B. 00072 * 00073 * B (output) COMPLEX*16 array, dimension (LDA, N). 00074 * On exit B N-by-N is initialized according to TYPE. 00075 * 00076 * X (output) COMPLEX*16 array, dimension (LDX, N). 00077 * On exit X is the N-by-N matrix of right eigenvectors. 00078 * 00079 * LDX (input) INTEGER 00080 * The leading dimension of X. 00081 * 00082 * Y (output) COMPLEX*16 array, dimension (LDY, N). 00083 * On exit Y is the N-by-N matrix of left eigenvectors. 00084 * 00085 * LDY (input) INTEGER 00086 * The leading dimension of Y. 00087 * 00088 * ALPHA (input) COMPLEX*16 00089 * BETA (input) COMPLEX*16 00090 * Weighting constants for matrix A. 00091 * 00092 * WX (input) COMPLEX*16 00093 * Constant for right eigenvector matrix. 00094 * 00095 * WY (input) COMPLEX*16 00096 * Constant for left eigenvector matrix. 00097 * 00098 * S (output) DOUBLE PRECISION array, dimension (N) 00099 * S(i) is the reciprocal condition number for eigenvalue i. 00100 * 00101 * DIF (output) DOUBLE PRECISION array, dimension (N) 00102 * DIF(i) is the reciprocal condition number for eigenvector i. 00103 * 00104 * ===================================================================== 00105 * 00106 * .. Parameters .. 00107 DOUBLE PRECISION RONE, TWO, THREE 00108 PARAMETER ( RONE = 1.0D+0, TWO = 2.0D+0, THREE = 3.0D+0 ) 00109 COMPLEX*16 ZERO, ONE 00110 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 00111 $ ONE = ( 1.0D+0, 0.0D+0 ) ) 00112 * .. 00113 * .. Local Scalars .. 00114 INTEGER I, INFO, J 00115 * .. 00116 * .. Local Arrays .. 00117 DOUBLE PRECISION RWORK( 50 ) 00118 COMPLEX*16 WORK( 26 ), Z( 8, 8 ) 00119 * .. 00120 * .. Intrinsic Functions .. 00121 INTRINSIC CDABS, DBLE, DCMPLX, DCONJG, SQRT 00122 * .. 00123 * .. External Subroutines .. 00124 EXTERNAL ZGESVD, ZLACPY, ZLAKF2 00125 * .. 00126 * .. Executable Statements .. 00127 * 00128 * Generate test problem ... 00129 * (Da, Db) ... 00130 * 00131 DO 20 I = 1, N 00132 DO 10 J = 1, N 00133 * 00134 IF( I.EQ.J ) THEN 00135 A( I, I ) = DCMPLX( I ) + ALPHA 00136 B( I, I ) = ONE 00137 ELSE 00138 A( I, J ) = ZERO 00139 B( I, J ) = ZERO 00140 END IF 00141 * 00142 10 CONTINUE 00143 20 CONTINUE 00144 IF( TYPE.EQ.2 ) THEN 00145 A( 1, 1 ) = DCMPLX( RONE, RONE ) 00146 A( 2, 2 ) = DCONJG( A( 1, 1 ) ) 00147 A( 3, 3 ) = ONE 00148 A( 4, 4 ) = DCMPLX( DBLE( ONE+ALPHA ), DBLE( ONE+BETA ) ) 00149 A( 5, 5 ) = DCONJG( A( 4, 4 ) ) 00150 END IF 00151 * 00152 * Form X and Y 00153 * 00154 CALL ZLACPY( 'F', N, N, B, LDA, Y, LDY ) 00155 Y( 3, 1 ) = -DCONJG( WY ) 00156 Y( 4, 1 ) = DCONJG( WY ) 00157 Y( 5, 1 ) = -DCONJG( WY ) 00158 Y( 3, 2 ) = -DCONJG( WY ) 00159 Y( 4, 2 ) = DCONJG( WY ) 00160 Y( 5, 2 ) = -DCONJG( WY ) 00161 * 00162 CALL ZLACPY( 'F', N, N, B, LDA, X, LDX ) 00163 X( 1, 3 ) = -WX 00164 X( 1, 4 ) = -WX 00165 X( 1, 5 ) = WX 00166 X( 2, 3 ) = WX 00167 X( 2, 4 ) = -WX 00168 X( 2, 5 ) = -WX 00169 * 00170 * Form (A, B) 00171 * 00172 B( 1, 3 ) = WX + WY 00173 B( 2, 3 ) = -WX + WY 00174 B( 1, 4 ) = WX - WY 00175 B( 2, 4 ) = WX - WY 00176 B( 1, 5 ) = -WX + WY 00177 B( 2, 5 ) = WX + WY 00178 A( 1, 3 ) = WX*A( 1, 1 ) + WY*A( 3, 3 ) 00179 A( 2, 3 ) = -WX*A( 2, 2 ) + WY*A( 3, 3 ) 00180 A( 1, 4 ) = WX*A( 1, 1 ) - WY*A( 4, 4 ) 00181 A( 2, 4 ) = WX*A( 2, 2 ) - WY*A( 4, 4 ) 00182 A( 1, 5 ) = -WX*A( 1, 1 ) + WY*A( 5, 5 ) 00183 A( 2, 5 ) = WX*A( 2, 2 ) + WY*A( 5, 5 ) 00184 * 00185 * Compute condition numbers 00186 * 00187 S( 1 ) = RONE / SQRT( ( RONE+THREE*CDABS( WY )*CDABS( WY ) ) / 00188 $ ( RONE+CDABS( A( 1, 1 ) )*CDABS( A( 1, 1 ) ) ) ) 00189 S( 2 ) = RONE / SQRT( ( RONE+THREE*CDABS( WY )*CDABS( WY ) ) / 00190 $ ( RONE+CDABS( A( 2, 2 ) )*CDABS( A( 2, 2 ) ) ) ) 00191 S( 3 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) / 00192 $ ( RONE+CDABS( A( 3, 3 ) )*CDABS( A( 3, 3 ) ) ) ) 00193 S( 4 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) / 00194 $ ( RONE+CDABS( A( 4, 4 ) )*CDABS( A( 4, 4 ) ) ) ) 00195 S( 5 ) = RONE / SQRT( ( RONE+TWO*CDABS( WX )*CDABS( WX ) ) / 00196 $ ( RONE+CDABS( A( 5, 5 ) )*CDABS( A( 5, 5 ) ) ) ) 00197 * 00198 CALL ZLAKF2( 1, 4, A, LDA, A( 2, 2 ), B, B( 2, 2 ), Z, 8 ) 00199 CALL ZGESVD( 'N', 'N', 8, 8, Z, 8, RWORK, WORK, 1, WORK( 2 ), 1, 00200 $ WORK( 3 ), 24, RWORK( 9 ), INFO ) 00201 DIF( 1 ) = RWORK( 8 ) 00202 * 00203 CALL ZLAKF2( 4, 1, A, LDA, A( 5, 5 ), B, B( 5, 5 ), Z, 8 ) 00204 CALL ZGESVD( 'N', 'N', 8, 8, Z, 8, RWORK, WORK, 1, WORK( 2 ), 1, 00205 $ WORK( 3 ), 24, RWORK( 9 ), INFO ) 00206 DIF( 5 ) = RWORK( 8 ) 00207 * 00208 RETURN 00209 * 00210 * End of ZLATM6 00211 * 00212 END
1.7.3 
