1 SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INFO, LDA, LDB
9 DOUBLE PRECISION RESULT, SCALE, WI, WR
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGET53 checks the generalized eigenvalues computed by DLAG2.
19 *
20 * The basic test for an eigenvalue is:
21 *
22 * | det( s A - w B ) |
23 * RESULT = ---------------------------------------------------
24 * ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
25 *
26 * Two "safety checks" are performed:
27 *
28 * (1) ulp*max( s*norm(A), |w|*norm(B) ) must be at least
29 * safe_minimum. This insures that the test performed is
30 * not essentially det(0*A + 0*B)=0.
31 *
32 * (2) s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
33 * This insures that s*A - w*B will not overflow.
34 *
35 * If these tests are not passed, then s and w are scaled and
36 * tested anyway, if this is possible.
37 *
38 * Arguments
39 * =========
40 *
41 * A (input) DOUBLE PRECISION array, dimension (LDA, 2)
42 * The 2x2 matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of A. It must be at least 2.
46 *
47 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
48 * The 2x2 upper-triangular matrix B.
49 *
50 * LDB (input) INTEGER
51 * The leading dimension of B. It must be at least 2.
52 *
53 * SCALE (input) DOUBLE PRECISION
54 * The "scale factor" s in the formula s A - w B . It is
55 * assumed to be non-negative.
56 *
57 * WR (input) DOUBLE PRECISION
58 * The real part of the eigenvalue w in the formula
59 * s A - w B .
60 *
61 * WI (input) DOUBLE PRECISION
62 * The imaginary part of the eigenvalue w in the formula
63 * s A - w B .
64 *
65 * RESULT (output) DOUBLE PRECISION
66 * If INFO is 2 or less, the value computed by the test
67 * described above.
68 * If INFO=3, this will just be 1/ulp.
69 *
70 * INFO (output) INTEGER
71 * =0: The input data pass the "safety checks".
72 * =1: s*norm(A) + |w|*norm(B) > 1/safe_minimum.
73 * =2: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
74 * =3: same as INFO=2, but s and w could not be scaled so
75 * as to compute the test.
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE
81 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
82 * ..
83 * .. Local Scalars ..
84 DOUBLE PRECISION ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
85 $ CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
86 $ SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
87 * ..
88 * .. External Functions ..
89 DOUBLE PRECISION DLAMCH
90 EXTERNAL DLAMCH
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC ABS, MAX, SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Initialize
98 *
99 INFO = 0
100 RESULT = ZERO
101 SCALES = SCALE
102 WRS = WR
103 WIS = WI
104 *
105 * Machine constants and norms
106 *
107 SAFMIN = DLAMCH( 'Safe minimum' )
108 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
109 ABSW = ABS( WRS ) + ABS( WIS )
110 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
111 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
112 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
113 $ SAFMIN )
114 *
115 * Check for possible overflow.
116 *
117 TEMP = ( SAFMIN*BNORM )*ABSW + ( SAFMIN*ANORM )*SCALES
118 IF( TEMP.GE.ONE ) THEN
119 *
120 * Scale down to avoid overflow
121 *
122 INFO = 1
123 TEMP = ONE / TEMP
124 SCALES = SCALES*TEMP
125 WRS = WRS*TEMP
126 WIS = WIS*TEMP
127 ABSW = ABS( WRS ) + ABS( WIS )
128 END IF
129 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
130 $ SAFMIN*MAX( SCALES, ABSW ) )
131 *
132 * Check for W and SCALE essentially zero.
133 *
134 IF( S1.LT.SAFMIN ) THEN
135 INFO = 2
136 IF( SCALES.LT.SAFMIN .AND. ABSW.LT.SAFMIN ) THEN
137 INFO = 3
138 RESULT = ONE / ULP
139 RETURN
140 END IF
141 *
142 * Scale up to avoid underflow
143 *
144 TEMP = ONE / MAX( SCALES*ANORM+ABSW*BNORM, SAFMIN )
145 SCALES = SCALES*TEMP
146 WRS = WRS*TEMP
147 WIS = WIS*TEMP
148 ABSW = ABS( WRS ) + ABS( WIS )
149 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
150 $ SAFMIN*MAX( SCALES, ABSW ) )
151 IF( S1.LT.SAFMIN ) THEN
152 INFO = 3
153 RESULT = ONE / ULP
154 RETURN
155 END IF
156 END IF
157 *
158 * Compute C = s A - w B
159 *
160 CR11 = SCALES*A( 1, 1 ) - WRS*B( 1, 1 )
161 CI11 = -WIS*B( 1, 1 )
162 CR21 = SCALES*A( 2, 1 )
163 CR12 = SCALES*A( 1, 2 ) - WRS*B( 1, 2 )
164 CI12 = -WIS*B( 1, 2 )
165 CR22 = SCALES*A( 2, 2 ) - WRS*B( 2, 2 )
166 CI22 = -WIS*B( 2, 2 )
167 *
168 * Compute the smallest singular value of s A - w B:
169 *
170 * |det( s A - w B )|
171 * sigma_min = ------------------
172 * norm( s A - w B )
173 *
174 CNORM = MAX( ABS( CR11 )+ABS( CI11 )+ABS( CR21 ),
175 $ ABS( CR12 )+ABS( CI12 )+ABS( CR22 )+ABS( CI22 ), SAFMIN )
176 CSCALE = ONE / SQRT( CNORM )
177 DETR = ( CSCALE*CR11 )*( CSCALE*CR22 ) -
178 $ ( CSCALE*CI11 )*( CSCALE*CI22 ) -
179 $ ( CSCALE*CR12 )*( CSCALE*CR21 )
180 DETI = ( CSCALE*CR11 )*( CSCALE*CI22 ) +
181 $ ( CSCALE*CI11 )*( CSCALE*CR22 ) -
182 $ ( CSCALE*CI12 )*( CSCALE*CR21 )
183 SIGMIN = ABS( DETR ) + ABS( DETI )
184 RESULT = SIGMIN / S1
185 RETURN
186 *
187 * End of DGET53
188 *
189 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INFO, LDA, LDB
9 DOUBLE PRECISION RESULT, SCALE, WI, WR
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DGET53 checks the generalized eigenvalues computed by DLAG2.
19 *
20 * The basic test for an eigenvalue is:
21 *
22 * | det( s A - w B ) |
23 * RESULT = ---------------------------------------------------
24 * ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
25 *
26 * Two "safety checks" are performed:
27 *
28 * (1) ulp*max( s*norm(A), |w|*norm(B) ) must be at least
29 * safe_minimum. This insures that the test performed is
30 * not essentially det(0*A + 0*B)=0.
31 *
32 * (2) s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
33 * This insures that s*A - w*B will not overflow.
34 *
35 * If these tests are not passed, then s and w are scaled and
36 * tested anyway, if this is possible.
37 *
38 * Arguments
39 * =========
40 *
41 * A (input) DOUBLE PRECISION array, dimension (LDA, 2)
42 * The 2x2 matrix A.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of A. It must be at least 2.
46 *
47 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
48 * The 2x2 upper-triangular matrix B.
49 *
50 * LDB (input) INTEGER
51 * The leading dimension of B. It must be at least 2.
52 *
53 * SCALE (input) DOUBLE PRECISION
54 * The "scale factor" s in the formula s A - w B . It is
55 * assumed to be non-negative.
56 *
57 * WR (input) DOUBLE PRECISION
58 * The real part of the eigenvalue w in the formula
59 * s A - w B .
60 *
61 * WI (input) DOUBLE PRECISION
62 * The imaginary part of the eigenvalue w in the formula
63 * s A - w B .
64 *
65 * RESULT (output) DOUBLE PRECISION
66 * If INFO is 2 or less, the value computed by the test
67 * described above.
68 * If INFO=3, this will just be 1/ulp.
69 *
70 * INFO (output) INTEGER
71 * =0: The input data pass the "safety checks".
72 * =1: s*norm(A) + |w|*norm(B) > 1/safe_minimum.
73 * =2: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
74 * =3: same as INFO=2, but s and w could not be scaled so
75 * as to compute the test.
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE
81 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
82 * ..
83 * .. Local Scalars ..
84 DOUBLE PRECISION ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
85 $ CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
86 $ SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
87 * ..
88 * .. External Functions ..
89 DOUBLE PRECISION DLAMCH
90 EXTERNAL DLAMCH
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC ABS, MAX, SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Initialize
98 *
99 INFO = 0
100 RESULT = ZERO
101 SCALES = SCALE
102 WRS = WR
103 WIS = WI
104 *
105 * Machine constants and norms
106 *
107 SAFMIN = DLAMCH( 'Safe minimum' )
108 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
109 ABSW = ABS( WRS ) + ABS( WIS )
110 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ),
111 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN )
112 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ),
113 $ SAFMIN )
114 *
115 * Check for possible overflow.
116 *
117 TEMP = ( SAFMIN*BNORM )*ABSW + ( SAFMIN*ANORM )*SCALES
118 IF( TEMP.GE.ONE ) THEN
119 *
120 * Scale down to avoid overflow
121 *
122 INFO = 1
123 TEMP = ONE / TEMP
124 SCALES = SCALES*TEMP
125 WRS = WRS*TEMP
126 WIS = WIS*TEMP
127 ABSW = ABS( WRS ) + ABS( WIS )
128 END IF
129 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
130 $ SAFMIN*MAX( SCALES, ABSW ) )
131 *
132 * Check for W and SCALE essentially zero.
133 *
134 IF( S1.LT.SAFMIN ) THEN
135 INFO = 2
136 IF( SCALES.LT.SAFMIN .AND. ABSW.LT.SAFMIN ) THEN
137 INFO = 3
138 RESULT = ONE / ULP
139 RETURN
140 END IF
141 *
142 * Scale up to avoid underflow
143 *
144 TEMP = ONE / MAX( SCALES*ANORM+ABSW*BNORM, SAFMIN )
145 SCALES = SCALES*TEMP
146 WRS = WRS*TEMP
147 WIS = WIS*TEMP
148 ABSW = ABS( WRS ) + ABS( WIS )
149 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ),
150 $ SAFMIN*MAX( SCALES, ABSW ) )
151 IF( S1.LT.SAFMIN ) THEN
152 INFO = 3
153 RESULT = ONE / ULP
154 RETURN
155 END IF
156 END IF
157 *
158 * Compute C = s A - w B
159 *
160 CR11 = SCALES*A( 1, 1 ) - WRS*B( 1, 1 )
161 CI11 = -WIS*B( 1, 1 )
162 CR21 = SCALES*A( 2, 1 )
163 CR12 = SCALES*A( 1, 2 ) - WRS*B( 1, 2 )
164 CI12 = -WIS*B( 1, 2 )
165 CR22 = SCALES*A( 2, 2 ) - WRS*B( 2, 2 )
166 CI22 = -WIS*B( 2, 2 )
167 *
168 * Compute the smallest singular value of s A - w B:
169 *
170 * |det( s A - w B )|
171 * sigma_min = ------------------
172 * norm( s A - w B )
173 *
174 CNORM = MAX( ABS( CR11 )+ABS( CI11 )+ABS( CR21 ),
175 $ ABS( CR12 )+ABS( CI12 )+ABS( CR22 )+ABS( CI22 ), SAFMIN )
176 CSCALE = ONE / SQRT( CNORM )
177 DETR = ( CSCALE*CR11 )*( CSCALE*CR22 ) -
178 $ ( CSCALE*CI11 )*( CSCALE*CI22 ) -
179 $ ( CSCALE*CR12 )*( CSCALE*CR21 )
180 DETI = ( CSCALE*CR11 )*( CSCALE*CI22 ) +
181 $ ( CSCALE*CI11 )*( CSCALE*CR22 ) -
182 $ ( CSCALE*CI12 )*( CSCALE*CR21 )
183 SIGMIN = ABS( DETR ) + ABS( DETI )
184 RESULT = SIGMIN / S1
185 RETURN
186 *
187 * End of DGET53
188 *
189 END