C      ALGORITHM 343, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 11, NO. 12, December, 1968, PP.820--826.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Fortran/
#	Fortran/Dp/
#	Fortran/Dp/Drivers/
#	Fortran/Dp/Drivers/Makefile
#	Fortran/Dp/Drivers/driver.f
#	Fortran/Dp/Drivers/res
#	Fortran/Dp/Src/
#	Fortran/Dp/Src/src.f
# This archive created: Fri Mar 10 09:33:00 2006
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Dp'
then
	mkdir 'Dp'
fi
cd 'Dp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver
RESULTS= Res

Objs1= driver.o src.o
driver: $(Objs1)
	$(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS)
Res: driver 
	./driver >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c***********************************************************************
c
cc TOMS343_PRB tests TOMS343.
c
c  Modified:
c
c    20 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none


      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS343_PRB'
      write ( *, '(a)' ) '  Test TOMS algorithm 343, to compute'
      write ( *, '(a)' ) '  the eigenvalues and eigenvectors of'
      write ( *, '(a)' ) '  a real general matrix.'

      call test01

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS343_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'

      write ( *, '(a)' ) ' '

      stop
      end
      subroutine test01

c***********************************************************************
c
cc TEST01 tests EIGENP.
c
c  Modified:
c
c    20 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer nm
      integer n1
      integer n2
      integer n3

      parameter ( nm = 5 )
      parameter ( n1 = 5 )
      parameter ( n2 = 4 )
      parameter ( n3 = 3 )

      real a(nm,nm)
      real a1(5,5)
      real a2(4,4)
      real a3(3,3)
      real evi(nm)
      real evr(nm)
      integer i
      integer indic(nm)
      integer j
      integer k
      integer n
      real t
      real veci(nm,nm)
      real vecr(nm,nm)

      data a1 /
     &  -0.5, 1.0, 0.0, 0.0, 0.0,
     &  -1.0, 0.0, 1.0, 0.0, 0.0,
     &  -1.0, 0.0, 0.0, 1.0, 0.0,
     &  -0.5, 0.0, 0.0, 0.0, 1.0,
     &  -1.0, 0.0, 0.0, 0.0, 0.0 /

      data a2 /
     &  -2.0, -7.0,  0.0, -1.0,
     &   1.0, -5.0, -1.0,  0.0,
     &   1.0, -2.0, -3.0, -1.0,
     &   1.0, -4.0, -2.0,  0.0 /

      data a3 /
     &   1.00, 0.10, 0.00,
     &   0.00, 1.00, 1.00,
     &   0.01, 0.00, 1.00 /

      t = 24.0E+00

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' ) '  Test EIGENP, which computes eigenvalues and'
      write ( *, '(a)' ) '  eigenvectors of a general real matrix.'
      write ( *, '(a)' ) ' '

      do k = 1, 3

        if ( k .eq. 1 ) then

          n = n1

          do i = 1, n
            do j = 1, n
              a(i,j) = a1(i,j)
            end do
          end do
      
        else if ( k .eq. 2 ) then

          n = n2

          do i = 1, n
            do j = 1, n
              a(i,j) = a2(i,j)
            end do
          end do

        else if ( k .eq. 3 ) then

          n = n3

          do i = 1, n
            do j = 1, n
              a(i,j) = a3(i,j)
            end do
          end do
      
        end if

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Matrix A:'
        write ( *, '(a)' ) ' '

        do i = 1, n
          write ( *, '(5f10.4)' ) ( a(i,j), j = 1, n )
        end do

        call eigenp ( n, nm, a, t, evr, evi, vecr, veci, indic )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '       I   INDIC   Real part  Imag part'
        write ( *, '(a)' ) ' '

        do i = 1, n
          write ( *, '(2x,i6,2x,i6,2x,g14.6,2x,g14.6)' ) 
     &      i, indic(i), evr(i), evi(i)
        end do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Real parts of eigenvectors:'
        write ( *, '(a)' ) ' '

        do i = 1, n
          write ( *, '(5f10.4)' ) ( vecr(i,j), j = 1, n )
        end do

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  Imaginary parts of eigenvectors:'
        write ( *, '(a)' ) ' '

        do i = 1, n
          write ( *, '(5f10.4)' ) ( veci(i,j), j = 1, n )
        end do

      end do

      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS343_PRB
  Test TOMS algorithm 343, to compute
  the eigenvalues and eigenvectors of
  a real general matrix.
 
TEST01
  Test EIGENP, which computes eigenvalues and
  eigenvectors of a general real matrix.
 
 
  Matrix A:
 
   -0.5000   -1.0000   -1.0000   -0.5000   -1.0000
    1.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    1.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    1.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    1.0000    0.0000
 
       I   INDIC   Real part  Imag part
 
       1       2    -1.00000         0.00000    
       2       2   -0.250000        0.968246    
       3       2   -0.250000       -0.968246    
       4       2    0.500000        0.866026    
       5       2    0.500000       -0.866026    
 
  Real parts of eigenvectors:
 
   -0.4472    1.0000    1.0000    0.5000    0.5000
    0.4472   -0.2500   -0.2500    1.0000    1.0000
   -0.4472   -0.8750   -0.8750    0.5000    0.5000
    0.4472    0.6875    0.6875   -0.5000   -0.5000
   -0.4472    0.5312    0.5312   -1.0000   -1.0000
 
  Imaginary parts of eigenvectors:
 
    0.0000   -0.0000    0.0000    0.8660   -0.8660
    0.0000   -0.9682    0.9682   -0.0000    0.0000
    0.0000    0.4841   -0.4841   -0.8660    0.8660
    0.0000    0.7262   -0.7262   -0.8660    0.8660
    0.0000   -0.8472    0.8472    0.0000   -0.0000
 
  Matrix A:
 
   -2.0000    1.0000    1.0000    1.0000
   -7.0000   -5.0000   -2.0000   -4.0000
    0.0000   -1.0000   -3.0000   -2.0000
   -1.0000    0.0000   -1.0000    0.0000
 
       I   INDIC   Real part  Imag part
 
       1       2    -4.00000         2.00000    
       2       2    -4.00000        -2.00000    
       3       2    -2.41421         0.00000    
       4       2    0.414214         0.00000    
 
  Real parts of eigenvectors:
 
   -0.2000   -0.2000    0.0000   -0.0000
    1.0000    1.0000   -0.7941    0.4760
    0.2000    0.2000    0.5615    0.3366
   -0.0000   -0.0000    0.2326   -0.8125
 
  Imaginary parts of eigenvectors:
 
   -0.4000    0.4000    0.0000    0.0000
   -0.0000    0.0000    0.0000    0.0000
    0.4000   -0.4000    0.0000    0.0000
   -0.0000    0.0000    0.0000    0.0000
 
  Matrix A:
 
    1.0000    0.0000    0.0100
    0.1000    1.0000    0.0000
    0.0000    1.0000    1.0000
 
       I   INDIC   Real part  Imag part
 
       1       2    0.950000        0.866025E-01
       2       2    0.950000       -0.866025E-01
       3       2     1.10000         0.00000    
 
  Real parts of eigenvectors:
 
   -0.0500   -0.0500   -0.0990
   -0.0500   -0.0500   -0.0990
    1.0000    1.0000   -0.9901
 
  Imaginary parts of eigenvectors:
 
   -0.0866    0.0866    0.0000
    0.0866   -0.0866    0.0000
   -0.0000    0.0000    0.0000
 
TOMS343_PRB
  Normal end of execution.
 
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      SUBROUTINE EIGENP ( N, NM, A, T, EVR, EVI, VECR, VECI, INDIC )

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,*)
      DOUBLE PRECISION D1
      DOUBLE PRECISION D2
      DOUBLE PRECISION D3
      REAL ENORM
      REAL EPS
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER IVEC
      INTEGER IWORK(100)
      INTEGER J
      INTEGER K
      INTEGER K1
      INTEGER KON
      INTEGER L
      INTEGER L1
      INTEGER LOCAL(100)
      INTEGER M
      INTEGER N
      DOUBLE PRECISION PRFACT(100)
      REAL R
      REAL R1
      REAL SUBDIA(100)
      REAL T
      REAL VECI(NM,*)
      REAL VECR(NM,*)
      REAL WORK(100)
      REAL WORK1(100)
      REAL WORK2(100)
C
C  THIS SUBROUTINE FINDS ALL THE EIGENVALULES AND THE
C  EIGENVECTORS OF A REAL GENERAL MATRIX OF ORDER N.
C
C  FIRST IN THE SUBROUTINE SCALE THE MATRIX IS SCALED SO THAT
C  THE CORRESPONDING ROWS AND COLUMNS ARE APPROXIMATELY 
C  BALANCED, AND THEN THE MATRIX IS NORMALIZED SO THAT THE
C  VALUE OF THE EUCLIDEAN NORM OF THE MATRIX IS EQUAL TO ONE.
C
C  THE EIGENVALUES ARE COMPUTED BY THE QR DOUBLE-STEP METHOD
C  IN THE SUBROUTINE HESQR.
C
C  THE EIGENVECTORS ARE COMPUTED BY INVERSE ITERATION IN
C  THE SUBROUTINE REALVE, FOR THE REAL EIGENVALUES, OR IN THE
C  SUBROUTINE COMPVE, FOR THE COMPLEX EIGENVALUES.
C
C  THE ELEMENTS OF THE MATRIX ARE TO BE STORED IN THE FIRST N
C  ROWS AND COLUMNS OF THE TWO DIMENSIONAL ARRAY A.  THE
C  ORIGINAL MATRIX IS DESTROYED BY THE SUBROUTINE.
C
C  N IS THE ORDER OF THE MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A, VECR, VECI, AND THE DIMENSION OF THE ONE
C  DIMENSIONAL ARRAYS EVR, EVI AND INDIC.  THEREFORE, THE
C  CALLING PROGRAM SHOULD CONTAIN THE FOLLOWING DECLARATIONS:
C    REAL A(NM,NN)
C    REAL EVI(NM)
C    REAL EVR(NM)
C    INTEGER INDIC(NM)
C    REAL VECI(NM,NN)
C    REAL VECR(NM,NN)
C  WHERE NM AND NN ARE ANY NUMBERS EQUAL TO OR GREATER THAN N.
C  THE UPPER LIMIT FOR NM IS EQUAL TO 100, BUT MAY BE 
C  INCREASED TO THE VALUE MAX BY REPLACING THE DIMENSION
C  STATEMENTS:
C    INTEGER IWORK(100)
C    INTEGER LOCAL(100)
C    DOUBLE PRECISION PRFACT(100)
C    REAL SUBDIA(100)
C    REAL WORK(100)
C    REAL WORK1(100)
C    REAL WORK2(100)
C  IN THE SUBROUTINE EIGENP BY
C    INTEGER IWORK(MAX)
C    INTEGER LOCAL(MAX)
C    DOUBLE PRECISION PRFACT(MAX)
C    REAL SUBDIA(MAX)
C    REAL WORK(MAX)
C    REAL WORK1(MAX)
C    REAL WORK2(MAX)
C  NM AND NN ARE OF COURSE BOUNDED BY THE SIZE OF THE STORE.
C
C  THE REAL PARAMETER T MUST BE SET EQUAL TO THE NUMBER OF
C  BINARY DIGITS IN THE MANTISSA OF A SINGLE PRECISION
C  FLOATING-POINT NUMBER.
C
C  THE REAL PARTS OF THE N COMPUTED EIGENVALUES WILL BE FOUND
C  IN THE FIRST N PLACES OF THE ARRAY EVI.
C  THE REAL COMPONENTS OF THE NORMALIZED EIGENVECTOR I
C  (I=1,2,...,N) CORRESPONDING TO THE EIGENVALUE STORED IN
C  EVR(I) AND EVI(I) WILL BE FOUND IN THE FIRST N PLACES OF
C  THE COLUMN I OF THE TWO DIMENSIONAL ARRAY VECR AND THE
C  IMAGINARY COMPONENTS IN THE FIRST N PLACES OF THE COLUMN I
C  OF THE TWO DIMENSIONAL ARRAY VECI.
C
C  THE REAL EIGENVECTOR IS NORMALIZED SO THAT THE SUM OF THE
C  SQUARES OF THE COMPONENTS IS EQUAL TO ONE.
C
C  THE COMPLEX EIGENVECTOR IS NORMALIZED SO THAT THE
C  COMPONENT WITH THE LARGEST VALUE IN MODULUS HAS ITS REAL
C  PART EQUAL TO ONE AND THE IMAGINARY PART EQUAL TO ZERO.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE SUBROUTINE
C  EIGENP AS FOLLOWS:
C
C    VALUE OF INDIC(I)    EIGENVALUE I    EIGENVECTOR I
C
C           0             NOT FOUND       NOT FOUND
C           1                 FOUND       NOT FOUND
C           2                 FOUND           FOUND
C
      IF ( N .NE. 1 ) GO TO 1
      EVR(1) = A(1,1)
      EVI(1) = 0.0E+00
      VECR(1,1) = 1.0E+00
      VECI(1,1) = 0.0E+00
      INDIC(1) = 2
      GO TO 25

1     CALL SCALE ( N, NM, A, VECI, PRFACT, ENORM )
C
C  THE COMPUTATION OF THE EIGENVALUES OF THE NORMALIZED
C  MATRIX.
C
      EX = EXP ( - T * ALOG ( 2.0E+00 ) )
      CALL HESQR ( N, NM, A, VECI, EVR, EVI, SUBDIA, INDIC, EPS, EX )
C
C  THE POSSIBLE DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX
C  INTO THE SUBMATRICESS OF LOWER ORDER IS INDICATED IN THE
C  ARRAY LOCAL.  THE DECOMPOSITION OCCURS WHEN SOME
C  SUBDIAGONAL ELEMENTS ARE IN MODULUS LESS THAN A SMALL
C  POSITIVE NUMBER EPS DEFINED IN THE SUBROUTINE HESQR.  THE
C  AMOUNT OF WORK IN THE EIGENVECTOR PROBLEM MAY BE
C  DIMINISHED IN THIS WAY.
C  
      J = N
      I = 1
      LOCAL(1) = 1
      IF ( J .EQ. 1 ) GO TO 4
2     IF ( ABS ( SUBDIA(J-1) ) .GT. EPS ) GO TO 3
      I = I + 1
      LOCAL(I) = 0
3     J = J - 1
      LOCAL(I) = LOCAL(I) + 1
      IF ( J .NE. 1 ) GO TO 2
C
C  THE EIGENVECTOR PROBLEM.
C
4     K = 1
      KON = 0
      L = LOCAL(1)
      M = N
      DO 10 I = 1, N
        IVEC = N - I + 1
        IF ( I .LE. L ) GO TO 5
        K = K + 1
        M = N - L
        L = L + LOCAL(K)
5       IF ( INDIC(IVEC) .EQ. 0 ) GO TO 10
        IF ( EVI(IVEC) .NE. 0.0E+00 ) GO TO 8
C
C  TRANSFER OF AN UPPER-HESSENBERG MATRIX OF ORDER M FROM
C  THE ARRAYS VECI AND SUBDIA INTO THE ARRAY A.
C
        DO 7 K1 = 1, M
          DO 6 L1 = K1, M
            A(K1,L1) = VECI(K1,L1)
6     CONTINUE
          IF ( K1 .EQ. 1 ) GO TO 7
          A(K1,K1-1) = SUBDIA(K1-1)
7       CONTINUE
C
C  THE COMPUTATION OF THE REAL EIGENVECTOR IVEC OF THE UPPER-
C  HESSENBERG MATRIX CORRESPONDING TO THE REAL EIGENVALUE
C  EVR(IVEC).
C
        CALL REALVE ( N, NM, M, IVEC, A, VECR, EVR, EVI, IWORK,
     &    WORK, INDIC, EPS, EX )
        GO TO 10
C
C  THE COMPUTATION OF THE COMPLEX EIGENVECTOR IVEC OF THE
C  UPPER-HESSENBERG MATRIX CORRESPONDING TO THE COMPLEX
C  EIGENVALUE EVR(IVEC)+I*EVI(IVEC).  IF THE VALUE OF KON IS
C  NOT EQUAL TO ZERO, THEN THIS COMPLEX EIGENVECTOR HAS
C  ALREADY BEEN FOUND FROM ITS CONJUGATE.
C
8       IF ( KON .NE. 0 ) GO TO 9
        KON = 1
        CALL COMPVE ( N, NM, M, IVEC, A, VECR, VECI, EVR, EVI, INDIC,
     &    IWORK, SUBDIA, WORK1, WORK2, WORK, EPS, EX )
        GO TO 10
9       KON = 0
10    CONTINUE
C
C  THE RECONSTRUCTION OF THE MATRIX USED IN THE REDUCTION OF
C  MATRIX A TO AN UPPER-HESSENBERG FORM BY HOUSEHOLDER METHOD.
C
      DO 12 I = 1, N
        DO 11 J = I, N
          A(I,J) = 0.0E+00
          A(J,I) = 0.0E+00
11    CONTINUE
        A(I,I) = 1.0E+00
12    CONTINUE

      IF ( N .LE. 2 ) GO TO 15
      M = N - 2
      DO 27 K  =  1, M
        L = K + 1
      DO 26 J  =  2, N
          D1 = 0.0E+00
          DO 13 I = L, N
            D2 = VECI(I,K)
            D1 = D1 + D2 * A(J,I)
13    CONTINUE
          DO 14 I = L, N
            A(J,I) = A(J,I) - VECI(I,K) * D1
14    CONTINUE
26    CONTINUE
27    CONTINUE
C
C  THE COMPUTATION OF THE EIGENVECTORS OF THE ORIGINAL NON-
C  SCALED MATRIX.
C
15    KON = 1
      DO 24 I = 1, N
        L = 0
        IF ( EVI(I) .EQ. 0.0E+00 ) GO TO 16
        L = 1
        IF ( KON .EQ. 0 ) GO TO 16
        KON = 0
        GO TO 24
16      DO 18 J = 1, N
          D1 = 0.0E+00
          D2 = 0.0E+00
          DO 17 K = 1, N
            D3 = A(J,K)
            D1 = D1 + D3 * VECR(K,I)
            IF ( L .EQ. 0 ) GO TO 17
            D2 = D2 + D3 * VECR(K,I-1)
17        CONTINUE
          WORK(J) = D1 / PRFACT(J)
          IF ( L .EQ. 0 ) GO TO 18
          SUBDIA(J) = D2 / PRFACT(J)
18      CONTINUE
C
C  THE NORMALIZATION OF THE EIGENVECTORS AND THE COMPUTATION
C  OF THE EIGENVALUES OF THE ORIGINAL NON-NORMALIZED MATRIX.
C
        IF ( L .EQ. 1 ) GO TO 21
        D1 = 0.0E+00
        DO 19 M = 1, N
          D1 = D1 + WORK(M)**2
19    CONTINUE
        D1 = SQRT ( D1 )

        DO 20 M = 1, N
          VECI(M,I) = 0.0E+00
          VECR(M,I) = WORK(M) / D1
20    CONTINUE

        EVR(I) = EVR(I) * ENORM
        GO TO 24

21      KON = 1
        EVR(I) = EVR(I) * ENORM
        EVR(I-1) = EVR(I)
        EVI(I) = EVI(I) * ENORM
        EVI(I-1) = -EVI(I)

        R = 0.0E+00
        DO 22 J = 1, N
          R1 = WORK(J)**2 + SUBDIA(J)**2
          IF ( R .GE. R1 ) GO TO 22
          R = R1
          L = J
22      CONTINUE

        D3 = WORK(L)
        R1 = SUBDIA(L)
        DO 23 J = 1, N
          D1 = WORK(J)
          D2 = SUBDIA(J)
          VECR(J,I) = ( D1 * D3 + D2 * R1 ) / R
          VECI(J,I) = ( D2 * D3 - D1 * R1 ) / R
          VECR(J,I-1) = VECR(J,I)
          VECI(J,I-1) = -VECI(J,I)
23    CONTINUE

24    CONTINUE

25    RETURN
      END
      SUBROUTINE SCALE ( N, NM, A, H, PRFACT, ENORM )

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,*)
      REAL BOUND1
      REAL BOUND2
      DOUBLE PRECISION COLUMN
      REAL ENORM
      DOUBLE PRECISION FACTOR
      DOUBLE PRECISION FNORM
      REAL H(NM,*)
      INTEGER I
      INTEGER ITER
      INTEGER J
      INTEGER N
      INTEGER NCOUNT
      DOUBLE PRECISION PRFACT(NM)
      DOUBLE PRECISION Q
      DOUBLE PRECISION ROW
C
C  THIS SUBROUTINE STORES THE MATRIX OF THE ORDER N FROM THE
C  ARRAY A INTO THE ARRAY H.  AFTERWARD, THE MATRIX IN THE
C  ARRAY A IS SCALED SO THAT THE QUOTIENT OF THE ABSOLUTE SUM
C  OF THE OFF-DIAGONAL ELEMENTS OF COLUMN I AND THE ABSOLUTE
C  SUM OF THE OFF-DIAGONAL ELEMENTS OF ROW I LIES WITHIN THE
C  VALUES OF BOUND1 AND BOUND2.
C
C  THE COMPONENT I OF THE EIGENVECTOR OBTAINED BY USING THE
C  SCALED MATRIX MUST BE DIVIDED BY THE VALUE FOUND IN THE
C  PRFACT(I) OF THE ARRAY PRFACT.  IN THIS WAY, THE EIGENVECTOR
C  OF THE NON-SCALED MATRIX IS OBTAINED.
C
C  AFTER THE MATRIX IS SCALED, IT IS NORMALIZED SO THAT THE
C  VALUE OF THE EUCLIDEAN NORM IS EQUAL TO ONE.
C  IF THE PROCESS OF SCALING WAS NOT SUCCESSFUL, THE ORIGINAL
C  MATRIX FROM THE ARRAY H WOULD BE STORED BAK INTO A AND
C  THE EIGENPROBLEM WOULD BE SOLVED BY USING THIS MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE ARRAYS A AND H.  NM
C  MUST BE GREATER OR EQUAL TO N.
C
C  THE EIGENVALUES OF THE NORMALIZED MATRIX MUST BE 
C  MULTIPLIED BY THE SCALAR ENORM IN ORDER THAT THEY BECOME
C  THE EIGENVALUES OF THE NON-NORMALIZED MATRIX.
C
      DO 2 I = 1, N
        DO 1 J = 1, N
          H(I,J) = A(I,J)
1     CONTINUE
        PRFACT(I) = 1.0E+00
2     CONTINUE

      BOUND1 = 0.75E+00
      BOUND2 = 1.33E+00
      ITER = 0
3     NCOUNT = 0

      DO 8 I = 1, N

        COLUMN = 0.0E+00
        ROW = 0.0E+00
        DO 4 J = 1, N
          IF ( I .EQ. J ) GO TO 4
          COLUMN = COLUMN + ABS ( A(J,I) )
          ROW = ROW + ABS ( A(I,J) )
4       CONTINUE

        IF ( COLUMN .EQ. 0.0E+00 ) GO TO 5
        IF ( ROW .EQ. 0.0E+00 ) GO TO 5
        Q = COLUMN / ROW
        IF ( Q .LT. BOUND1 ) GO TO 6
        IF ( Q .GT. BOUND2 ) GO TO 6
5       NCOUNT = NCOUNT + 1
        GO TO 8
6       FACTOR = DSQRT ( Q )

        DO 7 J = 1, N
          IF ( I .EQ. J ) GO TO 7
          A(I,J) = A(I,J) * FACTOR
          A(J,I) = A(J,I) / FACTOR
7       CONTINUE

        PRFACT(I) = PRFACT(I) * FACTOR

8     CONTINUE

      ITER = ITER + 1
      IF ( ITER .GT. 30 ) GO TO 11
      IF ( NCOUNT .LT. N ) GO TO 3

      FNORM = 0.0E+00
      DO 14 I  =  1, N
        DO 9 J = 1, N
          Q = A(I,J)
          FNORM = FNORM + Q * Q
9     CONTINUE
14    CONTINUE
      FNORM = DSQRT ( FNORM )

      DO 16 I  =  1, N
        DO 10 J = 1, N
          A(I,J) = A(I,J) / FNORM
10    CONTINUE
16    CONTINUE

      ENORM = FNORM
      GO TO 13

11    DO 15 I  =  1, N
        DO 12 J = 1, N
          A(I,J) = H(I,J)
12    CONTINUE
15    CONTINUE

      ENORM = 1.0E+00

13    RETURN
      END
      SUBROUTINE HESQR ( N, NM, A, H, EVR, EVI, SUBDIA, INDIC, EPS, EX )

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,*)
      REAL EPS
      REAL EX
      REAL H(NM,*)
      REAL EVI(NM)
      REAL EVR(NM)
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER MAXST
      INTEGER M1
      INTEGER N
      INTEGER NS
      REAL R
      DOUBLE PRECISION S
      REAL SHIFT
      DOUBLE PRECISION SR
      DOUBLE PRECISION SR2
      REAL SUBDIA(NM)
      REAL T
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      DOUBLE PRECISION Z
C
C  THIS SUBROUTINE FINDS ALL THE EIGENVALUES OF A REAL
C  GENERAL MATRIX.  THE ORIGINAL MATRIX A OF ORDER N IS
C  REDUCED TO THE UPPER-HESSENBERG FORM H BY MEANS OF
C  SIMILARITY TRANSFORMATIONS (HOUSEHOLDER METHOD).  THE MATRIX
C  H IS PRESERVED IN THE UPPER HALF OF THE ARRAY H AND IN THE
C  ARRAY SUBDIA.  THE SPECIAL VECTORS USED IN THE DEFINITION
C  OF THE HOUSEHOLDER TRANSFORMATION MATRICES ARE STORED IN
C  THE LOWER PART OF THE ARRAY H.
C
C  NM IS THE FIRST DIMENSION OF THE ARRAYS A AND H.  NM MUST
C  BE EQUAL TO OR GREATER THAN N.
C
C  THE REAL PARTS OF THE N EIGENVALUES WILL BE FOUND IN THE
C  FIRST N PLACES OF THE ARRAY EVR, AND
C  THE IMAGINARY PARTS IN THE FIRST N PLACES OF THE ARRAY EVI.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C
C    VALUE OF INDIC(I)  EIGENVALUE I
C           0             NOT FOUND
C           1                 FOUND
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF H)*EX, WHERE
C  EX = 2**(-T).  T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
C  REDUCTION OF THE MATRIX A TO AN UPPER-HESSENBERG FORM H.
C  THERE ARE N-2 STEPS.
C
      IF ( N - 2 ) 14, 1, 2

1     SUBDIA(1) = A(2,1)
      GO TO 14

2     M = N - 2

      DO 12 K = 1, M

        L = K + 1

        S = 0.0E+00
        DO 3 I = L, N
          H(I,K) = A(I,K)
          S = S + ABS ( A(I,K) )
3     CONTINUE

        IF ( S .NE. ABS ( A(K+1,K) ) ) GO TO 4
        SUBDIA(K) = A(K+1,K)
        H(K+1,K) = 0.0E+00
        GO TO 12
4       SR2 = 0.0E+00

        DO 5 I = L, N
          SR = A(I,K)
          SR = SR / S
          A(I,K) = SR
          SR2 = SR2 + SR * SR
5     CONTINUE

        SR = DSQRT ( SR2 )
        IF ( A(L,K) .LT. 0.0 ) GO TO 6
        SR = -SR
6       SR2 = SR2 - SR * A(L,K)
        A(L,K) = A(L,K) - SR
        H(L,K) = H(L,K) - SR * S
        SUBDIA(K) = SR * S
        X = S * DSQRT ( SR2 )

        DO 7 I = L, N
          H(I,K) = H(I,K) / X
          SUBDIA(I) = A(I,K) / SR2
7     CONTINUE
C
C  PREMULTIPLICATION BY THE MATRIX PR.
C
      DO 39 J  =  L, N
          SR = 0.0E+00
          DO 8 I = L, N
            SR = SR + A(I,K) * A(I,J)
8     CONTINUE
          DO 9 I = L, N
            A(I,J) = A(I,J) - SUBDIA(I) * SR
9     CONTINUE
39    CONTINUE
C
C  POSTMULTIPLICATION BY THE MATRIX PR.
C
      DO 38 J  =  1, N
          SR = 0.0E+00
          DO 10 I = L, N
            SR = SR + A(J,I) * A(I,K)
10    CONTINUE
          DO 11 I = L, N
            A(J,I) = A(J,I) - SUBDIA(I) * SR
11    CONTINUE
38    CONTINUE

12    CONTINUE

      DO 13 K = 1, M
        A(K+1,K) = SUBDIA(K)
13    CONTINUE
C
C  TRANSFER OF THE UPPER HALF OF THE MATRIX A INTO THE
C  ARRAY H, AND THE CALCULATION OF THE SMALL POSITIVE NUMBER EPS.
C
      SUBDIA(N-1) = A(N,N-1)
14    EPS = 0.0E+00
      DO 40 K  =  1, N
        INDIC(K) = 0
        IF ( K .NE. N ) EPS = EPS + SUBDIA(K)**2
        DO 15 I = K, N
          H(K,I) = A(K,I)
          EPS = EPS + A(K,I)**2
15    CONTINUE
40    CONTINUE

      EPS = EX * SQRT ( EPS )
C
C  THE QR ITERATIVE PROCESS.  THE UPPER-HESSENBERG MATRIX H IS
C  REDUCED TO THE UPPER-MODIFIED TRIANGULAR FORM.
C
C  DETERMINATION OF THE SHIFT OF ORIGIN FOR THE FIRST STEP OF
C  THE QR ITERATIVE PROCESS.
C
      SHIFT = A(N,N-1)
      IF ( N .LE. 2 ) SHIFT = 0.0E+00
      IF ( A(N,N) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      IF ( A(N-1,N) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      IF ( A(N-1,N-1) .NE. 0.0E+00 ) SHIFT = 0.0E+00
      M = N
      NS = 0
      MAXST = N * 10
C
C  TESTING IF THE UPPER HALF OF THE MATRIX IS EQUAL TO ZERO.
C  IF IT IS EQUAL TO ZERO, THE QR PROCESS IS NOT NECESSARY.
C
      DO 41 I  =  2, N
        DO 16 K = I, N
          IF ( A(I-1,K) .NE. 0.0E+00 ) GO TO 18
16        CONTINUE
41    CONTINUE

      DO 17 I = 1, N
        INDIC(I) = 1
        EVR(I) = A(I,I)
        EVI(I) = 0.0E+00
17    CONTINUE

      GO TO 37
C
C  START THE MAIN LOOP OF THE QR PROCESS.
C
18    K = M - 1
      M1 = K
      I = K
C
C  FIND ANY DECOMPOSITIONS OF THE MATRIX.
C  JUMP TO 34 IF THE LAST SUBMATRIX OF THE DECOMPOSITION IS
C  OF THE ORDER ONE.
C  JUMP TO 35 IF THE LAST SUBMATRIX OF THE DECOMPOSITION IS
C  OF THE ORDER TWO.
C
      IF ( K ) 37, 34, 19
19    IF ( ABS ( A(M,K) ) .LE. EPS ) GO TO 34
      IF ( M - 2 .EQ. 0 ) GO TO 35
20    I = I - 1
      IF ( ABS ( A(K,I) ) .LE. EPS ) GO TO 21
      K = I
      IF ( K .GT. 1 ) GO TO 20
21    IF ( K .EQ. M1 ) GO TO 35
C
C  TRANSFORMATION OF THE MATRIX OF THE ORDER GREATER THAN TWO.
C
      S = A(M,M) + A(M1,M1) + SHIFT
      SR = A(M,M) * A(M1,M1) - A(M,M1)*A(M1,M) + 0.25E+00 * SHIFT**2
      A(K+2,K) = 0.0E+00
C
C  CALCULATE X1, Y1, Z1 FOR THE SUBMATRIX OBTAINED BY THE
C  DECOMPOSITION.
C
      X = A(K,K) * ( A(K,K) - S ) + A(K,K+1) * A(K+1,K) + SR
      Y = A(K+1,K) * ( A(K,K) + A(K+1,K+1) - S )
      R = DABS ( X ) + DABS ( Y )
      IF ( R .EQ. 0.0E+00 ) SHIFT = A(M,M-1)
      IF ( R .EQ. 0.0E+00 ) GO TO 21
      Z = A(K+2,K+1) * A(K+1,K)
      SHIFT = 0.0E+00
      NS = NS + 1
C
C  THE LOOP FOR ONE STEP OF THE QR PROCESS.
C
      DO 33 I = K, M1
        IF ( I .EQ. K ) GO TO 22
C
C  CALCULATE XR, YR, ZR.
C
        X = A(I,I-1)
        Y = A(I+1,I-1)
        Z = 0.0E+00
        IF ( I + 2 .GT. M ) GO TO 22
        Z = A(I+2,I-1)
22      SR2 = DABS ( X ) + DABS ( Y ) + DABS ( Z )
        IF ( SR2 .EQ. 0.0E+00 ) GO TO 23
        X = X / SR2
        Y = Y / SR2
        Z = Z / SR2
23      S = DSQRT ( X * X + Y * Y + Z * Z )
        IF ( X .LT. 0.0E+0 ) GO TO 24
        S = -S
24      IF ( I .EQ. K ) GO TO 25
        A(I,I-1) = S * SR2
25      IF ( SR2 .NE. 0.0E+00 ) GO TO 26
        IF ( I + 3 .GT. M ) GO TO 33
        GO TO 32
26      SR = 1.0E+00 - X / S
        S = X - S
        X = Y / S
        Y = Z / S
C
C  PREMULTIPLICATION BY THE MATRIX PR.
C
        DO 28 J = I, M
          S = A(I,J) + A(I+1,J) * X
          IF ( I + 2 .GT. M ) GO TO 27
          S = S + A(I+2,J) * Y
27        S = S * SR
          A(I,J) = A(I,J) - S
          A(I+1,J) = A(I+1,J) - S * X
          IF ( I + 2 .GT. M ) GO TO 28
          A(I+2,J) = A(I+2,J) - S * Y
28      CONTINUE
C
C  POSTMULTIPLICATION BY THE MATRIX PR.
C
        L = I + 2
        IF ( I .LT. M1 ) GO TO 29
        L = M

29      DO 31 J = K, L
          S = A(J,I) + A(J,I+1) * X
          IF ( I + 2 .GT. M ) GO TO 30
          S = S + A(J,I+2) * Y
30        S = S * SR
          A(J,I) = A(J,I) - S
          A(J,I+1) = A(J,I+1) - S * X
          IF ( I + 2 .GT. M ) GO TO 31
          A(J,I+2) = A(J,I+2) - S * Y
31      CONTINUE

        IF ( I + 3 .GT. M ) GO TO 33
        S = -A(I+3,I+2) * Y * SR
32      A(I+3,I) = S
        A(I+3,I+1) = S * X
        A(I+3,I+2) = S * Y + A(I+3,I+2)
33    CONTINUE

      IF ( NS .GT. MAXST ) GO TO 37
      GO TO 18
C
C  COMPUTE THE LAST EIGENVALUE.
C
34    EVR(M) = A(M,M)
      EVI(M) = 0.0E+00
      INDIC(M) = 1
      M = K
      GO TO 18
C
C  COMPUTE THE EIGENVALUES OF THE LAST 2X2 MATRIX OBTAINED BY
C  THE DECOMPOSITION.
C
35    R = 0.5E+00 * ( A(K,K) + A(M,M) )
      S = 0.5E+00 * ( A(M,M) - A(K,K) )
      S = S * S + A(K,M) * A(M,K)
      INDIC(K) = 1
      INDIC(M) = 1
      IF ( S .LT. 0.0E+00 ) GO TO 36
      T = DSQRT ( S )
      EVR(K) = R - T   
      EVR(M) = R + T
      EVI(K) = 0.0E+00
      EVI(M) = 0.0E+00
      M = M - 2
      GO TO 18

36    T = DSQRT ( -S )
      EVR(K) = R
      EVI(K) = T
      EVR(M) = R
      EVI(M) = -T
      M = M - 2
      GO TO 18

37    RETURN
      END
      SUBROUTINE REALVE ( N, NM, M, IVEC, A, VECR, EVR, EVI,
     &  IWORK, WORK, INDIC, EPS, EX )

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,*)
      REAL BOUND
      REAL EPS
      REAL EVALUE
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      INTEGER I
      INTEGER INDIC(NM)
      INTEGER ITER
      INTEGER IVEC
      INTEGER IWORK(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER N
      INTEGER NS
      REAL PREVIS
      REAL R
      REAL R1
      DOUBLE PRECISION S
      DOUBLE PRECISION SR
      REAL T
      REAL VECR(NM,*)
      REAL WORK(NM)
C
C  THIS SUBROUTINE FINDS THE REAL EIGENVECTOR OF THE REAL
C  UPPER-HESSENBERG MATRIX IN THE ARRAY A, CORRESPONDING TO
C  THE REAL EIGENVALUE STORED IN EVR(IVEC).  THE INVERSE
C  ITERATION METHOD IS USED.
C
C  NOTE THE MATRIX IN A IS DESTROYED BY THE SUBROUTINE.
C
C  N IS THE ORDER OF THE UPPER HESSENBERG MATRIX.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A AND VECR.  NM MUST BE EQUAL TO OR GREATER THAN N.
C
C  M IS THE ORDER OF THE SUBMATRIX OBTAINED BY A SUITABLE
C  DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX IF SOME
C  SUBDIAGONAL ELEMENTS ARE EQUAL TO ZERO.  THE VALUE OF M IS
C  CHOSEN SO THAT THE LAST N-M COMPONENTS OF THE EIGENVECTOR
C  ARE ZERO.
C
C  IVEC GIVES THE POSITION OF THE EIGENVALUE IN THE ARRAY EVR
C  FOR WHICH THE CORRESPONDING EIGENVECTOR IS COMPUTED.
C  
C  THE ARRAY EVI WOULD CONTAIN THE IMAGINARY PARTS OF THE N
C  EIGENVALUES IF THEY EXISTED.
C
C  THE M COMPONENTS OF THE COMPUTED REAL EIGENVECTOR WILL BE
C  FOUND IN THE FIRST M PLACES OF THE COLUMN IVEC OF THE TWO
C  DIMENSIONAL ARRAY VECR.
C
C  IWORK AND WORK ARE THE WORKING STORES USED DURING THE
C  GAUSSIAN ELIMINATION AND BACKSUBSTITUTION PROCESS.
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C    VALUE OF INDIC(I)  EIGENVECTOR I
C           1             NOT FOUND
C           2                 FOUND
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF A)*EX, WHERE
C  EX = 2**(-T), T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
      VECR(1,IVEC) = 1.0E+00
      IF ( M .EQ. 1 ) GO TO 24
C
C  SMALL PERTURBATION OF EQUAL EIGENVALUES TO OBTAIN A FULL
C  SET OF EIGENVECTORS.
C
      EVALUE = EVR(IVEC)
      IF ( IVEC .EQ. M ) GO TO 2
      K = IVEC + 1
      R = 0.0E+00
      DO 1 I = K, M
        IF ( EVALUE .NE. EVR(I) ) GO TO 1
        IF ( EVI(I) .NE. 0.0E+00 ) GO TO 1
        R = R + 3.0E+00
1     CONTINUE
      EVALUE = EVALUE + R * EX

2     DO 3 K = 1, M
        A(K,K) = A(K,K) - EVALUE
3     CONTINUE
C
C  GAUSSIAN ELIMINATION OF THE UPPER-HESSENBERG MATRIX A.  ALL
C  ROW INTERCHANGES ARE INDICATED IN THE ARRAY IWORK.  ALL THE
C  MULTIPLIERS ARE STORED AS THE SUBDIAGONAL ELEMENTS OF A.
C
      K = M - 1
      DO 8 I = 1, K
        L = I + 1
        IWORK(I) = 0
        IF ( A(I+1,I) .NE. 0.0E+00 ) GO TO 4
        IF ( A(I,I) .NE. 0.0E+00 ) GO TO 8
        A(I,I) = EPS
        GO TO 8
4       IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+1,I) ) ) GO TO 6
        IWORK(I) = 1
        DO 5 J = I, M
          R = A(I,J)
          A(I,J) = A(I+1,J)
          A(I+1,J) = R
