LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlag2.f
Go to the documentation of this file.
1 *> \brief \b DLAG2 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 DLAG2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlag2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlag2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlag2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
22 * WR2, WI )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER LDA, LDB
26 * DOUBLE PRECISION SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION A( LDA, * ), B( LDB, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> DLAG2 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
86 *> The smallest positive number s.t. 1/SAFMIN does not
87 *> overflow. (This should always be DLAMCH('S') -- it is an
88 *> argument in order to avoid having to call DLAMCH frequently.)
89 *> \endverbatim
90 *>
91 *> \param[out] SCALE1
92 *> \verbatim
93 *> SCALE1 is DOUBLE PRECISION
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 threshhold if the exact eigenvalue
103 *> is sufficiently large.
104 *> \endverbatim
105 *>
106 *> \param[out] SCALE2
107 *> \verbatim
108 *> SCALE2 is DOUBLE PRECISION
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 *> threshhold if the exact eigenvalue is sufficiently large.
116 *> \endverbatim
117 *>
118 *> \param[out] WR1
119 *> \verbatim
120 *> WR1 is DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 *> \date September 2012
152 *
153 *> \ingroup doubleOTHERauxiliary
154 *
155 * =====================================================================
156  SUBROUTINE dlag2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1,
157  $ wr2, wi )
158 *
159 * -- LAPACK auxiliary routine (version 3.4.2) --
160 * -- LAPACK is a software package provided by Univ. of Tennessee, --
161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
162 * September 2012
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 threshhold and handle numbers above that
270 * threshhold 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
379  END