program eigendemo
C
C  demonstrates eigenvalues
C
C  CWJ  SDSU March 2005
C
C
      implicit none

      integer size
      parameter(size=100)
      
      integer n
      integer i

      real a(size,size)  ! the matrix in question
      real e(size)       ! eigenvalues
      real work(size)
      real vec(size,size) ! eigenvectors

      character ychar*1
      integer nvec
      real vec1(size),vec2(size)

      print*,' enter dimension of matrix '
      read*,n

      call fillmatrix(size,n,a)
      print*,' '
      print*,' Original matrix : '
      call printnicematrix(size,n,a)
      
      call eig(a,n,size,e,vec,work)
      print*,' '
      print*,' eigenvalues: '
      do i =1,n
          write(6,101)e(i)
      enddo
  101 format(f10.5)
      print*,' '
      print*,' Notice not sorted; do you want to sort?'
      read*,ychar
      if(ychar.eq.'y')then
        call eigsrt(e,vec,n,size)
        do i =1,n
          write(6,101)e(i)
        enddo      
      endif
      print*,' '
      print*,' now to check eigenvectors '
      
      do while(ychar.eq.'y')
          print*,' Enter which eigenvector to check '
          read*,nvec
          print*,' eigenvalue = ',e(nvec)
          print*,' '
          call copyvec(n,size,vec,nvec,vec1)
          call multiplyvec(n,size,a,vec1,vec2)
          print*,' '
          print*,' Column 1 = original vector '
          print*,' Column 2 = original vector x eigenvalue '
          print*,' Column 3 = A x vector '
          do i = 1,n
            write(6,102)vec1(i),vec1(i)*e(nvec),vec2(i)
          enddo
  102  format(3f10.5)
          print*,' '
          print*,' Do you want to check another eigenvector?'
          read*,ychar
      enddo 

      end

C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine fillmatrix(size,n,a)
C
C  fills a matrix with random elements
C  makes it symmetric
C   
C  calls function GAUSSVAR (from RANLIB)

      implicit none
      integer size
      real a(size,size)
      integer n
      integer i,j
      integer seed
     

      real gaussvar
      print*,' give me an integer seed '
      read*,seed
      do i = 1,n
	do j = i,n
           a(i,j)=gaussvar(seed)
           a(j,i)=a(i,j)
        enddo
      enddo
      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine printnicematrix(size,n,array)

      implicit none

      integer size,n
      real array(size,size)
      integer i,j

      do i =1,n
           write(6,100)(array(i,j),j=1,n)
       enddo
  100 format(6f10.5)
       return
       end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       subroutine copyvec(n,np,array,i,vec)
       implicit none

       integer np
       real array(np,np),vec(np)
       integer n
       integer i,j
      
       do j =1,n
          vec(j) = array(j,i)
   
       enddo
       return
      end
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCccc

      subroutine multiplyvec(n,np,a,v1,v2)
      implicit none
      integer np,n
      real a(np,np)
      real v1(np),v2(np)
      integer i,j

      do i = 1,n
         v2(i) = 0.0
         do j = 1,n
            v2(i) = v2(i)+a(i,j)*v1(j)
         enddo
      enddo
      return
      end