5     CONTINUE
6       R = -A(I+1,I) / A(I,I)
        A(I+1,I) = R
        DO 7 J = L, M
          A(I+1,J) = A(I+1,J) + R * A(I,J)
7     CONTINUE
8     CONTINUE

      IF ( A(M,M) .NE. 0.0E+00 ) GO TO 9
      A(M,M) = EPS
C
C  THE VECTOR (1,1,...,1) IS STORED IN THE PLACE OF THE RIGHT
C  HAND SIDE COLUMN VECTOR.
C
9     DO 11 I = 1, N
        IF ( I .GE. M ) GO TO 10
        WORK(I) = 1.0E+00
        GO TO 11
10      WORK(I) = 0.0E+00
11    CONTINUE
C
C  THE INVERSE ITERATION IS PERFORMED ON THE MATRIX UNTIL THE
C  INFINITY NORM OF THE RIGHT-HAND SIDE VECTOR IS GREATER
C  THAN THE BOUND DEFINED AS 0.01 / ( N * EX ).
C
      BOUND = 0.01E+00 / ( EX * FLOAT ( N ) )
      NS = 0
      ITER = 1
C
C  THE BACKSUBSTITUTION.
C
12    R = 0.0E+00
      DO 15 I = 1, M
        J = M - I + 1
        S = WORK(J)
        IF ( J .EQ. M ) GO TO 14
        L = J + 1
        DO 13 K = L, M
          SR = WORK(K)
          S = S - SR * A(J,K)
13    CONTINUE
14      WORK(J) = S / A(J,J)
        T = ABS ( WORK(J) )
        IF ( R .GE. T ) GO TO 15
        R = T
15    CONTINUE
C
C  THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR FOR THE NEW
C  ITERATION STEP.
C
      DO 16 I = 1, M
        WORK(I) = WORK(I) / R
16    CONTINUE
C
C  THE COMPUTATION OF THE RESIDUALS AND COMPARISON OF THE
C  RESIDUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE
C  ITERATION.  IF THE INFINITY NORM OF THE RESIDUAL VECTOR
C  IS GREATER THAN THE INFINITY NORM OF THE PREVIOUS RESIDUAL
C  VECTOR, THE COMPUTED EIGENVECTOR OF THE PREVIOUS STEP IS
C  TAKEN AS THE FINAL EIGENVECTOR.
C
      R1 = 0.0E+00
      DO 18 I = 1, M
        T = 0.0E+00
        DO 17 J = I, M
          T = T + A(I,J) * WORK(J)
17    CONTINUE
        T = ABS ( T )
        IF ( R1 .GE. T ) GO TO 18
        R1 = T
18    CONTINUE

      IF ( ITER .EQ. 1 ) GO TO 19
      IF ( PREVIS .LE. R1 ) GO TO 24
19    DO 20 I = 1, M
        VECR(I,IVEC) = WORK(I)
20    CONTINUE
      PREVIS = R1
      IF ( NS .EQ. 1 ) GO TO 24
      IF ( ITER .GT. 6 ) GO TO 25
      ITER = ITER + 1
      IF ( R .LT. BOUND ) GO TO 21
      NS = 1
C
C  GAUSSIAN ELIMINATION OF THE RIGHT HAND SIDE VECTOR.
C
21    K = M - 1
      DO 23 I = 1, K
        R = WORK(I+1)
        IF ( IWORK(I) .EQ. 0 ) GO TO 22
        WORK(I+1) = WORK(I) + WORK(I+1) * A(I+1,I)
        WORK(I) = R
        GO TO 23
22      WORK(I+1) = WORK(I+1) + WORK(I) * A(I+1,I)
23    CONTINUE

      GO TO 12

24    INDIC(IVEC) = 2
25    IF ( M .EQ. N ) GO TO 27
      J = M + 1
      DO 26 I = J, N
        VECR(I,IVEC) = 0.0E+00
26    CONTINUE
27    RETURN
      END
      SUBROUTINE COMPVE ( N, NM, M, IVEC, A, VECR, H, EVR, EVI, INDIC,
     &  IWORK, SUBDIA, WORK1, WORK2, WORK, EPS, EX )

      IMPLICIT NONE

      INTEGER NM

      REAL A(NM,*)
      REAL B
      REAL BOUND
      DOUBLE PRECISION D
      DOUBLE PRECISION D1
      REAL EPS
      REAL ETA
      REAL EVI(NM)
      REAL EVR(NM)
      REAL EX
      REAL FKSI
      REAL H(NM,*)
      INTEGER I
      INTEGER I1
      INTEGER I2
      INTEGER INDIC(NM)
      INTEGER ITER
      INTEGER IVEC
      INTEGER IWORK(NM)
      INTEGER J
      INTEGER K
      INTEGER L
      INTEGER M
      INTEGER N
      INTEGER NS
      REAL PREVIS
      REAL R
      REAL S
      REAL SUBDIA(NM)
      REAL U
      REAL V
      REAL VECR(NM,*)
      REAL WORK(NM)
      REAL WORK1(NM)
      REAL WORK2(NM)
C
C  THIS SUBROUTINE FINDS THE COMPLEX EIGENVECTOR OF THE REAL
C  UPPER-HESSENBERG MATRIX OF ORDER N CORRESPONDING TO THE
C  COMPLEX EIGENVALUE WITH THE REAL PART IN EVR(IVEC) AND THE
C  CORRESPONDING IMAGINARY PART IN EVI(IVEC).  THE INVERSE
C  ITERATION METHOD IS USED, MODIFIED TO AVOID THE USE OF
C  COMPLEX ARITHMETIC.
C
C  THE MATRIX ON WHICH THE INVERSE ITERATION IS PERFORMED IS
C  BUILT UP INTO THE ARRAY A BY USING THE UPPER-HESSENBERG
C  MATRIX PRESERVED IN THE UPPER HALF OF THE ARRAY H AND IN
C  THE ARRAY SUBDIA.
C
C  NM DEFINES THE FIRST DIMENSION OF THE TWO DIMENSIONAL
C  ARRAYS A, VECR, AND H.  NM MUST BE EQUAL TO OR GREATER
C  THAN N.
C
C  M IS THE ORDER OF THE SUBMATRIX OBTAINED BY A SUITABLE
C  DECOMPOSITION OF THE UPPER-HESSENBERG MATRIX IF SOME
C  SUBDIAGONAL ELEMENTS ARE EQUAL TO ZERO.  THE VALUE OF M IS
C  CHOSEN SO THAT THE LAST N-M COMPONENTS OF THE COMPLEX
C  EIGENVECTOR ARE ZERO.
C
C  THE REAL PARTS OF THE FIRST M COMPONENTS OF THE COMPUTED
C  COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST M PLACES OF
C  THE COLUMN WHOSE TOP ELEMENT IS VECR(1,IVEC), AND THE
C  CORRESPONDING IMAGINARY PARTS OF THE FIRST M COMPONENTS OF
C  THE COMPLEX EIGENVECTOR WILL BE FOUND IN THE FIRST M
C  PLACES OF THE COLUMN WHOSE TOP ELEMENT IS VECR(1,IVEC-1).
C
C  THE ARRAY INDIC INDICATES THE SUCCESS OF THE ROUTINE AS
C  FOLLOWS:
C    VALUE OF INDIC(I)  EIGENVECTOR I
C           1             NOT FOUND
C           2                 FOUND
C
C  THE ARRAYS IWORK, WORK1, WORK2 AND WORK ARE THE WORKING
C  STORES USED DURING THE INVERSE ITERATION PROCESS.
C
C  EPS IS A SMALL POSITIVE NUMBER THAT NUMERICALLY REPRESENTS
C  ZERO IN THE PROGRAM.  EPS = (EUCLIDEAN NORM OF H)*EX, WHERE
C  EX = 2**(-T).  T IS THE NUMBER OF BINARY DIGITS IN THE
C  MANTISSA OF A FLOATING POINT NUMBER.
C
      FKSI = EVR(IVEC)
      ETA = EVI(IVEC)
C
C  THE MODIFICATION OF THE EIGENVALUE ( FKSI + I * ETA ) IF MORE
C  EIGENVALUES ARE EQUAL.
C
      IF ( IVEC .EQ. M ) GO TO 2
      K = IVEC + 1
      R = 0.0E+00
      DO 1 I = K, M
        IF ( FKSI .NE. EVR(I) ) GO TO 1
        IF ( ABS ( ETA ) .NE. ABS ( EVI(I) ) ) GO TO 1
        R = R + 3.0E+00
1     CONTINUE
      R = R * EX
      FKSI = FKSI + R
      ETA = ETA + R
C
C  THE MATRIX ((H-FKSI*I)*(H-FKSI*I)+(ETA*ETA)*I) IS
C  STORED INTO THE ARRAY A.
C
2     R = FKSI * FKSI + ETA * ETA
      S = 2.0E+00 * FKSI
      L = M - 1

      DO 5 I = 1, M
        DO 4 J = I, M
          D = 0.0E+00
          A(J,I) = 0.0E+00
          DO 3 K = I, J
            D = D + H(I,K) * H(K,J)
3     CONTINUE
          A(I,J) = D - S * H(I,J)
4     CONTINUE
        A(I,I) = A(I,I) + R
5     CONTINUE

      DO 9 I = 1, L
        R = SUBDIA(I)
        A(I+1,I) = -S * R
        I1 = I + 1
        DO 6 J = 1, I1
          A(J,I) = A(J,I) + R * H(J,I+1)
6     CONTINUE
        IF ( I .EQ. 1 ) GO TO 7
        A(I+1,I-1) = R * SUBDIA(I-1)
7       DO 8 J = I, M
          A(I+1,J) = A(I+1,J) + R * H(I,J)
8     CONTINUE
9     CONTINUE

C
C  THE GAUSSIAN ELIMINATION OF THE MATRIX
C  ((H-FKSI*I)*(H-FKSI*I)+(ETA*ETA)*I) IN THE ARRAY A.  THE
C  ROW INTERCHANGES THAT OCCUR ARE INDICATED IN THE ARRAY
C  IWORK.  ALL THE MULTIPLIERS ARE STORED IN THE FIRST AND IN
C  THE SECOND SUBDIAGONAL OF THE ARRAY A.
C
      K = M - 1
      DO 18 I = 1, K
        I1 = I + 1
        I2 = I + 2
        IWORK(I) = 0
        IF ( I .EQ. K ) GO TO 10
        IF ( A(I+2,I) .NE. 0.0E+00 ) GO TO 11
10      IF ( A(I+1,I) .NE. 0.0E+00 ) GO TO 11
        IF ( A(I,I) .NE. 0.0E+00 ) GO TO 18
        A(I,I) = EPS
        GO TO 18

11      IF ( I .EQ. K ) GO TO 12
        IF ( ABS ( A(I+1,I) ) .GE. ABS ( A(I+2,I) ) ) GO TO 12
        IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+2,I) ) ) GO TO 16
        L = I + 2
        IWORK(I) = 2
        GO TO 13

12      IF ( ABS ( A(I,I) ) .GE. ABS ( A(I+1,I) ) ) GO TO 15
        L = I + 1
        IWORK(I) = 1

13      DO 14 J = I, M
          R = A(I,J)
          A(I,J) = A(L,J)
          A(L,J) = R
14    CONTINUE

15      IF ( I .NE. K ) GO TO 16
        I2 = I1
16    DO 48 L  =  I1, I2
          R = -A(L,I) / A(I,I)
          A(L,I) = R
          DO 17 J = I1, M
            A(L,J) = A(L,J) + R * A(I,J)
17    CONTINUE
48    CONTINUE
18    CONTINUE

      IF ( A(M,M) .NE. 0.0E+00 ) GO TO 19
      A(M,M) = EPS
C
C  THE VECTOR (1,1,...,1) IS STORED INTO THE RIGHT-HAND SIDE
C  VECTORS VECR(*,IVEC) AND VECR(*,IVEC-1), REPRESENTING THE
C  COMPLEX RIGHT HAND SIDE VECTOR.
C
19    DO 21 I = 1, N
        IF ( I .GT. M ) GO TO 20
        VECR(I,IVEC) = 1.0E+00
        VECR(I,IVEC-1) = 1.0E+00
        GO TO 21
20      VECR(I,IVEC) = 0.0E+00
        VECR(I,IVEC-1) = 0.0E+00
21    CONTINUE
C
C  THE INVERSE ITERATION IS PERFORMED ON THE MATRIX UNTIL THE
C  INFINITY NORM OF THE RIGHT-HAND SIDE IS GREATER
C  THAN THE BOUND DEFINED AS 0.01/(N*EX).
C
      BOUND = 0.01E+00 / ( EX * FLOAT ( N ) )
      NS = 0
      ITER = 1
      DO 22 I = 1, M
        WORK(I) = H(I,I) - FKSI
22    CONTINUE
C
C  THE SEQUENCE OF THE COMPLEX VECTORS Z(S) = P(S) + I * Q(S) AND
C  W(S+1) = U(S+1)+I*V(S+1) IS GIVEN BY THE RELATIONS
C  (A-FKSI-I*ETA)*I)*W(S+1) = Z(S) AND
C  Z(S+1) = W(S+1)/MAX(W(S+1)).
C  THE FINAL W(S) IS TAKEN AS THE COMPUTED EIGENVECTOR.
C
C  THE COMPUTATION OF THE RIGHT-HAND SIDE VECTOR
C  (A-FKSI*I)*P(S)-ETA*Q(S).  A IS AN UPPER-HESSENBERG MATRIX.
C
23    DO 27 I = 1, M
        D = WORK(I) * VECR(I,IVEC)
        IF ( I .EQ. 1 ) GO TO 24
        D = D + SUBDIA(I-1) * VECR(I-1,IVEC)
24      L = I + 1
        IF ( L .GT. M ) GO TO 26
        DO 25 K = L, M
          D = D + H(I,K) * VECR(K,IVEC)
25    CONTINUE
26      VECR(I,IVEC-1) = D - ETA * VECR(I,IVEC-1)
27    CONTINUE
C
C  GAUSSIAN ELIMINATION OF THE RIGHT-HAND SIDE VECTOR.
C
      K = M - 1
      DO 28 I = 1, K
        L = I + IWORK(I)
        R = VECR(L,IVEC-1)
        VECR(L,IVEC-1) = VECR(I,IVEC-1)
        VECR(I,IVEC-1) = R
        VECR(I+1,IVEC-1) = VECR(I+1,IVEC-1) + A(I+1,I) * R
        IF ( I .EQ. K ) GO TO 28
        VECR(I+2,IVEC-1) = VECR(I+2,IVEC-1) + A(I+2,I) * R
28    CONTINUE
C
C  THE COMPUTATION OF THE REAL PART U(S+1) OF THE COMPLEX
C  VECTOR W(S+1).  THE VECTOR U(S+1) IS OBTAINED AFTER THE
C  BACKSUBSTITUTION.
C
      DO 31 I = 1, M
        J = M - I + 1
        D = VECR(J,IVEC-1)
        IF ( J .EQ. M ) GO TO 30
        L = J + 1
        DO 29 K = L, M
          D1 = A(J,K)
          D = D - D1 * VECR(K,IVEC-1)
29    CONTINUE
30      VECR(J,IVEC-1) = D / A(J,J)
31    CONTINUE
C
C  THE COMPUTATION OF THE IMAGINARY PART V(S+1) OF THE VECTOR
C  W(S+1), WHERE V(S+1) = (P(S)-(A-FKSI*I)*U(S+1))/ETA.
C
      DO 35 I = 1, M
        D = WORK(I) * VECR(I,IVEC-1)
        IF ( I .EQ. 1 ) GO TO 32
        D = D + SUBDIA(I-1) * VECR(I-1,IVEC-1)
32      L = I + 1
        IF ( L .GT. M ) GO TO 34
        DO 33 K = L, M
          D = D + H(I,K) * VECR(K,IVEC-1)
33    CONTINUE
34      VECR(I,IVEC) = ( VECR(I,IVEC) - D ) / ETA
35    CONTINUE
C
C  THE COMPUTATION OF (INFINITY NORM OF W(S+1))**2.
C
      L = 1
      S = 0.0E+00
      DO 36 I = 1, M
        R = VECR(I,IVEC)**2 + VECR(I,IVEC-1)**2
        IF ( R .LE. S ) GO TO 36
        S = R
        L = I
36    CONTINUE
C
C  THE COMPUTATION OF THE VECTOR Z(S+1), WHERE Z(S+1) = W(S+1)/
C  (COMPONENT OF W(S+1) WITH THE LARGEST ABSOLUTE VALUE).
C
      U = VECR(L,IVEC-1)
      V = VECR(L,IVEC)
      DO 37 I = 1, M
        B = VECR(I,IVEC)
        R = VECR(I,IVEC-1)
        VECR(I,IVEC) = ( R * U + B * V ) / S
        VECR(I,IVEC-1) = ( B * U - R * V ) / S
37    CONTINUE
C
C  THE COMPUTATION OF THE RESIDUALS AND COMPARISON OF THE
C  RESIDUALS OF THE TWO SUCCESSIVE STEPS OF THE INVERSE
C  ITERATION.  IF THE INFINITY NORM OF THE RESIDUAL VECTOR IS
C  GREATER THAN THE INFINITY NORM OF THE PREVIOUS RESIDUAL
C  VECTOR, THE COMPUTED VECTOR OF THE PREVIOUS STEP IS TAKEN
C  AS THE COMPUTED APPROXIMATION TO THE EIGENVECTOR.
C
      B = 0.0E+00
      DO 41 I = 1, M
        R = WORK(I) * VECR(I,IVEC-1) - ETA * VECR(I,IVEC)
        U = WORK(I) * VECR(I,IVEC) + ETA * VECR(I,IVEC-1)
        IF ( I .EQ. 1 ) GO TO 38
        R = R + SUBDIA(I-1) * VECR(I-1,IVEC-1)
        U = U + SUBDIA(I-1) * VECR(I-1,IVEC)
38      L = I + 1
        IF ( L .GT. M ) GO TO 40
        DO 39 J = L, M
          R = R + H(I,J) * VECR(J,IVEC-1)
          U = U + H(I,J) * VECR(J,IVEC)
39    CONTINUE
40      U = R * R + U * U
        IF ( B .GE. U ) GO TO 41
        B = U
41    CONTINUE

      IF ( ITER .EQ. 1 ) GO TO 42
      IF ( PREVIS .LE. B ) GO TO 44
42    DO 43 I = 1, N
        WORK1(I) = VECR(I,IVEC)
        WORK2(I) = VECR(I,IVEC-1)
43    CONTINUE

      PREVIS = B
      IF ( NS .EQ. 1 ) GO TO 46
      IF ( ITER .GT. 6 ) GO TO 47
      ITER = ITER + 1
      IF ( BOUND .GT. SQRT ( S ) ) GO TO 23
      NS = 1
      GO TO 23

44    DO 45 I = 1, N
        VECR(I,IVEC) = WORK1(I)
        VECR(I,IVEC-1) = WORK2(I)
45    CONTINUE
46    INDIC(IVEC-1) = 2
      INDIC(IVEC) = 2

47    RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0