LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dlag2()

 subroutine dlag2 ( double precision, dimension( lda, * ) a, integer lda, double precision, dimension( ldb, * ) b, integer ldb, double precision safmin, double precision scale1, double precision scale2, double precision wr1, double precision wr2, double precision wi )

DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow.

Purpose:
``` DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
problem  A - w B, with scaling as necessary to avoid over-/underflow.

The scaling factor "s" results in a modified eigenvalue equation

s A - w B

where  s  is a non-negative scaling factor chosen so that  w,  w B,
and  s A  do not overflow and, if possible, do not underflow, either.```
Parameters
 [in] A ``` A is DOUBLE PRECISION array, dimension (LDA, 2) On entry, the 2 x 2 matrix A. It is assumed that its 1-norm is less than 1/SAFMIN. Entries less than sqrt(SAFMIN)*norm(A) are subject to being treated as zero.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= 2.``` [in] B ``` B is DOUBLE PRECISION array, dimension (LDB, 2) On entry, the 2 x 2 upper triangular matrix B. It is assumed that the one-norm of B is less than 1/SAFMIN. The diagonals should be at least sqrt(SAFMIN) times the largest element of B (in absolute value); if a diagonal is smaller than that, then +/- sqrt(SAFMIN) will be used instead of that diagonal.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= 2.``` [in] SAFMIN ``` SAFMIN is DOUBLE PRECISION The smallest positive number s.t. 1/SAFMIN does not overflow. (This should always be DLAMCH('S') -- it is an argument in order to avoid having to call DLAMCH frequently.)``` [out] SCALE1 ``` SCALE1 is DOUBLE PRECISION A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the first eigenvalue. If the eigenvalues are complex, then the eigenvalues are ( WR1 +/- WI i ) / SCALE1 (which may lie outside the exponent range of the machine), SCALE1=SCALE2, and SCALE1 will always be positive. If the eigenvalues are real, then the first (real) eigenvalue is WR1 / SCALE1 , but this may overflow or underflow, and in fact, SCALE1 may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large.``` [out] SCALE2 ``` SCALE2 is DOUBLE PRECISION A scaling factor used to avoid over-/underflow in the eigenvalue equation which defines the second eigenvalue. If the eigenvalues are complex, then SCALE2=SCALE1. If the eigenvalues are real, then the second (real) eigenvalue is WR2 / SCALE2 , but this may overflow or underflow, and in fact, SCALE2 may be zero or less than the underflow threshold if the exact eigenvalue is sufficiently large.``` [out] WR1 ``` WR1 is DOUBLE PRECISION If the eigenvalue is real, then WR1 is SCALE1 times the eigenvalue closest to the (2,2) element of A B**(-1). If the eigenvalue is complex, then WR1=WR2 is SCALE1 times the real part of the eigenvalues.``` [out] WR2 ``` WR2 is DOUBLE PRECISION If the eigenvalue is real, then WR2 is SCALE2 times the other eigenvalue. If the eigenvalue is complex, then WR1=WR2 is SCALE1 times the real part of the eigenvalues.``` [out] WI ``` WI is DOUBLE PRECISION If the eigenvalue is real, then WI is zero. If the eigenvalue is complex, then WI is SCALE1 times the imaginary part of the eigenvalues. WI will always be non-negative.```

Definition at line 154 of file dlag2.f.

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
Here is the caller graph for this function: