C      ALGORITHM 384, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 13, NO. 6, June, 1970, PP.369--371.
#! /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/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Wed Jan 18 20:30:12 2006
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
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 TOMS384_PRB tests the TOMS384 routine SYMQR.
c
c  Modified:
c
c    09 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS384_PRB'
      write ( *, '(a)' ) '  Test SYMQR, which computes eigenvalues'
      write ( *, '(a)' ) '  and eigenvectors of a symmetric matrix.'

      call test01

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

      stop
      end
      subroutine test01

c*******************************************************************************
c
cc TEST01 tests SYMQR on a full matrix.
c
c  Modified:
c
c    09 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 5 )

      real a(n,n)
      logical abscnv
      real angle
      real d(n)
      real e(n)
      real eps
      real exact
      integer fail
      integer i
      integer j
      real k0
      integer na
      real pi
      logical trd
      logical vec

      pi = 3.14159265E+00

      k0 = 0.0E+00
      na = n
      eps = 0.00001E+00
      abscnv = .false.
      vec = .true.
      trd = .false.

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
c
c  Set the values of the matrix A.
c
      do i = 1, n
        do j = 1, n
          a(i,j) = min ( i, j )
        end do
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  The matrix A:'
      write ( *, '(a)' ) ' '

      do i = 1, n
        write ( *, '(2x,5f10.4)' ) ( a(i,j), j = 1, n )
      end do
c
c  Set the values of the right hand side vector B.
c
      call symqr ( a, d, e, k0, n, na, eps, abscnv, vec, trd,
     &  fail )

      if ( fail .ne. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  FAIL was returned nonzero.'
        write ( *, '(a,i6)' ) '  FAIL = ', fail
        write ( *, '(a)' ) '  Only eigendata FAIL+1, ..., N can be'
        write ( *, '(a)' ) '  relied on.'
      end if
c
c  Print the results.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Eigenvalues:'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '    Computed          Exact'
      write ( *, '(a)' ) ' '

      do i = 1, n
        angle = real ( 2 * i - 1 ) * pi / real ( 2 * n + 1 )
        exact = 0.5E+00 / ( 1.0E+00 - cos ( angle ) )
        write ( *, '(2x,g14.6,2x,g14.6)' ) d(i), exact
      end do

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Eigenvector matrix:'
      write ( *, '(a)' ) ' '

      do i = 1, n
        write ( *, '(2x,5f10.4)' ) ( a(i,j), j = 1, n )
      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'
 
TOMS384_PRB
  Test SYMQR, which computes eigenvalues
  and eigenvectors of a symmetric matrix.
 
TEST01
 
  The matrix A:
 
      1.0000    1.0000    1.0000    1.0000    1.0000
      1.0000    2.0000    2.0000    2.0000    2.0000
      1.0000    2.0000    3.0000    3.0000    3.0000
      1.0000    2.0000    3.0000    4.0000    4.0000
      1.0000    2.0000    3.0000    4.0000    5.0000
 
  Eigenvalues:
 
    Computed          Exact
 
     12.3435         12.3435    
     1.44869         1.44869    
    0.582964        0.582964    
    0.353254        0.353253    
    0.271554        0.271554    
 
  Eigenvector matrix:
 
      0.1699   -0.4557    0.5969    0.5485    0.3260
      0.3260   -0.5984    0.1653   -0.4538   -0.5455
      0.4559   -0.3183   -0.5255   -0.1797    0.5819
      0.5483    0.1560   -0.3675    0.6146   -0.4288
      0.5970    0.5563    0.4788   -0.3358    0.1549
 
TOMS384_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 SYMQR ( A, D, E, K0, N, NA, EPS, ABSCNV, VEC, TRD, 
     &  FAIL )
C
C  EXPLANATION OF THE PARAMETERS IN THE CALLING SEQUENCE:
C
C    A      A DOUBLE DIMENSIONED ARRAY.  IF THE MATRIX IS NOT
C           INITIALLY TRIDIAGONAL, IT IS CONTAINED IN THE LOWER
C           TRIANGLE OF A.  IF EIGENVECTORS ARE NOT REQUESTED,
C           THE LOWER TRIANGLE OF A IS DESTROYED WHILE THE
C           ELEMENTS ABOVE THE DIAGONAL ARE LEFT UNDISTURBED.
C           IF EIGENVECTORS ARE REQUESTED, THEY ARE RETURNED IN THE
C           COLUMNS OF A.
C
C    D      A SINGLY SUBSCRIPTED ARRAY.  IF THE MATRIX IS
C           INITIALLY TRIDIAGONAL, D CONTAINS ITS DIAGONAL 
C           ELEMENTS.  ON RETURN, D CONTAINS THE EIGENVALUES OF
C           THE MATRIX.
C
C    E      A SINGLY SUBSCRIPTED ARRAY.  IF THE MATRIX IS
C           INITIALLY TRIDIAGONAL, E CONTAINS ITS OFF-DIAGONAL
C           ELEMENTS.  UPON RETURN, E(I) CONTAINS THE NUMBER OF
C           ITERATIONS REQUIRED TO COMPUTE THE APPROXIMATE
C           EIGENVALUE D(I).
C
C    K0     A REAL VARIABLE CONTAINING AN INITIAL ORIGIN SHIFT TO
C           BE USED UNTIL THE COMPUTED SHIFTS SETTLE DOWN.
C
C    N      AN INTEGER VARIABLE CONTAINING THE FIRST DIMENSION
C           OF THE ARRAY A.
C
C    EPS    A REAL VARIABLE CONTAINING A CONVERGENCE TOLERANCE.
C
C    ABSCNV A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE. IF
C           THE ABSOLUTE CONVERGENCE CRITERION IS TO BE USED
C           OR THE VALUE .FALSE. IF THE RELATIVE CRITERION
C           IS TO BE USED.
C
C    VEC    A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE. IF
C           EIGENVECTORS ARE TO BE COMPUTED AND RETURNED IN
C           THE ARRAY A AND OTHERWISE CONTAINING THE VALUE
C           .FALSE..
C
C    TRD    A LOGICAL VARIABLE CONTAINING THE VALUE .TRUE.
C           IF THE MATRIX IS TRIDIAGONAL AND LOCATED IN THE ARRAYS
C           D AND E AND OTHERWISE CONTAINING THE VALUE .FALSE..
C
C    FAIL   AN INTEGER VARIABLE CONTAINING AN ERROR SIGNAL.
C           ON RETURN, THE EIGENVALUES IN D(FAIL+1), ..., D(N)
C           AND THEIR CORRESPONDING EIGENVECTORS MAY BE PRESUMED
C           ACCURATE.
C
      IMPLICIT NONE

      INTEGER N
      INTEGER NA

      REAL A(NA,N)
      LOGICAL ABSCNV
      REAL C
      REAL CB
      REAL CC
      REAL CD
      REAL CON
      REAL D(N)
      REAL E(N)
      REAL EPS
      INTEGER FAIL
      INTEGER I
      INTEGER I1
      INTEGER J
      REAL K
      REAL K0
      REAL K1
      REAL K2
      INTEGER L
      INTEGER L1
      INTEGER LL
      INTEGER LL1
      REAL RMAX
      REAL NINF
      INTEGER NL
      INTEGER NM1
      INTEGER NM2
      INTEGER NNL
      REAL NORM
      INTEGER NU
      INTEGER NUM1
      REAL P
      REAL Q
      REAL R
      REAL S
      REAL S2
      LOGICAL SHFT
      REAL RSUM
      REAL SUM1
      REAL TEMP
      REAL TEST
      REAL TITTER
      LOGICAL TRD
      LOGICAL VEC

      TITTER = 50.0E+00
      NM1 = N - 1
      NM2 = N - 2
      NINF = 0.0E+00
      FAIL = 0
C
C  SIGNAL ERROR IF N IS NOT POSITIVE.
C
      IF ( N .GT. 0 ) GO TO 1
      FAIL = -1
      RETURN
C
C  SPECIAL TREATMENT FOR A MATRIX OF ORDER ONE.
C
1     IF ( N .GT. 1 ) GO TO 5
      IF ( .NOT. TRD ) D(1) = A(1,1)
      IF ( VEC ) A(1,1) = 1.0E+00
      FAIL = 0
      RETURN
C
C  IF THE MATRIX IS TRIDIAGONAL, SKIP THE REDUCTION.C
C
5     IF ( TRD ) GO TO 100
      IF ( N .EQ. 2 ) GO TO 80
C
C  REDUCE THE MATRIX TO TRIDIAGONAL FORM BY HOUSEHOLDER'S METHOD.
C
      DO 70 L = 1, NM2

        L1 = L + 1
        D(L) = A(L,L)
        RMAX = 0.0E+00

        DO 10 I = L1, N
          RMAX = MAX ( RMAX, ABS ( A(I,L) ) )
10    CONTINUE

        IF ( RMAX .NE. 0.0 ) GO TO 13
        E(L) = 0.0E+00
        A(L,L) = 1.0E+00
        GO TO 70
13      RSUM = 0.0E+00

        DO 17 I = L1, N
          A(I,L) = A(I,L) / RMAX
          RSUM = RSUM + A(I,L)**2
17    CONTINUE

        S2 = RSUM
        S2 = SQRT ( S2 )
        IF ( A(L1,L) .LT. 0.0 ) S2 = -S2
        E(L) = -S2 * RMAX
        A(L1,L) = A(L1,L) + S2
        A(L,L) = S2 * A(L1,L)

        SUM1 = 0.0E+00
        DO 50 I = L1, N
          RSUM = 0.0E+00
          DO 20 J = L1, I
            RSUM = RSUM + A(I,J) * A(J,L)
20    CONTINUE
          IF ( I .EQ. N ) GO TO 40
          I1 = I + 1
          DO 30 J = I1, N
            RSUM = RSUM + A(J,L) * A(J,I)
30    CONTINUE
40        E(I) = RSUM / A(L,L)
          SUM1 = SUM1 + A(I,L) * E(I)
50    CONTINUE

        CON = 0.5E+00 * SUM1 / A(L,L)
      DO 59 I  =  L1, N
          E(I) = E(I) - CON * A(I,L)
          DO 60 J = L1, I
            A(I,J) = A(I,J) - A(I,L) * E(J) - A(J,L) * E(I)
60    CONTINUE
59    CONTINUE

70    CONTINUE

80    D(NM1) = A(NM1,NM1)
      D(N) = A(N,N)
      E(NM1) = A(N,NM1)
C
C  IF EIGENVECTORS ARE REQUIRED, INITIALIZE A.
C
100   IF ( .NOT. VEC ) GO TO 180
C
C  IF THE MATRIX WAS TRIDIAGONAL, SET A EQUAL TO THE IDENTITY MATRIX.
C
      IF ( .NOT. TRD .AND. N .NE. 2 ) GO TO 130

      DO 120 I = 1, N
        DO 110 J = 1, N
          A(I,J) = 0.0E+00
110   CONTINUE
        A(I,I) = 1.0E+00
120   CONTINUE

      GO TO 180
C
C  IF THE MATRIX WAS NOT TRIDIAGONAL, MULTIPLY OUT THE
C  TRANSFORMATIONS OBTAINED IN THE HOUSEHOLDER REDUCTION.C
C
130   A(N,N) = 1.0E+00
      A(NM1,NM1) = 1.0E+00
      A(NM1,N) = 0.0E+00

      DO 170 L = 1, NM2
        LL = NM2 - L + 1
        LL1 = LL + 1

        DO 140 I = LL1, N
          RSUM = 0.0E+00
          DO 135 J = LL1, N
            RSUM = RSUM + A(J,LL) * A(J,I)
135   CONTINUE
          A(LL,I) = RSUM / A(LL,LL)
140   CONTINUE

      DO 149 I  =  LL1, N
            DO 150 J = LL1, N
              A(I,J) = A(I,J) - A(I,LL) * A(LL,J)
150   CONTINUE
149   CONTINUE

          DO 160 I = LL1, N
            A(I,LL) = 0.0E+00
            A(LL,I) = 0.0E+00
160   CONTINUE

          A(LL,LL) = 1.0E+00
170   CONTINUE
C
C  IF AN ABSOLUTE CONVERGENCE CRITERION IS REQUESTED,
C  ( ABSCNV = .TRUE. ), COMPUTE THE INFINITY NORM OF THE MATRIX.
C
180   IF ( .NOT. ABSCNV ) GO TO 200
      NINF = MAX ( 
     &  NINF, ABS ( D(1) ) + ABS ( E(1) ) + ABS ( E(I-1) ) )
      IF ( N .EQ. 2 ) GO TO 200

      DO 190 I = 2, NM1
        NINF = MAX ( NINF, 
     &    ABS ( D(I) ) + ABS ( E(I) ) + ABS ( E(I-1) ) )
190   CONTINUE
C
C  START THE QR ITERATION.
C
200   NU = N
      NUM1 = N - 1
      SHFT = .FALSE.
      K1 = K0
      TEST = NINF * EPS
      E(N) = 0.0E+00
C
C  CHECK FOR CONVERGENCE AND LOCATE THE SUBMATRIX IN WHICH THE
C  QR STEP IS TO BE PERFORMED.
C
210   DO 220 NNL = 1, NUM1
        NL = NUM1 - NNL + 1
        IF ( .NOT. ABSCNV ) 
     &    TEST = EPS * MIN ( ABS ( D(NL) ), ABS ( D(NL+1) ) )
        IF ( ABS ( E(NL) ) .LE. TEST ) GO TO 230
220   CONTINUE

      GO TO 240
230   E(NL) = 0.0E+00
      NL = NL + 1
      IF ( NL .NE. NU ) GO TO 240
      IF ( NUM1 .EQ. 1 ) RETURN
      NU = NUM1
      NUM1 = NU - 1
      GO TO 210

240   E(NU) = E(NU) + 1.0E+00
      IF ( E(NU) .LE. TITTER ) GO TO 250
      FAIL = NU
      RETURN
C
C  CALCULATE THE SHIFT.
C
250   CB = ( D(NUM1) - D(NU) ) / 2.0E+00
      RMAX = MAX ( ABS ( CB ), ABS ( E(NUM1) ) )
      CB = CB / RMAX
      CC = ( E(NUM1) / RMAX )**2
      CD = SQRT ( CB**2 + CC )
      IF ( CB .NE. 0.0E+00 ) CD = SIGN ( CD, CB )
      K2 = D(NU) - RMAX * CC / ( CB + CD )
      IF ( SHFT ) GO TO 270
      IF ( ABS ( K2 - K1 ) .LT. 0.5E+00 * ABS ( K2 ) ) GO TO 260
      K1 = K2
      K = K0
      GO TO 300
260   SHFT = .TRUE.
270   K = K2
C
C  PERFORM ONE QR STEP WITH SHIFT K ON ROWS AND COLUMNS
C  NL THROUGH NU.
C
300   P = D(NL) - K
      Q = E(NL)
      CALL SINCOS ( P, Q, C, S, NORM )

310   DO 380 I = NL, NUM1
C
C  IF REQUIRED, ROTATE THE EIGENVECTORS.
C
        IF ( .NOT. VEC ) GO TO 330

        DO 320 J = 1, N
          TEMP = C * A(J,I) + S * A(J,I+1)
          A(J,I+1) = -S * A(J,I) + C*A(J,I+1)
          A(J,I) = TEMP
320   CONTINUE
C
C  PERFORM THE SIMILARITY TRANSFORMATION AND CALCULATE THE NEXT
C  ROTATION.
C
330     D(I) = C * D(I) + S * E(I)
        TEMP = C * E(I) + S * D(I+1)
        D(I+1) = -S * E(I) + C * D(I+1)
        E(I) = -S * K
        D(I) = C * D(I) + S * TEMP
        IF ( I .EQ. NUM1 ) GO TO 380
        IF ( ABS ( S ) .GT. ABS ( C ) ) GO TO 350
        R = S / C
        D(I+1) = -S * E(I) + C * D(I+1)
        P = D(I+1) - K
        Q = C * E(I+1)
        CALL SINCOS ( P, Q, C, S, NORM )
340     E(I) = R * NORM
        E(I+1) = Q
        GO TO 380
350     P = C * E(I) + S * D(I+1)
        Q = S * E(I+1)
        D(I+1) = C * P / S + K
        E(I+1) = C * E(I+1)
        CALL SINCOS ( P, Q, C, S, NORM )
360     E(I) = NORM

380   CONTINUE

      TEMP = C * E(NUM1) + S * D(NU)
      D(NU) = -S * E(NUM1) + C * D(NU)
      E(NUM1) = TEMP
      GO TO 210

      END
      SUBROUTINE SINCOS ( P, Q, C, S, NORM )
C
C  SINCOS CALCULATES THE ROTATION CORRESPONDING TO THE VECTOR (P,Q).
C
      IMPLICIT NONE

      REAL C
      REAL NORM
      REAL P
      REAL PP
      REAL Q
      REAL QQ
      REAL S

500   PP = ABS ( P )
      QQ = ABS ( Q )
      IF ( QQ .GT. PP ) GO TO 510
      NORM = PP * SQRT ( 1.0E+00 + ( QQ / PP )**2 )
      GO TO 520
510   IF ( QQ .EQ. 0.0E+00 ) GO TO 530
      NORM = QQ * SQRT ( 1.0E+00 + ( PP / QQ )**2 )
520   C = P / NORM
      S = Q / NORM
      RETURN
530   C = 1.0E+00
      S = 0.0E+00
      NORM = 0.0E+00
      RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0