Mercurial > hg > octave-nkf

annotate libcruft/ordered-qz/ssubsp.f @ 10828:322f43e0e170

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Grammarcheck .txi documentation files.
author Rik <octave@nomad.inbox5.com>
date Wed, 28 Jul 2010 12:45:04 -0700
parents 96ba591be50f
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7793
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 SUBROUTINE SSUBSP(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
2 INTEGER NMAX, N, FTEST, NDIM, IND(N)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 LOGICAL FAIL
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 REAL A(NMAX,N), B(NMAX,N), Z(NMAX,N), EPS
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
5 C*
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
6 C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 C* WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS ROUTINE REORDERS THE DIAGONAL
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
8 C* BLOCKS ALONG WITH THEIR GENERALIZED EIGENVALUES BY CONSTRUCTING EQUI-
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 C* VALENCE TRANSFORMATIONS QT AND ZT. THE ROW TRANSFORMATION ZT IS ALSO
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 C* PERFORMED ON THE GIVEN (INITIAL) TRANSFORMATION Z (RESULTING FROM A
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 C* POSSIBLE PREVIOUS STEP OR INITIALIZED WITH THE IDENTITY MATRIX).
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 C* AFTER REORDERING, THE EIGENVALUES INSIDE THE REGION SPECIFIED BY THE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
13 C* FUNCTION FTEST APPEAR AT THE TOP. IF NDIM IS THEIR NUMBER THEN THE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 C* NDIM FIRST COLUMNS OF Z SPAN THE REQUESTED SUBSPACE. DSUBSP REQUIRES
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 C* THE SUBROUTINE EXCHQZ AND THE INTEGER FUNCTION FTEST WHICH HAS TO BE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 C* PROVIDED BY THE USER. THE PARAMETERS IN THE CALLING SEQUENCE ARE :
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 C* (STARRED PARAMETERS ARE ALTERED BY THE SUBROUTINE)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
18 C*
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 C* NMAX THE FIRST DIMENSION OF A, B AND Z
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 C* N THE ORDER OF A, B AND Z
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 C* *A,*B THE MATRIX PAIR WHOSE BLOCKS ARE TO BE REORDERED.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
22 C* *Z UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 C* TRANSFORMATION ZT.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
24 C* FTEST(LS,ALPHA,BETA,S,P) AN INTEGER FUNCTION DESCRIBING THE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 C* SPECTRUM OF THE DEFLATING SUBSPACE TO BE COMPUTED:
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 C* WHEN LS=1 FTEST CHECKS IF ALPHA/BETA IS IN THAT SPECTRUM
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 C* WHEN LS=2 FTEST CHECKS IF THE TWO COMPLEX CONJUGATE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
28 C* ROOTS WITH SUM S AND PRODUCT P ARE IN THAT SPECTRUM
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 C* IF THE ANSWER IS POSITIVE, FTEST=1, OTHERWISE FTEST=-1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
30 C* EPS THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 C* *NDIM AN INTEGER GIVING THE DIMENSION OF THE COMPUTED
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 C* DEFLATING SUBSPACE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 C* *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 C* TRUE OTHERWISE (WHEN SEXCHQZ FAILS)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 C* *IND AN INTEGER WORKING ARRAY OF DIMENSION AT LEAST N
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
36 C*
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
37 INTEGER L, LS, LS1, LS2, L1, LL, NUM, IS, L2I, L2K, I, K, II,
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 * ISTEP, IFIRST
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 REAL S, P, D, ALPHA, BETA
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 FAIL = .TRUE.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 NDIM = 0
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 NUM = 0
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 L = 0
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 LS = 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 C*** CONSTRUCT ARRAY IND(I) WHERE :
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 C*** IABS(IND(I)) IS THE SIZE OF THE BLOCK I
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 C*** SIGN(IND(I)) INDICATES THE LOCATION OF ITS EIGENVALUES
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 C*** (AS DETERMINED BY FTEST).
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 C*** NUM IS THE NUMBER OF ELEMENTS IN THIS ARRAY
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 DO 30 LL=1,N
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 L = L + LS
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 IF (L.GT.N) GO TO 40
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 L1 = L + 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 IF (L1.GT.N) GO TO 10
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
55 IF (A(L1,L).EQ.0.) GO TO 10
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 C* HERE A 2X2 BLOCK IS CHECKED *
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
57 LS = 2
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
58 D = B(L,L)*B(L1,L1)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
59 S = (A(L,L)*B(L1,L1)+A(L1,L1)*B(L,L)-A(L1,L)*B(L,L1))/D
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 P = (A(L,L)*A(L1,L1)-A(L,L1)*A(L1,L))/D
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
61 IS = FTEST(LS,ALPHA,BETA,S,P)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
62 GO TO 20
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
63 C* HERE A 1X1 BLOCK IS CHECKED *
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 10 LS = 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
65 IS = FTEST(LS,A(L,L),B(L,L),S,P)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
66 20 NUM = NUM + 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 IF (IS.EQ.1) NDIM = NDIM + LS
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
68 IND(NUM) = LS*IS
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 30 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
70 C*** REORDER BLOCKS SUCH THAT THOSE WITH POSITIVE VALUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 C*** OF IND(.) APPEAR FIRST.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
72 40 L2I = 1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 DO 100 I=1,NUM
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
74 IF (IND(I).GT.0) GO TO 90
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
75 C* IF A NEGATIVE IND(I) IS ENCOUNTERED, THEN SEARCH FOR THE FIRST
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
76 C* POSITIVE IND(K) FOLLOWING ON IT
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 L2K = L2I
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 DO 60 K=I,NUM
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
79 IF (IND(K).LT.0) GO TO 50
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 GO TO 70
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 50 L2K = L2K - IND(K)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 60 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
83 C* IF THERE ARE NO POSITIVE INDICES FOLLOWING ON A NEGATIVE ONE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 C* THEN STOP
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 GO TO 110
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 C* IF A POSITIVE IND(K) FOLLOWS ON A NEGATIVE IND(I) THEN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 C* INTERCHANGE BLOCK K BEFORE BLOCK I BY PERFORMING K-I SWAPS
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 70 ISTEP = K - I
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 LS2 = IND(K)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 L = L2K
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 DO 80 II=1,ISTEP
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 IFIRST = K - II
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 LS1 = -IND(IFIRST)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 L = L - LS1
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 CALL SEXCHQZ(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
96 IF (FAIL) RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 IND(IFIRST+1) = IND(IFIRST)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 80 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 IND(I) = LS2
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 90 L2I = L2I + IND(I)
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 100 CONTINUE
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 110 FAIL = .FALSE.
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
103 RETURN
96ba591be50f Add some more support for single precision to libcruft functions
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 END