LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 *> \date November 2011
123 *
124 *> \ingroup double_eig
125 *
126 * =====================================================================
127  SUBROUTINE dget53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO )
128 *
129 * -- LAPACK test routine (version 3.4.0) --
130 * -- LAPACK is a software package provided by Univ. of Tennessee, --
131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
132 * November 2011
133 *
134 * .. Scalar Arguments ..
135  INTEGER info, lda, ldb
136  DOUBLE PRECISION result, scale, wi, wr
137 * ..
138 * .. Array Arguments ..
139  DOUBLE PRECISION a( lda, * ), b( ldb, * )
140 * ..
141 *
142 * =====================================================================
143 *
144 * .. Parameters ..
145  DOUBLE PRECISION zero, one
146  parameter( zero = 0.0d0, one = 1.0d0 )
147 * ..
148 * .. Local Scalars ..
149  DOUBLE PRECISION absw, anorm, bnorm, ci11, ci12, ci22, cnorm,
150  $ cr11, cr12, cr21, cr22, cscale, deti, detr, s1,
151  $ safmin, scales, sigmin, temp, ulp, wis, wrs
152 * ..
153 * .. External Functions ..
154  DOUBLE PRECISION dlamch
155  EXTERNAL dlamch
156 * ..
157 * .. Intrinsic Functions ..
158  INTRINSIC abs, max, sqrt
159 * ..
160 * .. Executable Statements ..
161 *
162 * Initialize
163 *
164  info = 0
165  result = zero
166  scales = scale
167  wrs = wr
168  wis = wi
169 *
170 * Machine constants and norms
171 *
172  safmin = dlamch( 'Safe minimum' )
173  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
174  absw = abs( wrs ) + abs( wis )
175  anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
176  $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
177  bnorm = max( abs( b( 1, 1 ) ), abs( b( 1, 2 ) )+abs( b( 2, 2 ) ),
178  $ safmin )
179 *
180 * Check for possible overflow.
181 *
182  temp = ( safmin*bnorm )*absw + ( safmin*anorm )*scales
183  IF( temp.GE.one ) THEN
184 *
185 * Scale down to avoid overflow
186 *
187  info = 1
188  temp = one / temp
189  scales = scales*temp
190  wrs = wrs*temp
191  wis = wis*temp
192  absw = abs( wrs ) + abs( wis )
193  END IF
194  s1 = max( ulp*max( scales*anorm, absw*bnorm ),
195  $ safmin*max( scales, absw ) )
196 *
197 * Check for W and SCALE essentially zero.
198 *
199  IF( s1.LT.safmin ) THEN
200  info = 2
201  IF( scales.LT.safmin .AND. absw.LT.safmin ) THEN
202  info = 3
203  result = one / ulp
204  return
205  END IF
206 *
207 * Scale up to avoid underflow
208 *
209  temp = one / max( scales*anorm+absw*bnorm, safmin )
210  scales = scales*temp
211  wrs = wrs*temp
212  wis = wis*temp
213  absw = abs( wrs ) + abs( wis )
214  s1 = max( ulp*max( scales*anorm, absw*bnorm ),
215  $ safmin*max( scales, absw ) )
216  IF( s1.LT.safmin ) THEN
217  info = 3
218  result = one / ulp
219  return
220  END IF
221  END IF
222 *
223 * Compute C = s A - w B
224 *
225  cr11 = scales*a( 1, 1 ) - wrs*b( 1, 1 )
226  ci11 = -wis*b( 1, 1 )
227  cr21 = scales*a( 2, 1 )
228  cr12 = scales*a( 1, 2 ) - wrs*b( 1, 2 )
229  ci12 = -wis*b( 1, 2 )
230  cr22 = scales*a( 2, 2 ) - wrs*b( 2, 2 )
231  ci22 = -wis*b( 2, 2 )
232 *
233 * Compute the smallest singular value of s A - w B:
234 *
235 * |det( s A - w B )|
236 * sigma_min = ------------------
237 * norm( s A - w B )
238 *
239  cnorm = max( abs( cr11 )+abs( ci11 )+abs( cr21 ),
240  $ abs( cr12 )+abs( ci12 )+abs( cr22 )+abs( ci22 ), safmin )
241  cscale = one / sqrt( cnorm )
242  detr = ( cscale*cr11 )*( cscale*cr22 ) -
243  $ ( cscale*ci11 )*( cscale*ci22 ) -
244  $ ( cscale*cr12 )*( cscale*cr21 )
245  deti = ( cscale*cr11 )*( cscale*ci22 ) +
246  $ ( cscale*ci11 )*( cscale*cr22 ) -
247  $ ( cscale*ci12 )*( cscale*cr21 )
248  sigmin = abs( detr ) + abs( deti )
249  result = sigmin / s1
250  return
251 *
252 * End of DGET53
253 *
254  END