***QR-METHOD*************************************************************
*     PROGRAM FOR COMPUTING REAL EIGENVALUES OF AN ARBITRARY MATRIX     *
*************************************************************************
*
      program PQR
      implicit doubleprecision (a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd)
*
      print*, 'Enter the order of the square matrix:'
      read(*,*) n
      print*, 'Enter the matrix row-wise:'
      read(*,*) ((A(i,j), j = 1,n), i = 1,n)
      epsil = 1d-6
      maxim = 50
*
      CALL HESSN(A,n)
*
      CALL QR(A,n,maxim,epsil)
*
      write(6,*)'-------------------------------------------------------
     $-----------------------'
      write(6,*) '                     The Eigenvalues of the matrix are
     $:'
      write(6,*)
      write(6,101) (A(i,i),i = 1,n)
      write(6,*)'-------------------------------------------------------
     $-----------------------'
101   format('',29x,f12.4)
*
      stop
      end
*
***FORMATION OF HESSENBERG MATRIX***************************************
*
      subroutine HESSN(A,n)
      implicit doubleprecision (a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),V(nd)
*
      n1  = n - 1
      n2  = n - 2
      np1 = n + 1
      np2 = n + 2
*
      do k = 1, n2
         kk = k + 1
         eta = dabs(A(kk,k))
         do i = kk, n
            if (eta < dabs(A(i,k))) then
               eta = dabs(A(i,k))
            endif
         enddo
         if (eta.ne.0d0) then
            sum = 0
            do i = kk, n
               V(i) = A(i,k) /eta
               A(i,k) = V(i)
               sum = sum + V(i)**2
            enddo
            if (V(kk).ne.0d0) then
               sigma = dsqrt(sum) * V(kk) /dabs(V(kk))
            endif
            V(kk) = V(kk) + sigma
            pi = sigma * V(kk)
            A(np1,k) = pi
*
            do j = k+1, n
               sum = 0
               do i = kk, n
                  sum = sum + V(i)*A(i,j)
               enddo
               rho = sum /pi
               do i = kk, n
                  A(i,j) = A(i,j) - rho*V(i)
               enddo
            enddo
*
            do i = 1, n
               sum = 0
               do j = k+1, n
                  sum = sum + A(i,j)*V(j)
               enddo
               rho = sum /pi
               do j = kk, n
                  A(i,j) = A(i,j) - rho*V(j)
               enddo
            enddo
*
            A(np2,k) = V(kk)
            A(kk,k) = -eta * sigma
         else
            A(np1,k) = 0
         endif
      enddo
*
      do j = 1, n2
      do i = j+2, n
         A(i,j) = 0
      enddo
      enddo
*
      return
      end
*
***QR ALGORITHM (Determination of eigenvalues of Hessenberg matrix)******
*
      subroutine QR(A,n,maxim,epsil)
      implicit doubleprecision (a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),GAMA(nd),SIGMA(nd),kappa
*
      nst = n
      do i = 1, maxim
         kappa = A(n,n)
         A(1,1) = A(1,1) - kappa
         do k = 1, n
            kk = k + 1
            if (k.ne.n) then
               eta = dmax1(dabs(A(k,k)), dabs(A(kk,k)))
               alfa = A(k,k) /eta
               beta = A(kk,k) /eta
               delta = sqrt(alfa**2 + beta**2)
               GAMA(k) = alfa /delta
               SIGMA(k) = beta /delta
               xnu = eta * delta
*
               A(k,k) = xnu
               A(kk,k) = 0
               A(kk,kk) = A(kk,kk) - kappa
*
               do j = kk, n
                  t = A(k,j)
                  A(k,j) = GAMA(k) * A(k,j) + SIGMA(k)* A(kk,j)
                  A(kk,j) = GAMA(k) * A(kk,j) - SIGMA(k)* t
               enddo
            endif
*
            if (k.ne.1) then
               do j = 1, k
                  k1 = k - 1
                  t = A(j,k1)
                  A(j,k1) = GAMA(k1) * A(j,k1) + SIGMA(k1)*A(j,k)
                  A(j,k) = GAMA(k1)*A(j,k) - SIGMA(k1)* t
               enddo
               A(k1,k1) = A(k1,k1) + kappa
            endif
         enddo
         A(n,n) = A(n,n) + kappa
*
         if (dabs(A(n,n-1)) < epsil) then
            n = n - 1
         endif
         if (n == 1) then
            n = nst
            return
         endif
         count = i
         if (count >= maxim) then
            write(*,*)'-------------------------------------------------
     $-------------------------'
            write(*,*) 'The tolerance is not met after 50 iterations'
         endif
      enddo
*
      return
      end