	PROGRAM polint
	USE types

        integer, parameter :: ngx=20  ! size of the grid
        integer, parameter :: n=3     ! order of interpolation
	
        REAL(SP), DIMENSION(ngx) :: grid,coord 
        REAL(SP), DIMENSION(ngx) :: x,y,dy,nx
        INTEGER(i4B), DIMENSION(ngx) :: ns
        REAL(SP), DIMENSION(ngx,n) :: xa,ya
	REAL(SP), DIMENSION(ngx,n) :: c,d,ho,den
	INTEGER(I4B) :: m,i,k

!HPF$ PROCESSORS PROC(4)
!HPF$ DISTRIBUTE grid(BLOCK) ONTO PROC
!HPF$ ALIGN (:) WITH grid(:) :: coord,x,y,dy,nx,ns
!HPF$ ALIGN (:,*) WITH grid(:) :: xa,ya,c,d,ho,den


!       initialize coordinates and function values 
        forall(i=1:ngx) coord(i)=i-1
        forall(i=1:ngx:2) grid(i)=3*coord(i)+1 !initalize every second
!                                          the rest is to be inerpolated

write(6,*) coord
write(6,*) '=============================='
forall(i=1:ngx-n:2) x(i)=coord(i+1)
do k=1,n
forall(i=1:ngx-2*n:2) xa(i,k)=coord(i+2*(k-1))
write(6,*) xa(:,k)
write(6,*) '=============================='
forall(i=1:ngx-2*n:2) ya(i,k)=grid(i+2*(k-1))
forall(i=1:ngx:2) c(i,k)=ya(i,k)
forall(i=1:ngx:2) d(i,k)=ya(i,k)
forall(i=1:ngx-2*k:2) ho(i,k)=xa(i,k)-x(i)
enddo
ns=1
forall(i=1:ngx-n:2) y(i)=ya(i,ns(i))
ns=ns-1

do m=1,n-1
do k=1,n-m
forall(i=1:ngx-2*n:2) den(i,k)=ho(i,k)-ho(i,k+m)
forall(i=1:ngx-2*n:2) den(i,k)=(c(i,k+1)-d(i,k))/den(i,k)
forall(i=1:ngx-2*n:2) d(i,k)=ho(i,k+m)*den(i,k)
forall(i=1:ngx-2*n:2) c(i,k)=ho(i,k)*den(i,k)
enddo
forall(i=1:ngx-2*n:2) dy(i)=c(i,ns(i)+1)
forall(i=1:ngx-2*n:2) y(i)=y(i)+dy(i)
enddo

do i=1,ngx-2*n,2
write(6,*) 'interpolation polynomial'
write(6,*) xa(i,:)
write(6,*) ya(i,:)
write(6,*) 'x-value',x(i)
write(6,*) 'solution',y(i)

enddo

	END program polint
