LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slag2()

subroutine slag2 ( real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real  SAFMIN,
real  SCALE1,
real  SCALE2,
real  WR1,
real  WR2,
real  WI 
)

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

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

Purpose:
 SLAG2 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 REAL 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 REAL 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 REAL
          The smallest positive number s.t. 1/SAFMIN does not
          overflow.  (This should always be SLAMCH('S') -- it is an
          argument in order to avoid having to call SLAMCH frequently.)
[out]SCALE1
          SCALE1 is REAL
          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 REAL
          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 REAL
          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 REAL
          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 REAL
          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 slag2.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  REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
164 * ..
165 * .. Array Arguments ..
166  REAL A( LDA, * ), B( LDB, * )
167 * ..
168 *
169 * =====================================================================
170 *
171 * .. Parameters ..
172  REAL ZERO, ONE, TWO
173  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
174  REAL HALF
175  parameter( half = one / two )
176  REAL FUZZY1
177  parameter( fuzzy1 = one+1.0e-5 )
178 * ..
179 * .. Local Scalars ..
180  REAL 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 SLAG2
374 *
375  RETURN
Here is the caller graph for this function: