program diag
c ****************** Program DIAG ******************
c This program is designed to read in a real
c symmetric matrix (but as a square matrix on
c disk), diagonalize it, and return all
c eigenvalues and corresponding eigenvectors.
c
c Copyright (C) 1992 Jeff Nichols and Jack Simons
c
c **************************************************
      implicit real (a-h,o-z), integer(i-n)
      include "limits.h"
c
c     Set core needed to:
c         1 square array (eigenvectors) of order mmax
c         1 symmetric array (input) of order mmax
c         1 vector (eigenvalues) of length mmax
c         1 scratch vector of length mmax
c     total = mmax*mmax + mmax*(mmax+1)/2 + 2*mmax
c
      dimension pcore(mmaxsq + mmaxtr + 2*mmax), ipntr(10)
      itotal = mmax*mmax + mmax*(mmax+1)/2 + 2*mmax + 10
      call modinfo('diag        ')
      call shwmem('diag        ',itotal)
      call goon
   10 continue
      write(*,*)
      write(*,*)' Input the order of the real symmetric '
      write(*,*)' array you would like diagonalized. '
      read(*,*)msize
      if(msize.gt.mmax)then
        write(*,1010)mmax
        goto 10
      endif
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + (msize*(msize+1))/2
      ipntr(3) = ipntr(2) + msize*msize
      ipntr(4) = ipntr(3) + msize
      ipntr(5) = ipntr(4) + msize
      call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &            pcore(ipntr(4)),msize)
      call qmexit
      stop
 1010 format(' Maximum order possible is less than ',I3,/,
     &       ' Please try again. ')
      end
      subroutine driver(asym,avects,avals,scr1,msize)
      implicit real (a-h,o-z), integer(i-n)
      logical yesno
      character*30 dfname, fname, dscrpt
      character*1 ans
      dimension asym((msize*(msize+1))/2),
     &          avects(msize,msize),
     &          avals(msize),
     &          scr1(msize)
   10 continue
      dscrpt = 'matrix to be diagonalized'
      ids = 25
      dfname = 'diag.in'
      idf = 7
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,avects,msize,msize,
     &              iflnum)
c
c       Close input file.
c
        close(unit=32)
      else
        call rdmatu(dscrpt,ids,avects,msize,msize)
      endif
c
c     Write input array back to user.
c
      call wtmatu(dscrpt,ids,avects,msize,msize)
      write(*,*)' Is this correct? <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
c
c       Pack the square array into a 1-dimensional symmetric
c       array.
c
        ij = 0
        do 40 i = 1,msize
          do 30 j = 1,i
            ij = ij + 1
            asym(ij) = avects(i,j)
            avects(i,j) = 0.0e0
            avects(j,i) = 0.0e0
   30     continue
   40   continue
c
c       Diagonalize the array.
c
        call diaghe(asym,msize,avals,msize,avects,msize,scr1)
c
c       Write the eigenvalues to user.
c
        write(*,*)' The eigenvalues of this matrix are: '
        call matout(avals,msize,1,msize)
c
c       Write the eigenvectors to user.
c
        write(*,*)' their corresponding eigenvectors are: '
        call matout(avects,msize,msize,msize)
      else
        write(*,*)' Please try again. '
        goto 10
      endif
      dscrpt = 'eigenvectors'
      ids = 12
      dfname = 'diag_vectors.out'
      idf = 16
      fname = '                              '
      ifn = 30
      iflnum = 33
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,avects,msize,msize,
     &              iflnum)
c
c       Close output file.
c
        close(unit=33)
      endif
      dscrpt = 'eigenvalues'
      ids = 11
      dfname = 'diag_values.out'
      idf = 15
      fname = '                              '
      ifn = 30
      iflnum = 34
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,avals,msize,1,
     &              iflnum)
c
c       Close output file.
c
        close(unit=34)
      endif
      return
      end