program trsg

c----------------------------------------------------------------------
c Program to test the subroutine RSG
c----------------------------------------------------------------------

c----------------------------------------------------------------------
c necessary variables
c----------------------------------------------------------------------
      parameter (nm=5)

      real*8 pi 
      real*8 a(0:nm,0:nm),b(0:nm,0:nm),eigenvectors(0:nm,0:nm)
      real*8 eigenvalues(0:nm),fv1(0:nm),fv2(0:nm)
      integer ierr

      pi = 4d0*atan(1d0)
c----------------------------------------------------------------------
c fill the matrices a and b
c----------------------------------------------------------------------
      
      do i=0,nm
         do j=0,nm
c b is just the unit matrix, i.e. solve an ordinary eigenvalue problem
            if (i .eq. j) then
               b(i,j) = 1d0
            else
               b(i,j) = 0d0
            endif
c a has on the diagonal everywhere 2's and on the off-diagonals 
c everywhere -1's
            if (i .eq. j) then
               a(i,j) = 2d0
            else
               if ( (i .eq. (j-1)) .or. (i .eq. (j+1)) )  then
                  a(i,j) = -1d0
               else
                  a(i,j) = 0d0   
               endif
            endif
         enddo
      enddo


c----------------------------------------------------------------------
c write matrices a and b on screen
c----------------------------------------------------------------------
      write (6,*) ''
      write (6,*) 'Matrix a'
      do j=0,nm
         write(6,4015) (a(i,j),i=0,nm)
      enddo

      write (6,*) ''
      write (6,*) 'Matrix b'
      do j=0,nm
         write(6,4015) (b(i,j),i=0,nm)
      enddo


      call RSG(nm+1,nm+1,a,b,eigenvalues,1,eigenvectors,fv1,fv2,ierr)

c----------------------------------------------------------------------
c write eigenvectors on screen
c----------------------------------------------------------------------

      write (6,*) ''
      write (6,*) 'Eigenvectors'
      do j=0,nm
         write(6,4015) (eigenvectors(i,j),i=0,nm)
      enddo


c----------------------------------------------------------------------
c write eigenvalues on screen and compare with analytic values
c----------------------------------------------------------------------

      write (6,*) ''
      write (6,*) 'Eigenvalues'
      write(6,4015) (eigenvalues(i),i=0,nm)

      write (6,*) ''
      write (6,*) 'analytic values'
      write (6,4015) (4d0*sin(i*pi/(2d0*(nm+2)) )**2,i=1,nm+1)
         

 4015 format(/,(20F10.6))


      stop
      end