	PROGRAM polint
	USE types

        integer, parameter :: ngx=15  ! 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

!HPF$ DISTRIBUTE grid(BLOCK)
!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


forall(i=1:ngx-n:2) x(i)=coord(i+1)
forall(i=1:ngx-2:2) xa(i,1:n)=coord(i:i+2*n-1:2)
forall(i=1:ngx-2:2) ya(i,1:n)=grid(i:i+2*n-1:2)
forall(i=1:ngx-2:2) c(i,:)=ya(i,:)
forall(i=1:ngx-2:2) d(i,:)=ya(i,:)
forall(i=1:ngx-2:2) ho(i,:)=xa(i,:)-x(i)

ns=1
forall(i=1:ngx-2:2) y(i)=ya(i,ns(i))
ns=ns-1


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

do i=1,ngx-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
