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