1: SUBROUTINE ZLAQR1( N, H, LDH, S1, S2, V )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
5: * November 2006
6: *
7: * .. Scalar Arguments ..
8: COMPLEX*16 S1, S2
9: INTEGER LDH, N
10: * ..
11: * .. Array Arguments ..
12: COMPLEX*16 H( LDH, * ), V( * )
13: * ..
14: *
15: * Given a 2-by-2 or 3-by-3 matrix H, ZLAQR1 sets v to a
16: * scalar multiple of the first column of the product
17: *
18: * (*) K = (H - s1*I)*(H - s2*I)
19: *
20: * scaling to avoid overflows and most underflows.
21: *
22: * This is useful for starting double implicit shift bulges
23: * in the QR algorithm.
24: *
25: *
26: * N (input) integer
27: * Order of the matrix H. N must be either 2 or 3.
28: *
29: * H (input) COMPLEX*16 array of dimension (LDH,N)
30: * The 2-by-2 or 3-by-3 matrix H in (*).
31: *
32: * LDH (input) integer
33: * The leading dimension of H as declared in
34: * the calling procedure. LDH.GE.N
35: *
36: * S1 (input) COMPLEX*16
37: * S2 S1 and S2 are the shifts defining K in (*) above.
38: *
39: * V (output) COMPLEX*16 array of dimension N
40: * A scalar multiple of the first column of the
41: * matrix K in (*).
42: *
43: * ================================================================
44: * Based on contributions by
45: * Karen Braman and Ralph Byers, Department of Mathematics,
46: * University of Kansas, USA
47: *
48: * ================================================================
49: *
50: * .. Parameters ..
51: COMPLEX*16 ZERO
52: PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ) )
53: DOUBLE PRECISION RZERO
54: PARAMETER ( RZERO = 0.0d0 )
55: * ..
56: * .. Local Scalars ..
57: COMPLEX*16 CDUM, H21S, H31S
58: DOUBLE PRECISION S
59: * ..
60: * .. Intrinsic Functions ..
61: INTRINSIC ABS, DBLE, DIMAG
62: * ..
63: * .. Statement Functions ..
64: DOUBLE PRECISION CABS1
65: * ..
66: * .. Statement Function definitions ..
67: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
68: * ..
69: * .. Executable Statements ..
70: IF( N.EQ.2 ) THEN
71: S = CABS1( H( 1, 1 )-S2 ) + CABS1( H( 2, 1 ) )
72: IF( S.EQ.RZERO ) THEN
73: V( 1 ) = ZERO
74: V( 2 ) = ZERO
75: ELSE
76: H21S = H( 2, 1 ) / S
77: V( 1 ) = H21S*H( 1, 2 ) + ( H( 1, 1 )-S1 )*
78: $ ( ( H( 1, 1 )-S2 ) / S )
79: V( 2 ) = H21S*( H( 1, 1 )+H( 2, 2 )-S1-S2 )
80: END IF
81: ELSE
82: S = CABS1( H( 1, 1 )-S2 ) + CABS1( H( 2, 1 ) ) +
83: $ CABS1( H( 3, 1 ) )
84: IF( S.EQ.ZERO ) THEN
85: V( 1 ) = ZERO
86: V( 2 ) = ZERO
87: V( 3 ) = ZERO
88: ELSE
89: H21S = H( 2, 1 ) / S
90: H31S = H( 3, 1 ) / S
91: V( 1 ) = ( H( 1, 1 )-S1 )*( ( H( 1, 1 )-S2 ) / S ) +
92: $ H( 1, 2 )*H21S + H( 1, 3 )*H31S
93: V( 2 ) = H21S*( H( 1, 1 )+H( 2, 2 )-S1-S2 ) + H( 2, 3 )*H31S
94: V( 3 ) = H31S*( H( 1, 1 )+H( 3, 3 )-S1-S2 ) + H21S*H( 3, 2 )
95: END IF
96: END IF
97: END
CVSweb interface <joel.bertrand@systella.fr>