      SUBROUTINE HQRII(A,N,M,E,V)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension A(*), E(N), V(N,M)
*************************************************************
*
* HQRII IS A DIAGONALISATION ROUTINE, WRITTEN BY YOSHITAKA BEPPU OF
*       NAGOYA UNIVERSITY, JAPAN.
*       FOR DETAILS SEE 'COMPUTERS & CHEMISTRY' VOL.6 1982. PAGE 000.
*
* ON INPUT    A       = MATRIX TO BE DIAGONALISED (PACKED CANONICAL)
*             N       = SIZE OF MATRIX TO BE DIAGONALISED.
*             M       = NUMBER OF EIGENVECTORS NEEDED.
*             E       = ARRAY OF SIZE AT LEAST N
*             V       = ARRAY OF SIZE AT LEAST NMAX*M
*
* ON OUTPUT   E       = EIGENVALUES
*             V       = EIGENVECTORS IN ARRAY OF SIZE NMAX*M
*
************************************************************************
      dimension W(5,MAXPAR)
      if (N.LE.1 .OR. M .LE.1 .OR. M .GT. N) then
         if (N.EQ.1 .AND. M.EQ.1) then
            E(1)=A(1)
            V(1,1)=1.D0
            return
         end if
         write(MIOB,'(////10X,''IN HQRII, N ='',I4,'' M ='',I4)') N,M
         call filclo ! J/M
         stop ' in HQRII - nonsense in size of matrix'
      end if
*
* EPS3 AND EPS ARE MACHINE-PRECISION DEPENDENT
*
      EPS3=1.D-30
      ZERO=0.D0
      LL=(N*(N+1))/2+1
      EPS=1.D-8
      IORD=-1
      NM1=N-1
      if (N.EQ.2) go to 90
      NM2=N-2
      KRANK=0
C     HOUSEHOLDER TRANSFORMATION
      do 80 K=1,NM2
         KP1=K+1
         KRANK=KRANK+K
         W(2,K)=A(KRANK)
         SUM=0.
         JRANK=KRANK
         do 10 J=KP1,N
            W(2,J)=A(JRANK+K)
            JRANK=JRANK+J
   10    SUM=W(2,J)**2+SUM
         S=SIGN(SQRT(SUM),W(2,KP1))
         W(1,K)=-S
         W(2,KP1)=W(2,KP1)+S
         A(K+KRANK)=W(2,KP1)
         H=W(2,KP1)*S
         if (ABS(H).LT.EPS3) go to 80
         SUMM=0.D0
         IRANK=KRANK
         do 50 I=KP1,N
            SUM=0.D0
            do 20 J=KP1,I
   20       SUM=SUM+A(J+IRANK)*W(2,J)
            if (I.GE.N) go to 40
            IP1=I+1
            JRANK=I*(I+3)/2
            do 30 J=IP1,N
               SUM=SUM+A(JRANK)*W(2,J)
   30       JRANK=JRANK+J
   40       W(1,I)=SUM/H
            IRANK=IRANK+I
   50    SUMM=W(1,I)*W(2,I)+SUMM
         U=SUMM*0.5D0/H
         JRANK=KRANK
         do 70 J=KP1,N
            W(1,J)=W(2,J)*U-W(1,J)
            do 60 I=KP1,J
   60       A(I+JRANK)=W(1,I)*W(2,J)+W(1,J)*W(2,I)+A(I+JRANK)
   70    JRANK=JRANK+J
   80 A(KRANK)=H
   90 W(2,NM1)=A((NM1*(NM1+1))/2)
      W(2,N)=A((N*(N+1))/2)
      W(1,NM1)=A(NM1+(N*(N-1))/2)
      W(1,N)=0.D0
      GERSCH=ABS(W(2,1))+ABS(W(1,1))
      do 100 I=1,NM1
  100 GERSCH=MAX(ABS(W(2,I+1))+ABS(W(1,I))+ABS(W(1,I+1)),GERSCH)
      DEL=EPS*GERSCH
      do 110 I=1,N
         W(3,I)=W(1,I)
         E(I)=W(2,I)
  110 V(I,M)=E(I)
      if (ABS(DEL).LT.EPS3)  go to  220
C     QR-METHOD WITH ORIGIN SHIFT
      K=N
  120 L=K
  130 if (ABS(W(3,L-1)).LT.DEL) go to 140
      L=L-1
      if (L.GT.1)  go to 130
  140 if (L.EQ.K)  go to 170
      WW=(E(K-1)+E(K))*0.5D0
      R=E(K)-WW
      Z=SIGN(SQRT(W(3,K-1)**2+R*R),R)+WW
      EE=E(L)-Z
      E(L)=EE
      FF=W(3,L)
      R=SQRT(EE*EE+FF*FF)
      J=L
      go to 160
  150 R=SQRT(E(J)**2+W(3,J)**2)
      W(3,J-1)=S*R
      EE=E(J)*C
      FF=W(3,J)*C
  160 C=E(J)/R
      S=W(3,J)/R
      WW=E(J+1)-Z
      E(J)=(FF*C+WW*S)*S+EE+Z
      E(J+1)=C*WW-S*FF
      J=J+1
      if (J.LT.K) go to 150
      W(3,K-1)=E(K)*S
      E(K)=E(K)*C+Z
      go to 120
  170 K=K-1
      if (K.GT.1) go to 120
*    *    *    *    *    *    *    *    *    *    *    *    *
*
*   AT THIS POINT THE ARRAY 'E' CONTAINS THE UN-ORDERED EIGENVALUES
*
*    *    *    *    *    *    *    *    *    *    *    *    *
C     STRAIGHT SELECTION SORT OF EIGENVALUES
      SORTER=1.D0
      if (IORD.LT.0) SORTER=-1.D0
      J=N
  180 L=1
      II=1
      LL=1
      do 200 I=2,J
         if ((E(I)-E(L))*SORTER .GT. 0.D0) go to 190
         L=I
         go to 200
  190    II=I
         LL=L
  200 continue
      if (II.EQ.LL) go to 210
      WW=E(LL)
      E(LL)=E(II)
      E(II)=WW
  210 J=II-1
      if (J.GE.2) go to 180
  220 if (M.EQ.0) return
***************
*  ORDERING OF EIGENVALUES COMPLETE.
***************
C      INVERSE-ITERATION FOR EIGENVECTORS
      FN=FLOAT(N)
      EPS1=1.D-5
      SEPS=SQRT(EPS)
      EPS2=0.05D0
      RN=0.D0
      RA=EPS*0.6180339887485D0
C    0.618... IS THE FIBONACCI NUMBER (-1+SQRT(5))/2.
      IG=1
      do 450 I=1,M
         IM1=I-1
         do 230 J=1,N
            W(3,J)=0.D0
            W(4,J)=W(1,J)
            W(5,J)=V(J,M)-E(I)
            RN=RN+RA
            if (RN.GE.EPS) RN=RN-EPS
  230    V(J,I)=RN
         do 260 J=1,NM1
            if (ABS(W(5,J)).GE.ABS(W(1,J))) go to 240
            W(2,J)=-W(5,J)/W(1,J)
            W(5,J)=W(1,J)
            T=W(5,J+1)
            W(5,J+1)=W(4,J)
            W(4,J)=T
            W(3,J)=W(4,J+1)
            if (ABS(W(3,J)).LT.EPS3) W(3,J)=DEL
            W(4,J+1)=0.D0
            go to 250
  240       if (ABS(W(5,J)).LT.EPS3) W(5,J)=DEL
            W(2,J)=-W(1,J)/W(5,J)
  250       W(4,J+1)=W(3,J)*W(2,J)+W(4,J+1)
  260    W(5,J+1)=W(4,J)*W(2,J)+W(5,J+1)
         if (ABS(W(5,N)) .LT. EPS3) W(5,N)=DEL
         do 320 ITERE=1,5
            if (ITERE.EQ.1) go to 280
            do 270 J=1,NM1
               if (ABS(W(3,J)).LT.EPS3) go to 270
               T=V(J,I)
               V(J,I)=V(J+1,I)
               V(J+1,I)=T
  270       V(J+1,I)=V(J,I)*W(2,J)+V(J+1,I)
  280       V(N,I)=V(N,I)/W(5,N)
            V(NM1,I)=(V(NM1,I)-V(N,I)*W(4,NM1))/W(5,NM1)
            VN=MAX(ABS(V(N,I)),ABS(V(NM1,I)),1.D-20)
            if (N.EQ.2) go to 300
            K=NM2
  290       V(K,I)=(V(K,I)-V(K+1,I)*W(4,K)-V(K+2,I)*W(3,K))/W(5,K)
            VN=MAX(ABS(V(K,I)),VN,1.D-20)
            K=K-1
            if (K.GE.1) go to 290
  300       S=EPS1/VN
            do 310 J=1,N
  310       V(J,I)=V(J,I)*S
            if (ITERE.GT.1 .AND. VN.GT.1) go to 330
  320    continue
C     TRANSFORMATION OF EIGENVECTORS
  330    if (N.EQ.2) go to 380
         KRANK=NM2*(N+1)/2
         KPIV=NM2*NM1/2
         do 370 K=NM2,1,-1
            KP1=K+1
            if (ABS(A(KPIV)).LE.EPS3) go to 360
            SUM=0.D0
            do 340 KK=KP1,N
               SUM=SUM+A(KRANK)*V(KK,I)
  340       KRANK=KRANK+KK
            S=-SUM/A(KPIV)
            do 350 KK=N,KP1,-1
               KRANK=KRANK-KK
  350       V(KK,I)=A(KRANK)*S+V(KK,I)
  360       KPIV=KPIV-K
  370    KRANK=KRANK-KP1
  380    do 390 J=IG,I
            if (ABS(E(J)-E(I)) .LT. EPS2) go to 400
  390    continue
         J=I
  400    IG=J
         if (IG .EQ. I) go to 430
C     RE-ORTHOGONALISATION
         do 420 K=IG,IM1
            SUM=0.D0
            do 410 J=1,N
  410       SUM=V(J,K)*V(J,I)+SUM
            S=-SUM
            do 420 J=1,N
  420    V(J,I)=V(J,K)*S+V(J,I)
C     NORMALISATION
  430    SUM=1.D-24
         do 440 J=1,N
  440    SUM=SUM+V(J,I)**2
         SINV=1.D0/SQRT(SUM)
         do 450 J=1,N
  450 V(J,I)=V(J,I)*SINV
      return
      end
