1: *> \brief \b DLAQZ1
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAQZ1 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz1.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz1.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz1.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAQZ1( A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2,
22: * $ V )
23: * IMPLICIT NONE
24: *
25: * Arguments
26: * INTEGER, INTENT( IN ) :: LDA, LDB
27: * DOUBLE PRECISION, INTENT( IN ) :: A( LDA, * ), B( LDB, * ), SR1,
28: * $ SR2, SI, BETA1, BETA2
29: * DOUBLE PRECISION, INTENT( OUT ) :: V( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> Given a 3-by-3 matrix pencil (A,B), DLAQZ1 sets v to a
39: *> scalar multiple of the first column of the product
40: *>
41: *> (*) K = (A - (beta2*sr2 - i*si)*B)*B^(-1)*(beta1*A - (sr2 + i*si2)*B)*B^(-1).
42: *>
43: *> It is assumed that either
44: *>
45: *> 1) sr1 = sr2
46: *> or
47: *> 2) si = 0.
48: *>
49: *> This is useful for starting double implicit shift bulges
50: *> in the QZ algorithm.
51: *> \endverbatim
52: *
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] A
58: *> \verbatim
59: *> A is DOUBLE PRECISION array, dimension (LDA,N)
60: *> The 3-by-3 matrix A in (*).
61: *> \endverbatim
62: *>
63: *> \param[in] LDA
64: *> \verbatim
65: *> LDA is INTEGER
66: *> The leading dimension of A as declared in
67: *> the calling procedure.
68: *> \endverbatim
69: *
70: *> \param[in] B
71: *> \verbatim
72: *> B is DOUBLE PRECISION array, dimension (LDB,N)
73: *> The 3-by-3 matrix B in (*).
74: *> \endverbatim
75: *>
76: *> \param[in] LDB
77: *> \verbatim
78: *> LDB is INTEGER
79: *> The leading dimension of B as declared in
80: *> the calling procedure.
81: *> \endverbatim
82: *>
83: *> \param[in] SR1
84: *> \verbatim
85: *> SR1 is DOUBLE PRECISION
86: *> \endverbatim
87: *>
88: *> \param[in] SR2
89: *> \verbatim
90: *> SR2 is DOUBLE PRECISION
91: *> \endverbatim
92: *>
93: *> \param[in] SI
94: *> \verbatim
95: *> SI is DOUBLE PRECISION
96: *> \endverbatim
97: *>
98: *> \param[in] BETA1
99: *> \verbatim
100: *> BETA1 is DOUBLE PRECISION
101: *> \endverbatim
102: *>
103: *> \param[in] BETA2
104: *> \verbatim
105: *> BETA2 is DOUBLE PRECISION
106: *> \endverbatim
107: *>
108: *> \param[out] V
109: *> \verbatim
110: *> V is DOUBLE PRECISION array, dimension (N)
111: *> A scalar multiple of the first column of the
112: *> matrix K in (*).
113: *> \endverbatim
114: *
115: * Authors:
116: * ========
117: *
118: *> \author Thijs Steel, KU Leuven
119: *
120: *> \date May 2020
121: *
122: *> \ingroup doubleGEcomputational
123: *>
124: * =====================================================================
125: SUBROUTINE DLAQZ1( A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2,
126: $ V )
127: IMPLICIT NONE
128: *
129: * Arguments
130: INTEGER, INTENT( IN ) :: LDA, LDB
131: DOUBLE PRECISION, INTENT( IN ) :: A( LDA, * ), B( LDB, * ), SR1,
132: $ SR2, SI, BETA1, BETA2
133: DOUBLE PRECISION, INTENT( OUT ) :: V( * )
134: *
135: * Parameters
136: DOUBLE PRECISION :: ZERO, ONE, HALF
137: PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
138: *
139: * Local scalars
140: DOUBLE PRECISION :: W( 2 ), SAFMIN, SAFMAX, SCALE1, SCALE2
141: *
142: * External Functions
143: DOUBLE PRECISION, EXTERNAL :: DLAMCH
144: LOGICAL, EXTERNAL :: DISNAN
145: *
146: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
147: SAFMAX = ONE/SAFMIN
148: *
149: * Calculate first shifted vector
150: *
151: W( 1 ) = BETA1*A( 1, 1 )-SR1*B( 1, 1 )
152: W( 2 ) = BETA1*A( 2, 1 )-SR1*B( 2, 1 )
153: SCALE1 = SQRT( ABS( W( 1 ) ) ) * SQRT( ABS( W( 2 ) ) )
154: IF( SCALE1 .GE. SAFMIN .AND. SCALE1 .LE. SAFMAX ) THEN
155: W( 1 ) = W( 1 )/SCALE1
156: W( 2 ) = W( 2 )/SCALE1
157: END IF
158: *
159: * Solve linear system
160: *
161: W( 2 ) = W( 2 )/B( 2, 2 )
162: W( 1 ) = ( W( 1 )-B( 1, 2 )*W( 2 ) )/B( 1, 1 )
163: SCALE2 = SQRT( ABS( W( 1 ) ) ) * SQRT( ABS( W( 2 ) ) )
164: IF( SCALE2 .GE. SAFMIN .AND. SCALE2 .LE. SAFMAX ) THEN
165: W( 1 ) = W( 1 )/SCALE2
166: W( 2 ) = W( 2 )/SCALE2
167: END IF
168: *
169: * Apply second shift
170: *
171: V( 1 ) = BETA2*( A( 1, 1 )*W( 1 )+A( 1, 2 )*W( 2 ) )-SR2*( B( 1,
172: $ 1 )*W( 1 )+B( 1, 2 )*W( 2 ) )
173: V( 2 ) = BETA2*( A( 2, 1 )*W( 1 )+A( 2, 2 )*W( 2 ) )-SR2*( B( 2,
174: $ 1 )*W( 1 )+B( 2, 2 )*W( 2 ) )
175: V( 3 ) = BETA2*( A( 3, 1 )*W( 1 )+A( 3, 2 )*W( 2 ) )-SR2*( B( 3,
176: $ 1 )*W( 1 )+B( 3, 2 )*W( 2 ) )
177: *
178: * Account for imaginary part
179: *
180: V( 1 ) = V( 1 )+SI*SI*B( 1, 1 )/SCALE1/SCALE2
181: *
182: * Check for overflow
183: *
184: IF( ABS( V( 1 ) ).GT.SAFMAX .OR. ABS( V( 2 ) ) .GT. SAFMAX .OR.
185: $ ABS( V( 3 ) ).GT.SAFMAX .OR. DISNAN( V( 1 ) ) .OR.
186: $ DISNAN( V( 2 ) ) .OR. DISNAN( V( 3 ) ) ) THEN
187: V( 1 ) = ZERO
188: V( 2 ) = ZERO
189: V( 3 ) = ZERO
190: END IF
191: *
192: * End of DLAQZ1
193: *
194: END SUBROUTINE
CVSweb interface <joel.bertrand@systella.fr>