LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlag2.f
Go to the documentation of this file.
1*> \brief \b DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLAG2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlag2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlag2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlag2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
22* WR2, WI )
23*
24* .. Scalar Arguments ..
25* INTEGER LDA, LDB
26* DOUBLE PRECISION SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), B( LDB, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
39*> problem A - w B, with scaling as necessary to avoid over-/underflow.
40*>
41*> The scaling factor "s" results in a modified eigenvalue equation
42*>
43*> s A - w B
44*>
45*> where s is a non-negative scaling factor chosen so that w, w B,
46*> and s A do not overflow and, if possible, do not underflow, either.
47*> \endverbatim
48*
49* Arguments:
50* ==========
51*
52*> \param[in] A
53*> \verbatim
54*> A is DOUBLE PRECISION array, dimension (LDA, 2)
55*> On entry, the 2 x 2 matrix A. It is assumed that its 1-norm
56*> is less than 1/SAFMIN. Entries less than
57*> sqrt(SAFMIN)*norm(A) are subject to being treated as zero.
58*> \endverbatim
59*>
60*> \param[in] LDA
61*> \verbatim
62*> LDA is INTEGER
63*> The leading dimension of the array A. LDA >= 2.
64*> \endverbatim
65*>
66*> \param[in] B
67*> \verbatim
68*> B is DOUBLE PRECISION array, dimension (LDB, 2)
69*> On entry, the 2 x 2 upper triangular matrix B. It is
70*> assumed that the one-norm of B is less than 1/SAFMIN. The
71*> diagonals should be at least sqrt(SAFMIN) times the largest
72*> element of B (in absolute value); if a diagonal is smaller
73*> than that, then +/- sqrt(SAFMIN) will be used instead of
74*> that diagonal.
75*> \endverbatim
76*>
77*> \param[in] LDB
78*> \verbatim
79*> LDB is INTEGER
80*> The leading dimension of the array B. LDB >= 2.
81*> \endverbatim
82*>
83*> \param[in] SAFMIN
84*> \verbatim
85*> SAFMIN is DOUBLE PRECISION
86*> The smallest positive number s.t. 1/SAFMIN does not
87*> overflow. (This should always be DLAMCH('S') -- it is an
88*> argument in order to avoid having to call DLAMCH frequently.)
89*> \endverbatim
90*>
91*> \param[out] SCALE1
92*> \verbatim
93*> SCALE1 is DOUBLE PRECISION
94*> A scaling factor used to avoid over-/underflow in the
95*> eigenvalue equation which defines the first eigenvalue. If
96*> the eigenvalues are complex, then the eigenvalues are
97*> ( WR1 +/- WI i ) / SCALE1 (which may lie outside the
98*> exponent range of the machine), SCALE1=SCALE2, and SCALE1
99*> will always be positive. If the eigenvalues are real, then
100*> the first (real) eigenvalue is WR1 / SCALE1 , but this may
101*> overflow or underflow, and in fact, SCALE1 may be zero or
102*> less than the underflow threshold if the exact eigenvalue
103*> is sufficiently large.
104*> \endverbatim
105*>
106*> \param[out] SCALE2
107*> \verbatim
108*> SCALE2 is DOUBLE PRECISION
109*> A scaling factor used to avoid over-/underflow in the
110*> eigenvalue equation which defines the second eigenvalue. If
111*> the eigenvalues are complex, then SCALE2=SCALE1. If the
112*> eigenvalues are real, then the second (real) eigenvalue is
113*> WR2 / SCALE2 , but this may overflow or underflow, and in
114*> fact, SCALE2 may be zero or less than the underflow
115*> threshold if the exact eigenvalue is sufficiently large.
116*> \endverbatim
117*>
118*> \param[out] WR1
119*> \verbatim
120*> WR1 is DOUBLE PRECISION
121*> If the eigenvalue is real, then WR1 is SCALE1 times the
122*> eigenvalue closest to the (2,2) element of A B**(-1). If the
123*> eigenvalue is complex, then WR1=WR2 is SCALE1 times the real
124*> part of the eigenvalues.
125*> \endverbatim
126*>
127*> \param[out] WR2
128*> \verbatim
129*> WR2 is DOUBLE PRECISION
130*> If the eigenvalue is real, then WR2 is SCALE2 times the
131*> other eigenvalue. If the eigenvalue is complex, then
132*> WR1=WR2 is SCALE1 times the real part of the eigenvalues.
133*> \endverbatim
134*>
135*> \param[out] WI
136*> \verbatim
137*> WI is DOUBLE PRECISION
138*> If the eigenvalue is real, then WI is zero. If the
139*> eigenvalue is complex, then WI is SCALE1 times the imaginary
140*> part of the eigenvalues. WI will always be non-negative.
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup lag2
152*
153* =====================================================================
154 SUBROUTINE dlag2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
155 $ WR2, WI )
156*
157* -- LAPACK auxiliary routine --
158* -- LAPACK is a software package provided by Univ. of Tennessee, --
159* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160*
161* .. Scalar Arguments ..
162 INTEGER LDA, LDB
163 DOUBLE PRECISION SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
164* ..
165* .. Array Arguments ..
166 DOUBLE PRECISION A( LDA, * ), B( LDB, * )
167* ..
168*
169* =====================================================================
170*
171* .. Parameters ..
172 DOUBLE PRECISION ZERO, ONE, TWO
173 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
174 DOUBLE PRECISION HALF
175 parameter( half = one / two )
176 DOUBLE PRECISION FUZZY1
177 parameter( fuzzy1 = one+1.0d-5 )
178* ..
179* .. Local Scalars ..
180 DOUBLE PRECISION A11, A12, A21, A22, ABI22, ANORM, AS11, AS12,
181 $ as22, ascale, b11, b12, b22, binv11, binv22,
182 $ bmin, bnorm, bscale, bsize, c1, c2, c3, c4, c5,
183 $ diff, discr, pp, qq, r, rtmax, rtmin, s1, s2,
184 $ safmax, shift, ss, sum, wabs, wbig, wdet,
185 $ wscale, wsize, wsmall
186* ..
187* .. Intrinsic Functions ..
188 INTRINSIC abs, max, min, sign, sqrt
189* ..
190* .. Executable Statements ..
191*
192 rtmin = sqrt( safmin )
193 rtmax = one / rtmin
194 safmax = one / safmin
195*
196* Scale A
197*
198 anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
199 $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
200 ascale = one / anorm
201 a11 = ascale*a( 1, 1 )
202 a21 = ascale*a( 2, 1 )
203 a12 = ascale*a( 1, 2 )
204 a22 = ascale*a( 2, 2 )
205*
206* Perturb B if necessary to insure non-singularity
207*
208 b11 = b( 1, 1 )
209 b12 = b( 1, 2 )
210 b22 = b( 2, 2 )
211 bmin = rtmin*max( abs( b11 ), abs( b12 ), abs( b22 ), rtmin )
212 IF( abs( b11 ).LT.bmin )
213 $ b11 = sign( bmin, b11 )
214 IF( abs( b22 ).LT.bmin )
215 $ b22 = sign( bmin, b22 )
216*
217* Scale B
218*
219 bnorm = max( abs( b11 ), abs( b12 )+abs( b22 ), safmin )
220 bsize = max( abs( b11 ), abs( b22 ) )
221 bscale = one / bsize
222 b11 = b11*bscale
223 b12 = b12*bscale
224 b22 = b22*bscale
225*
226* Compute larger eigenvalue by method described by C. van Loan
227*
228* ( AS is A shifted by -SHIFT*B )
229*
230 binv11 = one / b11
231 binv22 = one / b22
232 s1 = a11*binv11
233 s2 = a22*binv22
234 IF( abs( s1 ).LE.abs( s2 ) ) THEN
235 as12 = a12 - s1*b12
236 as22 = a22 - s1*b22
237 ss = a21*( binv11*binv22 )
238 abi22 = as22*binv22 - ss*b12
239 pp = half*abi22
240 shift = s1
241 ELSE
242 as12 = a12 - s2*b12
243 as11 = a11 - s2*b11
244 ss = a21*( binv11*binv22 )
245 abi22 = -ss*b12
246 pp = half*( as11*binv11+abi22 )
247 shift = s2
248 END IF
249 qq = ss*as12
250 IF( abs( pp*rtmin ).GE.one ) THEN
251 discr = ( rtmin*pp )**2 + qq*safmin
252 r = sqrt( abs( discr ) )*rtmax
253 ELSE
254 IF( pp**2+abs( qq ).LE.safmin ) THEN
255 discr = ( rtmax*pp )**2 + qq*safmax
256 r = sqrt( abs( discr ) )*rtmin
257 ELSE
258 discr = pp**2 + qq
259 r = sqrt( abs( discr ) )
260 END IF
261 END IF
262*
263* Note: the test of R in the following IF is to cover the case when
264* DISCR is small and negative and is flushed to zero during
265* the calculation of R. On machines which have a consistent
266* flush-to-zero threshold and handle numbers above that
267* threshold correctly, it would not be necessary.
268*
269 IF( discr.GE.zero .OR. r.EQ.zero ) THEN
270 sum = pp + sign( r, pp )
271 diff = pp - sign( r, pp )
272 wbig = shift + sum
273*
274* Compute smaller eigenvalue
275*
276 wsmall = shift + diff
277 IF( half*abs( wbig ).GT.max( abs( wsmall ), safmin ) ) THEN
278 wdet = ( a11*a22-a12*a21 )*( binv11*binv22 )
279 wsmall = wdet / wbig
280 END IF
281*
282* Choose (real) eigenvalue closest to 2,2 element of A*B**(-1)
283* for WR1.
284*
285 IF( pp.GT.abi22 ) THEN
286 wr1 = min( wbig, wsmall )
287 wr2 = max( wbig, wsmall )
288 ELSE
289 wr1 = max( wbig, wsmall )
290 wr2 = min( wbig, wsmall )
291 END IF
292 wi = zero
293 ELSE
294*
295* Complex eigenvalues
296*
297 wr1 = shift + pp
298 wr2 = wr1
299 wi = r
300 END IF
301*
302* Further scaling to avoid underflow and overflow in computing
303* SCALE1 and overflow in computing w*B.
304*
305* This scale factor (WSCALE) is bounded from above using C1 and C2,
306* and from below using C3 and C4.
307* C1 implements the condition s A must never overflow.
308* C2 implements the condition w B must never overflow.
309* C3, with C2,
310* implement the condition that s A - w B must never overflow.
311* C4 implements the condition s should not underflow.
312* C5 implements the condition max(s,|w|) should be at least 2.
313*
314 c1 = bsize*( safmin*max( one, ascale ) )
315 c2 = safmin*max( one, bnorm )
316 c3 = bsize*safmin
317 IF( ascale.LE.one .AND. bsize.LE.one ) THEN
318 c4 = min( one, ( ascale / safmin )*bsize )
319 ELSE
320 c4 = one
321 END IF
322 IF( ascale.LE.one .OR. bsize.LE.one ) THEN
323 c5 = min( one, ascale*bsize )
324 ELSE
325 c5 = one
326 END IF
327*
328* Scale first eigenvalue
329*
330 wabs = abs( wr1 ) + abs( wi )
331 wsize = max( safmin, c1, fuzzy1*( wabs*c2+c3 ),
332 $ min( c4, half*max( wabs, c5 ) ) )
333 IF( wsize.NE.one ) THEN
334 wscale = one / wsize
335 IF( wsize.GT.one ) THEN
336 scale1 = ( max( ascale, bsize )*wscale )*
337 $ min( ascale, bsize )
338 ELSE
339 scale1 = ( min( ascale, bsize )*wscale )*
340 $ max( ascale, bsize )
341 END IF
342 wr1 = wr1*wscale
343 IF( wi.NE.zero ) THEN
344 wi = wi*wscale
345 wr2 = wr1
346 scale2 = scale1
347 END IF
348 ELSE
349 scale1 = ascale*bsize
350 scale2 = scale1
351 END IF
352*
353* Scale second eigenvalue (if real)
354*
355 IF( wi.EQ.zero ) THEN
356 wsize = max( safmin, c1, fuzzy1*( abs( wr2 )*c2+c3 ),
357 $ min( c4, half*max( abs( wr2 ), c5 ) ) )
358 IF( wsize.NE.one ) THEN
359 wscale = one / wsize
360 IF( wsize.GT.one ) THEN
361 scale2 = ( max( ascale, bsize )*wscale )*
362 $ min( ascale, bsize )
363 ELSE
364 scale2 = ( min( ascale, bsize )*wscale )*
365 $ max( ascale, bsize )
366 END IF
367 wr2 = wr2*wscale
368 ELSE
369 scale2 = ascale*bsize
370 END IF
371 END IF
372*
373* End of DLAG2
374*
375 RETURN
376 END
dlag2
subroutine dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:156