LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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.

Download DLAG2 + dependencies [TGZ] [ZIP] [TXT]

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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: