! Test program that solves the generalized eigenvalue problem Ax = eBx
program test
 
  implicit none 
  double precision, allocatable    :: Avals(:), Bvals(:), res(:), E(:), X(:,:)
  integer, allocatable             :: Arows(:), Acols(:), Brows(:), Bcols(:)
  integer                          :: fpm(128), info
  character*1                      :: UPLO
  parameter                        (UPLO='U') ! Store upper parts of matrices
  double precision                 :: epsout, Emin, Emax 
  integer                          :: loop, M, M0
  integer                          :: n, ierr, nnz_A, nnz_B
  integer                          :: i, j
  double precision, allocatable    :: Amat(:,:), Bmat(:,:)

  ! Size of problem
  n = 400
  nnz_A = 4484
  nnz_B = 4484

  ! Matrices
  allocate(Arows(n+1),stat=ierr)
  allocate(Avals(nnz_A),stat=ierr)
  allocate(Acols(nnz_A),stat=ierr)
  
  allocate(Brows(n+1),stat=ierr)
  allocate(Bvals(nnz_B),stat=ierr)
  allocate(Bcols(nnz_B),stat=ierr)
  
  open(1, file = 'Avals.dat', status = 'old')
  do i=1,nnz_A
      read(1,*) Avals(i)
  enddo
  close(1)
  open(2, file = 'Arows.dat', status = 'old')
  do i=1,n+1
      read(2,*) arows(i)
  enddo
  close(2)
  open(3, file = 'Acols.dat', status = 'old')
  do i=1,nnz_A
      read(3,*) acols(i)
  enddo
  close(3)
  open(4, file = 'Bvals.dat', status = 'old')
  do i=1,nnz_B
      read(4,*) bvals(i)
  enddo
  close(4)
  open(5, file = 'Brows.dat', status = 'old')
  do i=1,n+1
      read(5,*) brows(i)
  enddo
  close(5)
  open(10, file = 'Bcols.dat', status = 'old')
  do i=1,nnz_B
      read(10,*) bcols(i)
  enddo
  close(10)


  M0 = 400
  ! For FEAST
  allocate(res(M0),stat=ierr)
  allocate(E(M0),stat=ierr)
  allocate(X(n,n),stat=ierr)
  Emin = -1d0
  Emax = 1d10
  epsout = 0d0

  
  call feastinit(fpm)
  fpm(1)= 1 ! Print run-time output
  fpm(2) =8 
  fpm(27) = 1
  fpm(28) = 1
  write(*,*) 'M0', M0
  write(*,*) 'Emin, Emax', Emin, Emax
  write(*,*) 'n', n

  E = 0d0
  X = 0d0
  res = 0d0

  call dfeast_scsrgv(UPLO,n,Avals,Arows,Acols,Bvals,Brows,Bcols,fpm,epsout,loop,Emin,Emax,M0,E,X,M,res,info)
  print  *,' FEAST OUTPUT INFO ',info
  write(*,*) E
  
end program test











