GNU Octave  3.8.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
ssubsp.f
Go to the documentation of this file.
1  SUBROUTINE ssubsp(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
2  INTEGER nmax, n, ftest, ndim, ind(n)
3  LOGICAL fail
4  REAL a(nmax,n), b(nmax,n), z(nmax,n), eps
5 C*
6 C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A
7 C* WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS ROUTINE REORDERS THE DIAGONAL
8 C* BLOCKS ALONG WITH THEIR GENERALIZED EIGENVALUES BY CONSTRUCTING EQUI-
9 C* VALENCE TRANSFORMATIONS QT AND ZT. THE ROW TRANSFORMATION ZT IS ALSO
10 C* PERFORMED ON THE GIVEN (INITIAL) TRANSFORMATION Z (RESULTING FROM A
11 C* POSSIBLE PREVIOUS STEP OR INITIALIZED WITH THE IDENTITY MATRIX).
12 C* AFTER REORDERING, THE EIGENVALUES INSIDE THE REGION SPECIFIED BY THE
13 C* FUNCTION FTEST APPEAR AT THE TOP. IF NDIM IS THEIR NUMBER THEN THE
14 C* NDIM FIRST COLUMNS OF Z SPAN THE REQUESTED SUBSPACE. DSUBSP REQUIRES
15 C* THE SUBROUTINE EXCHQZ AND THE INTEGER FUNCTION FTEST WHICH HAS TO BE
16 C* PROVIDED BY THE USER. THE PARAMETERS IN THE CALLING SEQUENCE ARE :
17 C* (STARRED PARAMETERS ARE 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 REORDERED.
22 C* *Z UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN
23 C* TRANSFORMATION ZT.
24 C* FTEST(LS,ALPHA,BETA,S,P) AN INTEGER FUNCTION DESCRIBING THE
25 C* SPECTRUM OF THE DEFLATING SUBSPACE TO BE COMPUTED:
26 C* WHEN LS=1 FTEST CHECKS IF ALPHA/BETA IS IN THAT SPECTRUM
27 C* WHEN LS=2 FTEST CHECKS IF THE TWO COMPLEX CONJUGATE
28 C* ROOTS WITH SUM S AND PRODUCT P ARE IN THAT SPECTRUM
29 C* IF THE ANSWER IS POSITIVE, FTEST=1, OTHERWISE FTEST=-1
30 C* EPS THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT
31 C* *NDIM AN INTEGER GIVING THE DIMENSION OF THE COMPUTED
32 C* DEFLATING SUBSPACE
33 C* *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN,
34 C* TRUE OTHERWISE (WHEN SEXCHQZ FAILS)
35 C* *IND AN INTEGER WORKING ARRAY OF DIMENSION AT LEAST N
36 C*
37  INTEGER l, ls, ls1, ls2, l1, ll, num, is, l2i, l2k, i, k, ii,
38  * istep, ifirst
39  REAL s, p, d, alpha, beta
40  fail = .true.
41  ndim = 0
42  num = 0
43  l = 0
44  ls = 1
45 C*** CONSTRUCT ARRAY IND(I) WHERE :
46 C*** IABS(IND(I)) IS THE SIZE OF THE BLOCK I
47 C*** SIGN(IND(I)) INDICATES THE LOCATION OF ITS EIGENVALUES
48 C*** (AS DETERMINED BY FTEST).
49 C*** NUM IS THE NUMBER OF ELEMENTS IN THIS ARRAY
50  DO 30 ll=1,n
51  l = l + ls
52  IF (l.GT.n) go to 40
53  l1 = l + 1
54  IF (l1.GT.n) go to 10
55  IF (a(l1,l).EQ.0.) go to 10
56 C* HERE A 2X2 BLOCK IS CHECKED *
57  ls = 2
58  d = b(l,l)*b(l1,l1)
59  s = (a(l,l)*b(l1,l1)+a(l1,l1)*b(l,l)-a(l1,l)*b(l,l1))/d
60  p = (a(l,l)*a(l1,l1)-a(l,l1)*a(l1,l))/d
61  is = ftest(ls,alpha,beta,s,p)
62  go to 20
63 C* HERE A 1X1 BLOCK IS CHECKED *
64  10 ls = 1
65  is = ftest(ls,a(l,l),b(l,l),s,p)
66  20 num = num + 1
67  IF (is.EQ.1) ndim = ndim + ls
68  ind(num) = ls*is
69  30 CONTINUE
70 C*** REORDER BLOCKS SUCH THAT THOSE WITH POSITIVE VALUE
71 C*** OF IND(.) APPEAR FIRST.
72  40 l2i = 1
73  DO 100 i=1,num
74  IF (ind(i).GT.0) go to 90
75 C* IF A NEGATIVE IND(I) IS ENCOUNTERED, THEN SEARCH FOR THE FIRST
76 C* POSITIVE IND(K) FOLLOWING ON IT
77  l2k = l2i
78  DO 60 k=i,num
79  IF (ind(k).LT.0) go to 50
80  go to 70
81  50 l2k = l2k - ind(k)
82  60 CONTINUE
83 C* IF THERE ARE NO POSITIVE INDICES FOLLOWING ON A NEGATIVE ONE
84 C* THEN STOP
85  go to 110
86 C* IF A POSITIVE IND(K) FOLLOWS ON A NEGATIVE IND(I) THEN
87 C* INTERCHANGE BLOCK K BEFORE BLOCK I BY PERFORMING K-I SWAPS
88  70 istep = k - i
89  ls2 = ind(k)
90  l = l2k
91  DO 80 ii=1,istep
92  ifirst = k - ii
93  ls1 = -ind(ifirst)
94  l = l - ls1
95  CALL sexchqz(nmax, n, a, b, z, l, ls1, ls2, eps, fail)
96  IF (fail) RETURN
97  ind(ifirst+1) = ind(ifirst)
98  80 CONTINUE
99  ind(i) = ls2
100  90 l2i = l2i + ind(i)
101  100 CONTINUE
102  110 fail = .false.
103  RETURN
104  END