ssubsp.f

Go to the documentation of this file.
00001       SUBROUTINE SSUBSP(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
00002       INTEGER NMAX, N, FTEST, NDIM, IND(N)
00003       LOGICAL FAIL
00004       REAL A(NMAX,N), B(NMAX,N), Z(NMAX,N), EPS
00005 C*
00006 C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A
00007 C* WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS ROUTINE REORDERS THE DIAGONAL
00008 C* BLOCKS ALONG WITH THEIR GENERALIZED EIGENVALUES BY CONSTRUCTING EQUI-
00009 C* VALENCE TRANSFORMATIONS QT AND ZT. THE ROW TRANSFORMATION ZT IS ALSO
00010 C* PERFORMED ON THE GIVEN (INITIAL) TRANSFORMATION Z (RESULTING FROM A
00011 C* POSSIBLE PREVIOUS STEP OR INITIALIZED WITH THE IDENTITY MATRIX).
00012 C* AFTER REORDERING, THE EIGENVALUES INSIDE THE REGION SPECIFIED BY THE
00013 C* FUNCTION FTEST APPEAR AT THE TOP. IF NDIM IS THEIR NUMBER THEN THE
00014 C* NDIM FIRST COLUMNS OF Z SPAN THE REQUESTED SUBSPACE. DSUBSP REQUIRES
00015 C* THE SUBROUTINE EXCHQZ AND THE INTEGER FUNCTION FTEST WHICH HAS TO BE
00016 C* PROVIDED BY THE USER. THE PARAMETERS IN THE CALLING SEQUENCE ARE :
00017 C* (STARRED PARAMETERS ARE ALTERED BY THE SUBROUTINE)
00018 C*
00019 C*    NMAX     THE FIRST DIMENSION OF A, B AND Z
00020 C*    N        THE ORDER OF A, B AND Z
00021 C*   *A,*B     THE MATRIX PAIR WHOSE BLOCKS ARE TO BE REORDERED.
00022 C*   *Z        UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN
00023 C*             TRANSFORMATION ZT.
00024 C*    FTEST(LS,ALPHA,BETA,S,P) AN INTEGER FUNCTION DESCRIBING THE
00025 C*             SPECTRUM OF THE DEFLATING SUBSPACE TO BE COMPUTED:
00026 C*             WHEN LS=1 FTEST CHECKS IF ALPHA/BETA IS IN THAT SPECTRUM
00027 C*             WHEN LS=2 FTEST CHECKS IF THE TWO COMPLEX CONJUGATE
00028 C*             ROOTS WITH SUM S AND PRODUCT P ARE IN THAT SPECTRUM
00029 C*             IF THE ANSWER IS POSITIVE, FTEST=1, OTHERWISE FTEST=-1
00030 C*    EPS      THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT
00031 C*   *NDIM     AN INTEGER GIVING THE DIMENSION OF THE COMPUTED
00032 C*             DEFLATING SUBSPACE
00033 C*   *FAIL     A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN,
00034 C*             TRUE OTHERWISE (WHEN SEXCHQZ FAILS)
00035 C*   *IND      AN INTEGER WORKING ARRAY OF DIMENSION AT LEAST N
00036 C*
00037       INTEGER L, LS, LS1, LS2, L1, LL, NUM, IS, L2I, L2K, I, K, II,
00038      * ISTEP, IFIRST
00039       REAL S, P, D, ALPHA, BETA
00040       FAIL = .TRUE.
00041       NDIM = 0
00042       NUM = 0
00043       L = 0
00044       LS = 1
00045 C*** CONSTRUCT ARRAY IND(I) WHERE :
00046 C***     IABS(IND(I)) IS THE SIZE OF THE BLOCK I
00047 C***     SIGN(IND(I)) INDICATES THE LOCATION OF ITS EIGENVALUES
00048 C***                  (AS DETERMINED BY FTEST).
00049 C*** NUM IS THE NUMBER OF ELEMENTS IN THIS ARRAY
00050       DO 30 LL=1,N
00051         L = L + LS
00052         IF (L.GT.N) GO TO 40
00053         L1 = L + 1
00054         IF (L1.GT.N) GO TO 10
00055         IF (A(L1,L).EQ.0.) GO TO 10
00056 C* HERE A 2X2  BLOCK IS CHECKED *
00057         LS = 2
00058         D = B(L,L)*B(L1,L1)
00059         S = (A(L,L)*B(L1,L1)+A(L1,L1)*B(L,L)-A(L1,L)*B(L,L1))/D
00060         P = (A(L,L)*A(L1,L1)-A(L,L1)*A(L1,L))/D
00061         IS = FTEST(LS,ALPHA,BETA,S,P)
00062         GO TO 20
00063 C* HERE A 1X1  BLOCK IS CHECKED *
00064    10   LS = 1
00065         IS = FTEST(LS,A(L,L),B(L,L),S,P)
00066    20   NUM = NUM + 1
00067         IF (IS.EQ.1) NDIM = NDIM + LS
00068         IND(NUM) = LS*IS
00069    30 CONTINUE
00070 C***  REORDER BLOCKS SUCH THAT THOSE WITH POSITIVE VALUE
00071 C***    OF IND(.) APPEAR FIRST.
00072    40 L2I = 1
00073       DO 100 I=1,NUM
00074         IF (IND(I).GT.0) GO TO 90
00075 C* IF A NEGATIVE IND(I) IS ENCOUNTERED, THEN SEARCH FOR THE FIRST
00076 C* POSITIVE IND(K) FOLLOWING ON IT
00077         L2K = L2I
00078         DO 60 K=I,NUM
00079           IF (IND(K).LT.0) GO TO 50
00080           GO TO 70
00081    50     L2K = L2K - IND(K)
00082    60   CONTINUE
00083 C* IF THERE ARE NO POSITIVE INDICES FOLLOWING ON A NEGATIVE ONE
00084 C* THEN STOP
00085         GO TO 110
00086 C* IF A POSITIVE IND(K) FOLLOWS ON A NEGATIVE IND(I) THEN
00087 C* INTERCHANGE BLOCK K BEFORE BLOCK I BY PERFORMING K-I SWAPS
00088    70   ISTEP = K - I
00089         LS2 = IND(K)
00090         L = L2K
00091         DO 80 II=1,ISTEP
00092           IFIRST = K - II
00093           LS1 = -IND(IFIRST)
00094           L = L - LS1
00095           CALL SEXCHQZ(NMAX, N, A, B, Z, L, LS1, LS2, EPS, FAIL)
00096           IF (FAIL) RETURN
00097           IND(IFIRST+1) = IND(IFIRST)
00098    80   CONTINUE
00099         IND(I) = LS2
00100    90   L2I = L2I + IND(I)
00101   100 CONTINUE
00102   110 FAIL = .FALSE.
00103       RETURN
00104       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines