LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
slag2.f
Go to the documentation of this file.
1 *> \brief \b SLAG2 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 *> \htmlonly
9 *> Download SLAG2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slag2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slag2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slag2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
22 * WR2, WI )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER LDA, LDB
26 * REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
27 * ..
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), B( LDB, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> SLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue
39 *> problem A - w B, with scaling as necessary to avoid over-/underflow.
40 *>
41 *> The scaling factor "s" results in a modified eigenvalue equation
42 *>
43 *> s A - w B
44 *>
45 *> where s is a non-negative scaling factor chosen so that w, w B,
46 *> and s A do not overflow and, if possible, do not underflow, either.
47 *> \endverbatim
48 *
49 * Arguments:
50 * ==========
51 *
52 *> \param[in] A
53 *> \verbatim
54 *> A is REAL array, dimension (LDA, 2)
55 *> On entry, the 2 x 2 matrix A. It is assumed that its 1-norm
56 *> is less than 1/SAFMIN. Entries less than
57 *> sqrt(SAFMIN)*norm(A) are subject to being treated as zero.
58 *> \endverbatim
59 *>
60 *> \param[in] LDA
61 *> \verbatim
62 *> LDA is INTEGER
63 *> The leading dimension of the array A. LDA >= 2.
64 *> \endverbatim
65 *>
66 *> \param[in] B
67 *> \verbatim
68 *> B is REAL array, dimension (LDB, 2)
69 *> On entry, the 2 x 2 upper triangular matrix B. It is
70 *> assumed that the one-norm of B is less than 1/SAFMIN. The
71 *> diagonals should be at least sqrt(SAFMIN) times the largest
72 *> element of B (in absolute value); if a diagonal is smaller
73 *> than that, then +/- sqrt(SAFMIN) will be used instead of
74 *> that diagonal.
75 *> \endverbatim
76 *>
77 *> \param[in] LDB
78 *> \verbatim
79 *> LDB is INTEGER
80 *> The leading dimension of the array B. LDB >= 2.
81 *> \endverbatim
82 *>
83 *> \param[in] SAFMIN
84 *> \verbatim
85 *> SAFMIN is REAL
86 *> The smallest positive number s.t. 1/SAFMIN does not
87 *> overflow. (This should always be SLAMCH('S') -- it is an
88 *> argument in order to avoid having to call SLAMCH frequently.)
89 *> \endverbatim
90 *>
91 *> \param[out] SCALE1
92 *> \verbatim
93 *> SCALE1 is REAL
94 *> A scaling factor used to avoid over-/underflow in the
95 *> eigenvalue equation which defines the first eigenvalue. If
96 *> the eigenvalues are complex, then the eigenvalues are
97 *> ( WR1 +/- WI i ) / SCALE1 (which may lie outside the
98 *> exponent range of the machine), SCALE1=SCALE2, and SCALE1
99 *> will always be positive. If the eigenvalues are real, then
100 *> the first (real) eigenvalue is WR1 / SCALE1 , but this may
101 *> overflow or underflow, and in fact, SCALE1 may be zero or
102 *> less than the underflow threshold if the exact eigenvalue
103 *> is sufficiently large.
104 *> \endverbatim
105 *>
106 *> \param[out] SCALE2
107 *> \verbatim
108 *> SCALE2 is REAL
109 *> A scaling factor used to avoid over-/underflow in the
110 *> eigenvalue equation which defines the second eigenvalue. If
111 *> the eigenvalues are complex, then SCALE2=SCALE1. If the
112 *> eigenvalues are real, then the second (real) eigenvalue is
113 *> WR2 / SCALE2 , but this may overflow or underflow, and in
114 *> fact, SCALE2 may be zero or less than the underflow
115 *> threshold if the exact eigenvalue is sufficiently large.
116 *> \endverbatim
117 *>
118 *> \param[out] WR1
119 *> \verbatim
120 *> WR1 is REAL
121 *> If the eigenvalue is real, then WR1 is SCALE1 times the
122 *> eigenvalue closest to the (2,2) element of A B**(-1). If the
123 *> eigenvalue is complex, then WR1=WR2 is SCALE1 times the real
124 *> part of the eigenvalues.
125 *> \endverbatim
126 *>
127 *> \param[out] WR2
128 *> \verbatim
129 *> WR2 is REAL
130 *> If the eigenvalue is real, then WR2 is SCALE2 times the
131 *> other eigenvalue. If the eigenvalue is complex, then
132 *> WR1=WR2 is SCALE1 times the real part of the eigenvalues.
133 *> \endverbatim
134 *>
135 *> \param[out] WI
136 *> \verbatim
137 *> WI is REAL
138 *> If the eigenvalue is real, then WI is zero. If the
139 *> eigenvalue is complex, then WI is SCALE1 times the imaginary
140 *> part of the eigenvalues. WI will always be non-negative.
141 *> \endverbatim
142 *
143 * Authors:
144 * ========
145 *
146 *> \author Univ. of Tennessee
147 *> \author Univ. of California Berkeley
148 *> \author Univ. of Colorado Denver
149 *> \author NAG Ltd.
150 *
151 *> \ingroup realOTHERauxiliary
152 *
153 * =====================================================================
154  SUBROUTINE slag2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
155  $ WR2, WI )
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
376  END
subroutine slag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: slag2.f:156