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