dsubsp.f

Go to the documentation of this file.

      SUBROUTINE dsubsp(NMAX, N, A, B, Z, FTEST, EPS, NDIM, FAIL, IND)
      INTEGER nmax, n, ftest, ndim, ind(n)
      LOGICAL fail
      DOUBLE PRECISION a(nmax,n), b(nmax,n), z(nmax,n), eps
C*
C* GIVEN THE UPPER TRIANGULAR MATRIX B AND UPPER HESSENBERG MATRIX A
C* WITH 1X1 OR 2X2 DIAGONAL BLOCKS, THIS ROUTINE REORDERS THE DIAGONAL
C* BLOCKS ALONG WITH THEIR GENERALIZED EIGENVALUES BY CONSTRUCTING EQUI-
C* VALENCE TRANSFORMATIONS QT AND ZT. THE ROW TRANSFORMATION ZT IS ALSO
C* PERFORMED ON THE GIVEN (INITIAL) TRANSFORMATION Z (RESULTING FROM A
C* POSSIBLE PREVIOUS STEP OR INITIALIZED WITH THE IDENTITY MATRIX).
C* AFTER REORDERING, THE EIGENVALUES INSIDE THE REGION SPECIFIED BY THE
C* FUNCTION FTEST APPEAR AT THE TOP. IF NDIM IS THEIR NUMBER THEN THE
C* NDIM FIRST COLUMNS OF Z SPAN THE REQUESTED SUBSPACE. DSUBSP REQUIRES
C* THE SUBROUTINE EXCHQZ AND THE INTEGER FUNCTION FTEST WHICH HAS TO BE
C* PROVIDED BY THE USER. THE PARAMETERS IN THE CALLING SEQUENCE ARE :
C* (STARRED PARAMETERS ARE ALTERED BY THE SUBROUTINE)
C*
C*    NMAX     THE FIRST DIMENSION OF A, B AND Z
C*    N        THE ORDER OF A, B AND Z
C*   *A,*B     THE MATRIX PAIR WHOSE BLOCKS ARE TO BE REORDERED.
C*   *Z        UPON RETURN THIS ARRAY IS MULTIPLIED BY THE COLUMN
C*             TRANSFORMATION ZT.
C*    FTEST(LS,ALPHA,BETA,S,P) AN INTEGER FUNCTION DESCRIBING THE
C*             SPECTRUM OF THE DEFLATING SUBSPACE TO BE COMPUTED:
C*             WHEN LS=1 FTEST CHECKS IF ALPHA/BETA IS IN THAT SPECTRUM
C*             WHEN LS=2 FTEST CHECKS IF THE TWO COMPLEX CONJUGATE
C*             ROOTS WITH SUM S AND PRODUCT P ARE IN THAT SPECTRUM
C*             IF THE ANSWER IS POSITIVE, FTEST=1, OTHERWISE FTEST=-1
C*    EPS      THE REQUIRED ABSOLUTE ACCURACY OF THE RESULT
C*   *NDIM     AN INTEGER GIVING THE DIMENSION OF THE COMPUTED
C*             DEFLATING SUBSPACE
C*   *FAIL     A LOGICAL VARIABLE WHICH IS FALSE ON A NORMAL RETURN,
C*             TRUE OTHERWISE (WHEN EXCHQZ FAILS)
C*   *IND      AN INTEGER WORKING ARRAY OF DIMENSION AT LEAST N
C*
      INTEGER l, ls, ls1, ls2, l1, ll, num, is, l2i, l2k, i, k, ii,
     * istep, ifirst
      DOUBLE PRECISION s, p, d, alpha, beta
      fail = .true.
      ndim = 0
      num = 0
      l = 0
      ls = 1
C*** CONSTRUCT ARRAY IND(I) WHERE :
C***     IABS(IND(I)) IS THE SIZE OF THE BLOCK I
C***     SIGN(IND(I)) INDICATES THE LOCATION OF ITS EIGENVALUES
C***                  (AS DETERMINED BY FTEST).
C*** NUM IS THE NUMBER OF ELEMENTS IN THIS ARRAY
      DO 30 ll=1,n
        l = l + ls
        IF (l.GT.n) go to 40
        l1 = l + 1
        IF (l1.GT.n) go to 10
        IF (a(l1,l).EQ.0.) go to 10
C* HERE A 2X2  BLOCK IS CHECKED *
        ls = 2
        d = b(l,l)*b(l1,l1)
        s = (a(l,l)*b(l1,l1)+a(l1,l1)*b(l,l)-a(l1,l)*b(l,l1))/d
        p = (a(l,l)*a(l1,l1)-a(l,l1)*a(l1,l))/d
        is = ftest(ls,alpha,beta,s,p)
        go to 20
C* HERE A 1X1  BLOCK IS CHECKED *
   10   ls = 1
        is = ftest(ls,a(l,l),b(l,l),s,p)
   20   num = num + 1
        IF (is.EQ.1) ndim = ndim + ls
        ind(num) = ls*is
   30 CONTINUE
C***  REORDER BLOCKS SUCH THAT THOSE WITH POSITIVE VALUE
C***    OF IND(.) APPEAR FIRST.
   40 l2i = 1
      DO 100 i=1,num
        IF (ind(i).GT.0) go to 90
C* IF A NEGATIVE IND(I) IS ENCOUNTERED, THEN SEARCH FOR THE FIRST
C* POSITIVE IND(K) FOLLOWING ON IT
        l2k = l2i
        DO 60 k=i,num
          IF (ind(k).LT.0) go to 50
          go to 70
   50     l2k = l2k - ind(k)
   60   CONTINUE
C* IF THERE ARE NO POSITIVE INDICES FOLLOWING ON A NEGATIVE ONE
C* THEN STOP
        go to 110
C* IF A POSITIVE IND(K) FOLLOWS ON A NEGATIVE IND(I) THEN
C* INTERCHANGE BLOCK K BEFORE BLOCK I BY PERFORMING K-I SWAPS
   70   istep = k - i
        ls2 = ind(k)
        l = l2k
        DO 80 ii=1,istep
          ifirst = k - ii
          ls1 = -ind(ifirst)
          l = l - ls1
          CALL exchqz(nmax, n, a, b, z, l, ls1, ls2, eps, fail)
          IF (fail) RETURN
          ind(ifirst+1) = ind(ifirst)
   80   CONTINUE
        ind(i) = ls2
   90   l2i = l2i + ind(i)
  100 CONTINUE
  110 fail = .false.
      RETURN
      END