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