c  these subroutines are required for the DeFT program.  they were
c  previously public domain.  they also include domain routines that
c  replace existing functions on cray machines
c
      subroutine flush(iunit)

      return
      end
c***********************************************************************
      SUBROUTINE MATINV(A,B,N,L,IDIM)

      IMPLICIT REAL*8(A-H,O-Z)

      PARAMETER (NDMATI=1000)

      DIMENSION A(IDIM,IDIM),B(IDIM),IP(NDMATI),IN(NDMATI,2)

      DATA ZERO,ONE,SMALL/0.0D0,1.0D0,1.0D-30/

      IF(N.GT.NDMATI) stop 'increase size of tables in matinv'

      D=ONE
      DO 1002 I=1,N
 1002 IP(I)=0

      DO 1003 I=1,N

      AMAX=ZERO

      DO 1004 J=1,N

      IF(IP(J).GT.0) GOTO 1004
      IF(IP(J).LT.0) GOTO 1005

      DO 1006 K=1,N

      IF(IP(K).EQ.1) GOTO 1006
      IF(IP(K).GT.1) GOTO 1005

      IF(ABS(A(J,K)).LE.AMAX) GOTO 1006
      IR=J
      IC=K
      AMAX=ABS(A(J,K))

 1006 CONTINUE

 1004 CONTINUE

      IP(IC)=IP(IC)+1

      IF(AMAX.GT.SMALL) GOTO 1007

 1005 STOP 'MATINV'

 1007 IF(IR.EQ.IC) GOTO 1008

      D=-D

      DO 1009 K=1,N
      AMAX=A(IR,K)
      A(IR,K)=A(IC,K)
 1009 A(IC,K)=AMAX

      IF(L.EQ.0) GOTO 1008

      AMAX=B(IR)
      B(IR)=B(IC)
      B(IC)=AMAX

 1008 IN(I,1)=IR

      IN(I,2)=IC
      AMAX=A(IC,IC)
      D=D*AMAX
      A(IC,IC)=ONE

      DO 1010 K=1,N
 1010 A(IC,K)=A(IC,K)/AMAX

      IF(L.EQ.0) GOTO 1011

      B(IC)=B(IC)/AMAX

 1011 DO 1003 J=1,N

      IF(J.EQ.IC) GOTO 1003

      AMAX=A(J,IC)
      A(J,IC)=ZERO

      DO 1012 K=1,N
 1012 A(J,K)=A(J,K)-A(IC,K)*AMAX

      IF(L.EQ.0) GOTO 1003
      B(J)=B(J)-B(IC)*AMAX

 1003 CONTINUE

      IF(L.EQ.1) GOTO 1013

      DO 1014 I=1,N

      J=N+1-I

      IF(IN(J,1).EQ.IN(J,2)) GOTO 1014
      IR=IN(J,1)
      IC=IN(J,2)

      DO 1015 K=1,N
      AMAX=A(K,IR)
      A(K,IR)=A(K,IC)
 1015 A(K,IC)=AMAX

 1014 CONTINUE
 1013 CONTINUE

      return
      end
C***********************************************************************
      SUBROUTINE TQL9(NM,N,D,E,Z,AUX,AUX1,IERR)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION D(N),E(N),Z(NM,N),AUX(N),AUX1(N)
C
      DATA PRECSN,ZERO,ONE,TWO/1.0D-24,0.0D0,1.0D0,2.0D0/
C
      IERR=0
      IF(N.EQ.1) GOTO 1600
      DO 100 I=2,N
  100 E(I-1)=E(I)
      F=ZERO
      B=ZERO
      E(N)=ZERO
      DO 1100 L=1,N
      J=0
      H=PRECSN*(ABS(D(L))+ABS(E(L)))
      IF(B.LT.H) B=H
      DO 200 M=L,N
      IF(ABS(E(M)).LE.B) GOTO 300
  200 CONTINUE
  300 IF(M.EQ.L) GOTO 1000
  400 IF(J.EQ.30) GOTO 1500
      J=J+1
      L1=L+1
      G=D(L)
      P=(D(L1)-G)/(TWO*E(L))
      R=SQRT(P*P+ONE)
      D(L)=E(L)/(P+SIGN(R,P))
      H=G-D(L)
      DO 500 I=L1,N
  500 D(I)=D(I)-H
      F=F+H
      P=D(M)
      C=ONE
      S=ZERO
      MML=M-L
      DO 900 II=1,MML
      I=M-II
      G=C*E(I)
      H=C*P
      IF(ABS(P).LT.ABS(E(I))) GOTO 600
      C=E(I)/P
      R=SQRT(C*C+ONE)
      E(I+1)=S*P*R
      S=C/R
      C=ONE/R
      GOTO 700
  600 C=P/E(I)
      R=SQRT(C*C+ONE)
      E(I+1)=S*E(I)*R
      S=ONE/R
      C=C*S
  700 P=C*D(I)-S*G
      D(I+1)=H+S*(C*G+S*D(I))
      DO 800 K=1,N
      AUX(K)=C*Z(K,I+1)
      AUX1(K)=S*Z(K,I+1)
      Z(K,I+1)=S*Z(K,I)
  800 Z(K,I)=C*Z(K,I)
      DO 850 K=1,N
      Z(K,I+1)=Z(K,I+1)+AUX(K)
  850 Z(K,I)=Z(K,I)-AUX1(K)
  900 CONTINUE
      E(L)=S*P
      D(L)=C*P
      IF(ABS(E(L)).GT.B) GOTO 400
 1000 D(L)=D(L)+F
 1100 CONTINUE
      DO 1400 II=2,N
      I=II-1
      K=I
      P=D(I)
      DO 1200 J=II,N
      IF(D(J).GE.P) GOTO 1200
      K=J
      P=D(J)
 1200 CONTINUE
      IF(K.EQ.I) GOTO 1400
      D(K)=D(I)
      D(I)=P
      DO 1250 J=1,N
      AUX(J)=Z(J,I)
 1250 Z(J,I)=Z(J,K)
      DO 1300 J=1,N
 1300 Z(J,K)=AUX(J)
 1400 CONTINUE
      GOTO 1600
 1500 IERR=L
 1600 RETURN
      END
C***********************************************************************
      SUBROUTINE TRED9(NM,N,A,D,E,AUX,Z)
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      DIMENSION A(NM,N),D(N),E(N),Z(NM,N),AUX(N)
C
      DATA ZERO,ONE/0.0D0,1.0D0/
C
      DO 100 I=1,N
      DO 100 J=1,I
      Z(I,J)=A(I,J)
  100 CONTINUE
      IF(N.EQ.1) GOTO 1300
      DO 1200 II=2,N
      I=N+2-II
      L=I-1
      H=ZERO
      SCALE=ZERO
      IF(L.LT.2) GOTO 300
      DO 200 K=1,L
      AUX(K)=ABS(Z(I,K))
  200 SCALE=SCALE+AUX(K)
      IF(SCALE.NE.ZERO) GOTO 400
  300 E(I)=Z(I,L)
      GOTO 1100
  400 DO 500 K=1,L
      Z(I,K)=Z(I,K)/SCALE
      H=H+Z(I,K)*Z(I,K)
  500 CONTINUE
      F=Z(I,L)
      G=-SIGN(SQRT(H),F)
      E(I)=SCALE*G
      H=H-F*G
      Z(I,L)=F-G
      F=ZERO
      DO 900 J=1,L
      Z(J,I)=Z(I,J)/H
      G=ZERO
      DO 600 K=1,J
  600 G=G+Z(J,K)*Z(I,K)
      JP1=J+1
      IF(L.LT.JP1) GOTO 800
      DO 700 K=JP1,L
  700 G=G+Z(K,J)*Z(I,K)
  800 E(J)=G/H
      F=F+E(J)*Z(I,J)
  900 CONTINUE
      HH=F/(H+H)
      DO 1000 J=1,L
      F=Z(I,J)
      AUX(J)=E(J)-HH*F
      E(J)=AUX(J)
      DO 950 K=1,J
      AUX(K)=F*E(K)
  950 Z(J,K)=Z(J,K)-AUX(K)
      DO 1000 K=1,J
      AUX(K)=E(J)*Z(I,K)
 1000 Z(J,K)=Z(J,K)-AUX(K)
 1100 D(I)=H
 1200 CONTINUE
 1300 D(1)=ZERO
      E(1)=ZERO
      DO 1800 I=1,N
      L=I-1
      IF(D(I).EQ.ZERO) GOTO 1600
      DO 1500 J=1,L
      G=ZERO
      DO 1400 K=1,L
 1400 G=G+Z(I,K)*Z(K,J)
      DO 1450 K=1,L
 1450 AUX(K)=G*Z(K,I)
      DO 1500 K=1,L
 1500 Z(K,J)=Z(K,J)-AUX(K)
 1600 D(I)=Z(I,I)
      Z(I,I)=ONE
      IF(L.LT.1) GOTO 1800
      DO 1650 J=1,L
 1650 Z(I,J)=ZERO
      DO 1700 J=1,L
      Z(J,I)=ZERO
 1700 CONTINUE
 1800 CONTINUE
      RETURN
      END
