GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
exchqz.f
Go to the documentation of this file.
1  SUBROUTINE exchqz(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL)
2  INTEGER NMAX, N, L, LS1, LS2
3  DOUBLE PRECISION A(nmax,n), B(nmax,n), Z(nmax,n), EPS
4  LOGICAL FAIL
5 c modified july 9, 1998 a.s.hodel@eng.auburn.edu:
6 c REAL changed to DOUBLE PRECISION
7 c calls to AMAX1 changed to call MAX instead.
8 c calls to SROT changed to DROT (both in BLAS)
9 c calls to giv changed to dlartg (LAPACK); required new variable tempr
10 C*
11 C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A
12 C* WITH CONSECUTIVE LS1XLS1 AND LS2XLS2 DIAGONAL BLOCKS (LS1,LS2.LE.2)
13 C* STARTING AT ROW/COLUMN L, EXCHQZ PRODUCES EQUIVALENCE TRANSFORMA-
14 C* TIONS QT AND ZT THAT EXCHANGE THE BLOCKS ALONG WITH THEIR GENERALIZED
15 C* EIGENVALUES. EXCHQZ REQUIRES THE SUBROUTINES DROT (BLAS) AND GIV.
16 C* THE PARAMETERS IN THE CALLING SEQUENCE ARE (STARRED PARAMETERS ARE
17 C* ALTERED BY THE SUBROUTINE):
18 C*
19 C* NMAX THE FIRST DIMENSION OF A, B AND Z
20 C* N THE ORDER OF A, B AND Z
21 C* *A,*B THE MATRIX PAIR WHOSE BLOCKS ARE TO BE INTERCHANGED
22 C* *Z UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN
23 C* TRANSFORMATION ZT.
24 C* L THE POSITION OF THE BLOCKS
25 C* LS1 THE SIZE OF THE FIRST BLOCK
26 C* LS2 THE SIZE OF THE SECOND BLOCK
27 C* EPS THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT
28 C* *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN,
29 C* TRUE OTHERWISE.
30 C*
31  INTEGER I, J, L1, L2, L3, LI, LJ, LL, IT1, IT2
32  DOUBLE PRECISION U(3,3), D, E, F, G, SA, SB, A11B11, A21B11,
33  * a12b22, b12b22,
34  * a22b22, ammbmm, anmbmm, amnbnn, bmnbnn, annbnn, tempr
35  LOGICAL ALTB
36  fail = .false.
37  l1 = l + 1
38  ll = ls1 + ls2
39  IF (ll.GT.2) go to 10
40 C*** INTERCHANGE 1X1 AND 1X1 BLOCKS VIA AN EQUIVALENCE
41 C*** TRANSFORMATION A:=Q*A*Z , B:=Q*B*Z
42 C*** WHERE Q AND Z ARE GIVENS ROTATIONS
43  f = max(abs(a(l1,l1)),abs(b(l1,l1)))
44  altb = .true.
45  IF (abs(a(l1,l1)).GE.f) altb = .false.
46  sa = a(l1,l1)/f
47  sb = b(l1,l1)/f
48  f = sa*b(l,l) - sb*a(l,l)
49 C* CONSTRUCT THE COLUMN TRANSFORMATION Z
50  g = sa*b(l,l1) - sb*a(l,l1)
51  CALL dlartg(f, g, d, e,tempr)
52  CALL drot(l1, a(1,l), 1, a(1,l1), 1, e, -d)
53  CALL drot(l1, b(1,l), 1, b(1,l1), 1, e, -d)
54  CALL drot(n, z(1,l), 1, z(1,l1), 1, e, -d)
55 C* CONSTRUCT THE ROW TRANSFORMATION Q
56  IF (altb) CALL dlartg(b(l,l), b(l1,l), d, e,tempr)
57  IF (.NOT.altb) CALL dlartg(a(l,l), a(l1,l), d, e,tempr)
58  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
59  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
60  a(l1,l) = 0.
61  b(l1,l) = 0.
62  RETURN
63 C*** INTERCHANGE 1X1 AND 2X2 BLOCKS VIA AN EQUIVALENCE
64 C*** TRANSFORMATION A:=Q2*Q1*A*Z1*Z2 , B:=Q2*Q1*B*Z1*Z2
65 C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION
66  10 l2 = l + 2
67  IF (ls1.EQ.2) go to 60
68  g = max(abs(a(l,l)),abs(b(l,l)))
69  altb = .true.
70  IF (abs(a(l,l)).LT.g) go to 20
71  altb = .false.
72  CALL dlartg(a(l1,l1), a(l2,l1), d, e,tempr)
73  CALL drot(n-l, a(l1,l1), nmax, a(l2,l1), nmax, d, e)
74  CALL drot(n-l, b(l1,l1), nmax, b(l2,l1), nmax, d, e)
75 C** EVALUATE THE PENCIL AT THE EIGENVALUE CORRESPONDING
76 C** TO THE 1X1 BLOCK
77  20 sa = a(l,l)/g
78  sb = b(l,l)/g
79  DO 40 j=1,2
80  lj = l + j
81  DO 30 i=1,3
82  li = l + i - 1
83  u(i,j) = sa*b(li,lj) - sb*a(li,lj)
84  30 CONTINUE
85  40 CONTINUE
86  CALL dlartg(u(3,1), u(3,2), d, e,tempr)
87  CALL drot(3, u(1,1), 1, u(1,2), 1, e, -d)
88 C* PERFORM THE ROW TRANSFORMATION Q1
89  CALL dlartg(u(1,1), u(2,1), d, e,tempr)
90  u(2,2) = -u(1,2)*e + u(2,2)*d
91  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
92  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
93 C* PERFORM THE COLUMN TRANSFORMATION Z1
94  IF (altb) CALL dlartg(b(l1,l), b(l1,l1), d, e,tempr)
95  IF (.NOT.altb) CALL dlartg(a(l1,l), a(l1,l1), d, e,tempr)
96  CALL drot(l2, a(1,l), 1, a(1,l1), 1, e, -d)
97  CALL drot(l2, b(1,l), 1, b(1,l1), 1, e, -d)
98  CALL drot(n, z(1,l), 1, z(1,l1), 1, e, -d)
99 C* PERFORM THE ROW TRANSFORMATION Q2
100  CALL dlartg(u(2,2), u(3,2), d, e,tempr)
101  CALL drot(n-l+1, a(l1,l), nmax, a(l2,l), nmax, d, e)
102  CALL drot(n-l+1, b(l1,l), nmax, b(l2,l), nmax, d, e)
103 C* PERFORM THE COLUMN TRANSFORMATION Z2
104  IF (altb) CALL dlartg(b(l2,l1), b(l2,l2), d, e,tempr)
105  IF (.NOT.altb) CALL dlartg(a(l2,l1), a(l2,l2), d, e,tempr)
106  CALL drot(l2, a(1,l1), 1, a(1,l2), 1, e, -d)
107  CALL drot(l2, b(1,l1), 1, b(1,l2), 1, e, -d)
108  CALL drot(n, z(1,l1), 1, z(1,l2), 1, e, -d)
109  IF (altb) go to 50
110  CALL dlartg(b(l,l), b(l1,l), d, e,tempr)
111  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
112  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
113 C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO
114  50 a(l2,l) = 0.
115  a(l2,l1) = 0.
116  b(l1,l) = 0.
117  b(l2,l) = 0.
118  b(l2,l1) = 0.
119  RETURN
120 C*** INTERCHANGE 2X2 AND 1X1 BLOCKS VIA AN EQUIVALENCE
121 C*** TRANSFORMATION A:=Q2*Q1*A*Z1*Z2 , B:=Q2*Q1*B*Z1*Z2
122 C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION
123  60 IF (ls2.EQ.2) go to 110
124  g = max(abs(a(l2,l2)),abs(b(l2,l2)))
125  altb = .true.
126  IF (abs(a(l2,l2)).LT.g) go to 70
127  altb = .false.
128  CALL dlartg(a(l,l), a(l1,l), d, e,tempr)
129  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
130  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
131 C** EVALUATE THE PENCIL AT THE EIGENVALUE CORRESPONDING
132 C** TO THE 1X1 BLOCK
133  70 sa = a(l2,l2)/g
134  sb = b(l2,l2)/g
135  DO 90 i=1,2
136  li = l + i - 1
137  DO 80 j=1,3
138  lj = l + j - 1
139  u(i,j) = sa*b(li,lj) - sb*a(li,lj)
140  80 CONTINUE
141  90 CONTINUE
142  CALL dlartg(u(1,1), u(2,1), d, e,tempr)
143  CALL drot(3, u(1,1), 3, u(2,1), 3, d, e)
144 C* PERFORM THE COLUMN TRANSFORMATION Z1
145  CALL dlartg(u(2,2), u(2,3), d, e,tempr)
146  u(1,2) = u(1,2)*e - u(1,3)*d
147  CALL drot(l2, a(1,l1), 1, a(1,l2), 1, e, -d)
148  CALL drot(l2, b(1,l1), 1, b(1,l2), 1, e, -d)
149  CALL drot(n, z(1,l1), 1, z(1,l2), 1, e, -d)
150 C* PERFORM THE ROW TRANSFORMATION Q1
151  IF (altb) CALL dlartg(b(l1,l1), b(l2,l1), d, e,tempr)
152  IF (.NOT.altb) CALL dlartg(a(l1,l1), a(l2,l1), d, e,tempr)
153  CALL drot(n-l+1, a(l1,l), nmax, a(l2,l), nmax, d, e)
154  CALL drot(n-l+1, b(l1,l), nmax, b(l2,l), nmax, d, e)
155 C* PERFORM THE COLUMN TRANSFORMATION Z2
156  CALL dlartg(u(1,1), u(1,2), d, e,tempr)
157  CALL drot(l2, a(1,l), 1, a(1,l1), 1, e, -d)
158  CALL drot(l2, b(1,l), 1, b(1,l1), 1, e, -d)
159  CALL drot(n, z(1,l), 1, z(1,l1), 1, e, -d)
160 C* PERFORM THE ROW TRANSFORMATION Q2
161  IF (altb) CALL dlartg(b(l,l), b(l1,l), d, e,tempr)
162  IF (.NOT.altb) CALL dlartg(a(l,l), a(l1,l), d, e,tempr)
163  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
164  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
165  IF (altb) go to 100
166  CALL dlartg(b(l1,l1), b(l2,l1), d, e,tempr)
167  CALL drot(n-l, a(l1,l1), nmax, a(l2,l1), nmax, d, e)
168  CALL drot(n-l, b(l1,l1), nmax, b(l2,l1), nmax, d, e)
169 C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO
170  100 a(l1,l) = 0.
171  a(l2,l) = 0.
172  b(l1,l) = 0.
173  b(l2,1) = 0.
174  b(l2,l1) = 0.
175  RETURN
176 C*** INTERCHANGE 2X2 AND 2X2 BLOCKS VIA A SEQUENCE OF
177 C*** QZ-STEPS REALIZED BY THE EQUIVALENCE TRANSFORMATIONS
178 C*** A:=Q5*Q4*Q3*Q2*Q1*A*Z1*Z2*Z3*Z4*Z5
179 C*** B:=Q5*Q4*Q3*Q2*Q1*B*Z1*Z2*Z3*Z4*Z5
180 C*** WHERE EACH QI AND ZI IS A GIVENS ROTATION
181  110 l3 = l + 3
182 C* COMPUTE IMPLICIT SHIFT
183  ammbmm = a(l,l)/b(l,l)
184  anmbmm = a(l1,l)/b(l,l)
185  amnbnn = a(l,l1)/b(l1,l1)
186  annbnn = a(l1,l1)/b(l1,l1)
187  bmnbnn = b(l,l1)/b(l1,l1)
188  DO 130 it1=1,3
189  u(1,1) = 1.
190  u(2,1) = 1.
191  u(3,1) = 1.
192  DO 120 it2=1,10
193 C* PERFORM ROW TRANSFORMATIONS Q1 AND Q2
194  CALL dlartg(u(2,1), u(3,1), d, e,tempr)
195  CALL drot(n-l+1, a(l1,l), nmax, a(l2,l), nmax, d, e)
196  CALL drot(n-l, b(l1,l1), nmax, b(l2,l1), nmax, d, e)
197  u(2,1) = d*u(2,1) + e*u(3,1)
198  CALL dlartg(u(1,1), u(2,1), d, e,tempr)
199  CALL drot(n-l+1, a(l,l), nmax, a(l1,l), nmax, d, e)
200  CALL drot(n-l+1, b(l,l), nmax, b(l1,l), nmax, d, e)
201 C* PERFORM COLUMN TRANSFORMATIONS Z1 AND Z2
202  CALL dlartg(b(l2,l1), b(l2,l2), d, e,tempr)
203  CALL drot(l3, a(1,l1), 1, a(1,l2), 1, e, -d)
204  CALL drot(l2, b(1,l1), 1, b(1,l2), 1, e, -d)
205  CALL drot(n, z(1,l1), 1, z(1,l2), 1, e, -d)
206  CALL dlartg(b(l1,l), b(l1,l1), d, e,tempr)
207  CALL drot(l3, a(1,l), 1, a(1,l1), 1, e, -d)
208  CALL drot(l1, b(1,l), 1, b(1,l1), 1, e, -d)
209  CALL drot(n, z(1,l), 1, z(1,l1), 1, e, -d)
210 C* PERFORM TRANSFORMATIONS Q3,Z3,Q4,Z4,Q5 AND Z5 IN
211 C* ORDER TO REDUCE THE PENCIL TO HESSENBERG FORM
212  CALL dlartg(a(l2,l), a(l3,l), d, e,tempr)
213  CALL drot(n-l+1, a(l2,l), nmax, a(l3,l), nmax, d, e)
214  CALL drot(n-l1, b(l2,l2), nmax, b(l3,l2), nmax, d, e)
215  CALL dlartg(b(l3,l2), b(l3,l3), d, e,tempr)
216  CALL drot(l3, a(1,l2), 1, a(1,l3), 1, e, -d)
217  CALL drot(l3, b(1,l2), 1, b(1,l3), 1, e, -d)
218  CALL drot(n, z(1,l2), 1, z(1,l3), 1, e, -d)
219  CALL dlartg(a(l1,l), a(l2,l), d, e,tempr)
220  CALL drot(n-l+1, a(l1,l), nmax, a(l2,l), nmax, d, e)
221  CALL drot(n-l, b(l1,l1), nmax, b(l2,l1), nmax, d, e)
222  CALL dlartg(b(l2,l1), b(l2,l2), d, e,tempr)
223  CALL drot(l3, a(1,l1), 1, a(1,l2), 1, e, -d)
224  CALL drot(l2, b(1,l1), 1, b(1,l2), 1, e, -d)
225  CALL drot(n, z(1,l1), 1, z(1,l2), 1, e, -d)
226  CALL dlartg(a(l2,l1), a(l3,l1), d, e,tempr)
227  CALL drot(n-l, a(l2,l1), nmax, a(l3,l1), nmax, d, e)
228  CALL drot(n-l1, b(l2,l2), nmax, b(l3,l2), nmax, d, e)
229  CALL dlartg(b(l3,l2), b(l3,l3), d, e,tempr)
230  CALL drot(l3, a(1,l2), 1, a(1,l3), 1, e, -d)
231  CALL drot(l3, b(1,l2), 1, b(1,l3), 1, e, -d)
232  CALL drot(n, z(1,l2), 1, z(1,l3), 1, e, -d)
233 C* TEST OF CONVERGENCE ON THE ELEMENT SEPARATING THE BLOCKS
234  IF (abs(a(l2,l1)).LE.eps) go to 140
235 C* COMPUTE A NEW SHIFT IN CASE OF NO CONVERGENCE
236  a11b11 = a(l,l)/b(l,l)
237  a12b22 = a(l,l1)/b(l1,l1)
238  a21b11 = a(l1,l)/b(l,l)
239  a22b22 = a(l1,l1)/b(l1,l1)
240  b12b22 = b(l,l1)/b(l1,l1)
241  u(1,1) = ((ammbmm-a11b11)*(annbnn-a11b11)-amnbnn*
242  * anmbmm+anmbmm*bmnbnn*a11b11)/a21b11 + a12b22 - a11b11*b12b22
243  u(2,1) = (a22b22-a11b11) - a21b11*b12b22 - (ammbmm-a11b11) -
244  * (annbnn-a11b11) + anmbmm*bmnbnn
245  u(3,1) = a(l2,l1)/b(l1,l1)
246  120 CONTINUE
247  130 CONTINUE
248  fail = .true.
249  RETURN
250 C* PUT THE NEGLECTABLE ELEMENTS EQUAL TO ZERO IN
251 C* CASE OF CONVERGENCE
252  140 a(l2,l) = 0.
253  a(l2,l1) = 0.
254  a(l3,l) = 0.
255  a(l3,l1) = 0.
256  b(l1,l) = 0.
257  b(l2,l) = 0.
258  b(l2,l1) = 0.
259  b(l3,l) = 0.
260  b(l3,l1) = 0.
261  b(l3,l2) = 0.
262  RETURN
263  END
exchqz
subroutine exchqz(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL)
Definition: exchqz.f:1
max
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:233
abs
T abs(T x)
Definition: pr-output.cc:3062
Generated on Thu Jun 4 2015 23:30:24 for GNU Octave by   doxygen 1.8.8