Hi, Dr. Sibert:

I encountered a problem for which I was unable to find
a satisfied data parallel solution. It is a kind of interpolating
points on a curve over a grided plane. I wonder if it is
a problem not good for data parallel implementation.
I'll apprecaite it if you have time to look into it.
Included are the abstracted original problem in f77 and
two solutions I tried in CMfortran. 

Thanks a lot !

--- Gang

c     f77 sequential code --- the abstracted original problem
c     input array: index1, index2, v
c     output array: u
      subroutine interpolation(index1,index2,v,u)
      parameter (N=1000,L=500,M=800)
      parameter (c1=0.15,c2=0.22,c3=0.009,c4=0.91)
      real u(N),v(L,M),index1(N),index2(N)
      integer i,j,k

      do k=1,N
         i = mod(index1(k),L)
         j = mod(index2(k),M)
         u(k) = c1*v(i,j) + c2*v(i+1,j) +
     &          c3*v(i,j+1) + c4*v(i+1,j+1)
      enddo

      write(*,*)(u(k),k=1,N)
      return
      end


c     cmfortran code --- first version, straight forward solution
c     input array: index1, index2, v
c     output array: u
      subroutine interpolation(index1,index2,v,u)
      parameter (N=1000,L=500,M=800)
      parameter (c1=0.15,c2=0.22,c3=0.009,c4=0.91)
      real u(N),v(L,M),index1(N),index2(N)
      integer ii(N),j(N),k,i,j
      real v1(N),v2(N),v3(N),v4(N)
cmf$  layout u(:news),index1(:news),index2(:news)
cmf$  layout v(:news,:news)
cmf$  layout ii(:news),jj(:news)
cmf$  layout v1(:news),v2(:news),v3(:news),v4(:news)

      ii = mod(index1,L)
      jj = mod(index2,M)

c     this do-loop takes 97% of the total execution time
      do k=1,N
         i = ii(k)
         j = jj(k)
         v1(k) = v(i,j)
         v2(k) = v(i+1,j)
         v3(k) = v(i,j+1)
         v4(k) = v(i+1,j+1)
      enddo
         
      u = c1*v1 + c2*v2 + c3*v3 + c4*v4

      write(*,*) u
      return
      end


c     cmfortran code --- second version, much progress
c       but still unsatisfied (slower than the sequential one,
c       even when N=10000)
c     input array: index1, index2, v
c     output array: u
      subroutine interpolation(index1,index2,v,u)
      parameter (N=1000,L=500,M=800)
      parameter (c1=0.15,c2=0.22,c3=0.009,c4=0.91)
      real u(N),v(L,M),index1(N),index2(N)
      integer ii(N),j(N),k,i,j
      real v1,v2,v3,v4
      real v_1d(L*M)
cmf$  layout u(:news),index1(:news),index2(:news)
cmf$  layout v(:news,:news)
cmf$  layout ii(:news),jj(:news)
cmf$  layout v1(:news),v2(:news),v3(:news),v4(:news)
cmf$  layout v_1d(:news)

      forall (i=1:L,j=1:M) v_1d((j-1)*L+i) = v(i,j)
      ii = mod(index1,L)
      jj = mod(index2,M)

      v1 = (jj-1) * L + ii
      v2 = v1 + 1
      v3 = v1 + L
      v4 = v1 + L + 1

c     this step taks 90% of the total execution time
      u = c1*v_1d(v1) + c2*v_1d(v2) + c3*v_1d(v3) + c4*v_1d(v4)

      write(*,*) u
      return
      end
         
