LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dget53.f
Go to the documentation of this file.
1*> \brief \b DGET53
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
12*
13* .. Scalar Arguments ..
14* INTEGER INFO, LDA, LDB
15* DOUBLE PRECISION RESULT, SCALE, WI, WR
16* ..
17* .. Array Arguments ..
18* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> DGET53 checks the generalized eigenvalues computed by DLAG2.
28*>
29*> The basic test for an eigenvalue is:
30*>
31*> | det( s A - w B ) |
32*> RESULT = ---------------------------------------------------
33*> ulp max( s norm(A), |w| norm(B) )*norm( s A - w B )
34*>
35*> Two "safety checks" are performed:
36*>
37*> (1) ulp*max( s*norm(A), |w|*norm(B) ) must be at least
38*> safe_minimum. This insures that the test performed is
39*> not essentially det(0*A + 0*B)=0.
40*>
41*> (2) s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum.
42*> This insures that s*A - w*B will not overflow.
43*>
44*> If these tests are not passed, then s and w are scaled and
45*> tested anyway, if this is possible.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] A
52*> \verbatim
53*> A is DOUBLE PRECISION array, dimension (LDA, 2)
54*> The 2x2 matrix A.
55*> \endverbatim
56*>
57*> \param[in] LDA
58*> \verbatim
59*> LDA is INTEGER
60*> The leading dimension of A. It must be at least 2.
61*> \endverbatim
62*>
63*> \param[in] B
64*> \verbatim
65*> B is DOUBLE PRECISION array, dimension (LDB, N)
66*> The 2x2 upper-triangular matrix B.
67*> \endverbatim
68*>
69*> \param[in] LDB
70*> \verbatim
71*> LDB is INTEGER
72*> The leading dimension of B. It must be at least 2.
73*> \endverbatim
74*>
75*> \param[in] SCALE
76*> \verbatim
77*> SCALE is DOUBLE PRECISION
78*> The "scale factor" s in the formula s A - w B . It is
79*> assumed to be non-negative.
80*> \endverbatim
81*>
82*> \param[in] WR
83*> \verbatim
84*> WR is DOUBLE PRECISION
85*> The real part of the eigenvalue w in the formula
86*> s A - w B .
87*> \endverbatim
88*>
89*> \param[in] WI
90*> \verbatim
91*> WI is DOUBLE PRECISION
92*> The imaginary part of the eigenvalue w in the formula
93*> s A - w B .
94*> \endverbatim
95*>
96*> \param[out] RESULT
97*> \verbatim
98*> RESULT is DOUBLE PRECISION
99*> If INFO is 2 or less, the value computed by the test
100*> described above.
101*> If INFO=3, this will just be 1/ulp.
102*> \endverbatim
103*>
104*> \param[out] INFO
105*> \verbatim
106*> INFO is INTEGER
107*> =0: The input data pass the "safety checks".
108*> =1: s*norm(A) + |w|*norm(B) > 1/safe_minimum.
109*> =2: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum
110*> =3: same as INFO=2, but s and w could not be scaled so
111*> as to compute the test.
112*> \endverbatim
113*
114* Authors:
115* ========
116*
117*> \author Univ. of Tennessee
118*> \author Univ. of California Berkeley
119*> \author Univ. of Colorado Denver
120*> \author NAG Ltd.
121*
122*> \ingroup double_eig
123*
124* =====================================================================
125 SUBROUTINE dget53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
126*
127* -- LAPACK test routine --
128* -- LAPACK is a software package provided by Univ. of Tennessee, --
129* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
130*
131* .. Scalar Arguments ..
132 INTEGER INFO, LDA, LDB
133 DOUBLE PRECISION RESULT, SCALE, WI, WR
134* ..
135* .. Array Arguments ..
136 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
137* ..
138*
139* =====================================================================
140*
141* .. Parameters ..
142 DOUBLE PRECISION ZERO, ONE
143 parameter( zero = 0.0d0, one = 1.0d0 )
144* ..
145* .. Local Scalars ..
146 DOUBLE PRECISION ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM,
147 $ CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1,
148 $ SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS
149* ..
150* .. External Functions ..
151 DOUBLE PRECISION DLAMCH
152 EXTERNAL dlamch
153* ..
154* .. Intrinsic Functions ..
155 INTRINSIC abs, max, sqrt
156* ..
157* .. Executable Statements ..
158*
159* Initialize
160*
161 info = 0
162 result = zero
163 scales = scale
164 wrs = wr
165 wis = wi
166*
167* Machine constants and norms
168*
169 safmin = dlamch( 'Safe minimum' )
170 ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
171 absw = abs( wrs ) + abs( wis )
172 anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
173 $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
174 bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
175 $ safmin )
176*
177* Check for possible overflow.
178*
179 temp = ( safmin*bnorm )*absw + ( safmin*anorm )*scales
180 IF( temp.GE.one ) THEN
181*
182* Scale down to avoid overflow
183*
184 info = 1
185 temp = one / temp
186 scales = scales*temp
187 wrs = wrs*temp
188 wis = wis*temp
189 absw = abs( wrs ) + abs( wis )
190 END IF
191 s1 = max( ulp*max( scales*anorm, absw*bnorm ),
192 $ safmin*max( scales, absw ) )
193*
194* Check for W and SCALE essentially zero.
195*
196 IF( s1.LT.safmin ) THEN
197 info = 2
198 IF( scales.LT.safmin .AND. absw.LT.safmin ) THEN
199 info = 3
200 result = one / ulp
201 RETURN
202 END IF
203*
204* Scale up to avoid underflow
205*
206 temp = one / max( scales*anorm+absw*bnorm, safmin )
207 scales = scales*temp
208 wrs = wrs*temp
209 wis = wis*temp
210 absw = abs( wrs ) + abs( wis )
211 s1 = max( ulp*max( scales*anorm, absw*bnorm ),
212 $ safmin*max( scales, absw ) )
213 IF( s1.LT.safmin ) THEN
214 info = 3
215 result = one / ulp
216 RETURN
217 END IF
218 END IF
219*
220* Compute C = s A - w B
221*
222 cr11 = scales*a( 1, 1 ) - wrs*b( 1, 1 )
223 ci11 = -wis*b( 1, 1 )
224 cr21 = scales*a( 2, 1 )
225 cr12 = scales*a( 1, 2 ) - wrs*b( 1, 2 )
226 ci12 = -wis*b( 1, 2 )
227 cr22 = scales*a( 2, 2 ) - wrs*b( 2, 2 )
228 ci22 = -wis*b( 2, 2 )
229*
230* Compute the smallest singular value of s A - w B:
231*
232* |det( s A - w B )|
233* sigma_min = ------------------
234* norm( s A - w B )
235*
236 cnorm = max( abs( cr11 )+abs( ci11 )+abs( cr21 ),
237 $ abs( cr12 )+abs( ci12 )+abs( cr22 )+abs( ci22 ), safmin )
238 cscale = one / sqrt( cnorm )
239 detr = ( cscale*cr11 )*( cscale*cr22 ) -
240 $ ( cscale*ci11 )*( cscale*ci22 ) -
241 $ ( cscale*cr12 )*( cscale*cr21 )
242 deti = ( cscale*cr11 )*( cscale*ci22 ) +
243 $ ( cscale*ci11 )*( cscale*cr22 ) -
244 $ ( cscale*ci12 )*( cscale*cr21 )
245 sigmin = abs( detr ) + abs( deti )
246 result = sigmin / s1
247 RETURN
248*
249* End of DGET53
250*
251 END
subroutine dget53(a, lda, b, ldb, scale, wr, wi, result, info)
DGET53
Definition dget53.f:126