LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
June 2016

Definition at line 158 of file dlag2.f.

158 *
159 * -- LAPACK auxiliary routine (version 3.6.1) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * June 2016
163 *
164 * .. Scalar Arguments ..
165  INTEGER lda, ldb
166  DOUBLE PRECISION safmin, scale1, scale2, wi, wr1, wr2
167 * ..
168 * .. Array Arguments ..
169  DOUBLE PRECISION a( lda, * ), b( ldb, * )
170 * ..
171 *
172 * =====================================================================
173 *
174 * .. Parameters ..
175  DOUBLE PRECISION zero, one, two
176  parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
177  DOUBLE PRECISION half
178  parameter ( half = one / two )
179  DOUBLE PRECISION fuzzy1
180  parameter ( fuzzy1 = one+1.0d-5 )
181 * ..
182 * .. Local Scalars ..
183  DOUBLE PRECISION a11, a12, a21, a22, abi22, anorm, as11, as12,
184  $ as22, ascale, b11, b12, b22, binv11, binv22,
185  $ bmin, bnorm, bscale, bsize, c1, c2, c3, c4, c5,
186  $ diff, discr, pp, qq, r, rtmax, rtmin, s1, s2,
187  $ safmax, shift, ss, sum, wabs, wbig, wdet,
188  $ wscale, wsize, wsmall
189 * ..
190 * .. Intrinsic Functions ..
191  INTRINSIC abs, max, min, sign, sqrt
192 * ..
193 * .. Executable Statements ..
194 *
195  rtmin = sqrt( safmin )
196  rtmax = one / rtmin
197  safmax = one / safmin
198 *
199 * Scale A
200 *
201  anorm = max( abs( a( 1, 1 ) )+abs( a( 2, 1 ) ),
202  $ abs( a( 1, 2 ) )+abs( a( 2, 2 ) ), safmin )
203  ascale = one / anorm
204  a11 = ascale*a( 1, 1 )
205  a21 = ascale*a( 2, 1 )
206  a12 = ascale*a( 1, 2 )
207  a22 = ascale*a( 2, 2 )
208 *
209 * Perturb B if necessary to insure non-singularity
210 *
211  b11 = b( 1, 1 )
212  b12 = b( 1, 2 )
213  b22 = b( 2, 2 )
214  bmin = rtmin*max( abs( b11 ), abs( b12 ), abs( b22 ), rtmin )
215  IF( abs( b11 ).LT.bmin )
216  $ b11 = sign( bmin, b11 )
217  IF( abs( b22 ).LT.bmin )
218  $ b22 = sign( bmin, b22 )
219 *
220 * Scale B
221 *
222  bnorm = max( abs( b11 ), abs( b12 )+abs( b22 ), safmin )
223  bsize = max( abs( b11 ), abs( b22 ) )
224  bscale = one / bsize
225  b11 = b11*bscale
226  b12 = b12*bscale
227  b22 = b22*bscale
228 *
229 * Compute larger eigenvalue by method described by C. van Loan
230 *
231 * ( AS is A shifted by -SHIFT*B )
232 *
233  binv11 = one / b11
234  binv22 = one / b22
235  s1 = a11*binv11
236  s2 = a22*binv22
237  IF( abs( s1 ).LE.abs( s2 ) ) THEN
238  as12 = a12 - s1*b12
239  as22 = a22 - s1*b22
240  ss = a21*( binv11*binv22 )
241  abi22 = as22*binv22 - ss*b12
242  pp = half*abi22
243  shift = s1
244  ELSE
245  as12 = a12 - s2*b12
246  as11 = a11 - s2*b11
247  ss = a21*( binv11*binv22 )
248  abi22 = -ss*b12
249  pp = half*( as11*binv11+abi22 )
250  shift = s2
251  END IF
252  qq = ss*as12
253  IF( abs( pp*rtmin ).GE.one ) THEN
254  discr = ( rtmin*pp )**2 + qq*safmin
255  r = sqrt( abs( discr ) )*rtmax
256  ELSE
257  IF( pp**2+abs( qq ).LE.safmin ) THEN
258  discr = ( rtmax*pp )**2 + qq*safmax
259  r = sqrt( abs( discr ) )*rtmin
260  ELSE
261  discr = pp**2 + qq
262  r = sqrt( abs( discr ) )
263  END IF
264  END IF
265 *
266 * Note: the test of R in the following IF is to cover the case when
267 * DISCR is small and negative and is flushed to zero during
268 * the calculation of R. On machines which have a consistent
269 * flush-to-zero threshold and handle numbers above that
270 * threshold correctly, it would not be necessary.
271 *
272  IF( discr.GE.zero .OR. r.EQ.zero ) THEN
273  sum = pp + sign( r, pp )
274  diff = pp - sign( r, pp )
275  wbig = shift + sum
276 *
277 * Compute smaller eigenvalue
278 *
279  wsmall = shift + diff
280  IF( half*abs( wbig ).GT.max( abs( wsmall ), safmin ) ) THEN
281  wdet = ( a11*a22-a12*a21 )*( binv11*binv22 )
282  wsmall = wdet / wbig
283  END IF
284 *
285 * Choose (real) eigenvalue closest to 2,2 element of A*B**(-1)
286 * for WR1.
287 *
288  IF( pp.GT.abi22 ) THEN
289  wr1 = min( wbig, wsmall )
290  wr2 = max( wbig, wsmall )
291  ELSE
292  wr1 = max( wbig, wsmall )
293  wr2 = min( wbig, wsmall )
294  END IF
295  wi = zero
296  ELSE
297 *
298 * Complex eigenvalues
299 *
300  wr1 = shift + pp
301  wr2 = wr1
302  wi = r
303  END IF
304 *
305 * Further scaling to avoid underflow and overflow in computing
306 * SCALE1 and overflow in computing w*B.
307 *
308 * This scale factor (WSCALE) is bounded from above using C1 and C2,
309 * and from below using C3 and C4.
310 * C1 implements the condition s A must never overflow.
311 * C2 implements the condition w B must never overflow.
312 * C3, with C2,
313 * implement the condition that s A - w B must never overflow.
314 * C4 implements the condition s should not underflow.
315 * C5 implements the condition max(s,|w|) should be at least 2.
316 *
317  c1 = bsize*( safmin*max( one, ascale ) )
318  c2 = safmin*max( one, bnorm )
319  c3 = bsize*safmin
320  IF( ascale.LE.one .AND. bsize.LE.one ) THEN
321  c4 = min( one, ( ascale / safmin )*bsize )
322  ELSE
323  c4 = one
324  END IF
325  IF( ascale.LE.one .OR. bsize.LE.one ) THEN
326  c5 = min( one, ascale*bsize )
327  ELSE
328  c5 = one
329  END IF
330 *
331 * Scale first eigenvalue
332 *
333  wabs = abs( wr1 ) + abs( wi )
334  wsize = max( safmin, c1, fuzzy1*( wabs*c2+c3 ),
335  $ min( c4, half*max( wabs, c5 ) ) )
336  IF( wsize.NE.one ) THEN
337  wscale = one / wsize
338  IF( wsize.GT.one ) THEN
339  scale1 = ( max( ascale, bsize )*wscale )*
340  $ min( ascale, bsize )
341  ELSE
342  scale1 = ( min( ascale, bsize )*wscale )*
343  $ max( ascale, bsize )
344  END IF
345  wr1 = wr1*wscale
346  IF( wi.NE.zero ) THEN
347  wi = wi*wscale
348  wr2 = wr1
349  scale2 = scale1
350  END IF
351  ELSE
352  scale1 = ascale*bsize
353  scale2 = scale1
354  END IF
355 *
356 * Scale second eigenvalue (if real)
357 *
358  IF( wi.EQ.zero ) THEN
359  wsize = max( safmin, c1, fuzzy1*( abs( wr2 )*c2+c3 ),
360  $ min( c4, half*max( abs( wr2 ), c5 ) ) )
361  IF( wsize.NE.one ) THEN
362  wscale = one / wsize
363  IF( wsize.GT.one ) THEN
364  scale2 = ( max( ascale, bsize )*wscale )*
365  $ min( ascale, bsize )
366  ELSE
367  scale2 = ( min( ascale, bsize )*wscale )*
368  $ max( ascale, bsize )
369  END IF
370  wr2 = wr2*wscale
371  ELSE
372  scale2 = ascale*bsize
373  END IF
374  END IF
375 *
376 * End of DLAG2
377 *
378  RETURN

Here is the caller graph for this function: