program main 
! Jan Mandel September 2008
! Jan.Mandel@gmail.com

! a simple main program for cuda testing
! orignally based on Fotran_Cuda from nvidia.com

implicit none 

!define length of the array
!integer, parameter :: N=4, M=5, S=6, bw=10
!integer, parameter :: N=1000000, M=1000000, bw=100, S=bw*N
integer, parameter :: N=200000, M=200000, bw=100, S=bw*N
!integer, parameter :: N=4, M=5, bw=2, S=bw*N
!integer, parameter :: N=40000, M=50000, bw=200, S=bw*N
!integer, parameter:: pres=1
integer, parameter:: pres=0

integer*8 t1,t2,time_cuda,time_fort
integer,dimension(N+1)::x
integer,dimension(S)::r
real,dimension(S)::v
real,dimension(M)::a
real,dimension(N)::b,b_gold,b_cuda
real::t_gold,t_cuda
!data x/1,3,6,6,7/
!data r/1,2,2,3,1,5/

integer::i,j,k


x(1)=1
do i=1,N
x(i+1)=x(i)+bw
 do k=x(i),x(i+1)-1
  j=i+k-x(i)
  if (j>0.and. j<M)then
   r(k)=j
  else
   r(k)=0
  endif
 enddo
enddo

call cuda_device_init(0) ! print device name and make sure cupports cuda
print *,'Sparse matrix-vector multiplication'
print *,'Matrix size ',M,N,' nonzeros',x(N+1) 
print *,'Time in ms'
  
do i = 1, M
 a(i)=mod(i,1000)/1000.-.5
enddo
do i=1,S
 v(i)=mod(i,1000)/1000.-.5
enddo
 
call time(t1)
call spmtv(M,N,S,x,r,v,a,b)
call time(t2)
if(pres.gt.0) print *,'result from fortran:',b
print *,'time from fortran:',t2-t1
time_fort=t2-t1

call time(t1)
call spmtv_gold(M,N,S,x,r,v,a,b_gold)
call time(t2)
if(pres.gt.0) print *,'result from C      :',b_gold
print *,'time from C      :',t2-t1

do i=1,2
call time(t1)
call spmtv_cuda(1,M,N,S,x,r,v,a,b_cuda)
call time(t2)
print *,'time matrix in   :',t2-t1
call time(t1)
call spmtv_cuda(2,M,N,S,x,r,v,a,b_cuda)
call time(t2)
print *,'time vector in   :',t2-t1
time_cuda=t2-t1
call time(t1)
call spmtv_cuda(3,M,N,S,x,r,v,a,b_cuda)
call time(t2)
print *,'time from cuda   :',t2-t1
time_cuda=time_cuda+t2-t1
call time(t1)
call spmtv_cuda(4,M,N,S,x,r,v,a,b_cuda)
call time(t2)
print *,'time vector out  :',t2-t1
time_cuda=time_cuda+t2-t1
print *,'repeat cuda time :',time_cuda
print *,'cuda speedup:',time_fort/real(time_cuda)
if(pres.gt.0) print *,'result from cuda   :',b_cuda
call spmtv_cuda(5,M,N,S,x,r,v,a,b_cuda) ! free device memory
enddo
t_gold=0.
t_cuda=0.
do i=1,N
 t_gold=max(t_gold,abs(b(i)-b_gold(i)))
 t_cuda=max(t_cuda,abs(b(i)-b_cuda(i)))
enddo
print *,'difference C:',t_gold,' cuda: ',t_cuda

end program main

subroutine time(t)
integer*8 t
integer*8, dimension(8):: v
!integer count,rate,cmax
!call system_clock(count,rate,cmax)
!t=dble(count)/dble(rate)
CALL DATE_AND_TIME(VALUES=v) 
t=v(8)+1000*(v(7)+60*(v(6)+12*v(5)))
!print *,v(5),':',v(6),':',v(7),',',v(8)
end