
c correl.f -- A simple one-dimension correlation example
c f77 -O3 -autopar -explicitpar -loopinfo correl.f
c
	implicit real*8(a-h,o-z)
	parameter ( len = 30000, lsearch = 3000, period = 500d0 )
	dimension sig1(len), sig2(len), corrs(lsearch)
c Generate a test case
	center = len/2
	call sinc(sig1, len, 1.0d0, period, center)
	call sinc(sig2, len, 1.0d0, period, center+lsearch/3)
c
c Look for the offset `imx' between sig1 and sig2 at which they are
c maximally correlated.
c
c LOOP M-2. Compute into corrs(i) the correlation between sig1 and sig2
c with sample offset i.  Note that these computations are all
c independent, and can be done in parallel.
c$par doall
	do i = 1, lsearch
		call correl(corrs(i), sig1(1), sig2(1), 1, i, len-lsearch)
	end do
c Find the maximum correlation by looking for max in corrs(i)
	mxcorr = 0
c LOOP M-3. This loop will not parallelize because of dependencies on
c variables mxcorr and imx, but it is a very fast loop.
	do i = 1, lsearch
		if (corrs(i) .gt. mxcorr) then
			imx = i
			mxcorr = corrs(i)
		end if
	end do
	print 10, imx-1
10	format('The maximum offset is ', I5)
	end
c sinc -- fill array sig() with sin(x)/x
c	real*8 sig()				signal array							out
c	integer n				number of samples							in
c	real*8 amp				amplitude							in
c	real*8 period				number of samples							in
c	real*8 center				center of sync function							in
c
	subroutine sinc(sig, n, amp, period, center)
	implicit real*8(a-h,o-z)
	dimension sig(n)
	integer n
c LOOP S-1. This loop will automatically parallelize.
	do i = 1, n
	   a = 2.0 * 3.14159 *(i-center)/period
c Avoid division by zero
	   if (a .ne. 0) then
		  sig(i) = amp * sin(a)/a
	   else
		  sig(i) = amp
	   end if
	end do
	return
	end
c correl -- correlate two signals and return correlation
c	real*8 cor				correlation								out
c	real*8 sig1()				first signal								in
c	real*8 sig2()				second signal								in
c	integer beg1				index of first sig1 sample								in
c	integer beg2				index of first sig2 sample								in
c	integer n				number of samples to look at 								in
c
	subroutine correl(cor, sig1, sig2, beg1, beg2, n)
	implicit real*8(a-h,o-z)
	integer beg1, beg2, n
	dimension sig1(n), sig2(n)
	   sum = 0.0
c LOOP C-1. This loop will not parallelize because of the dependency
c on variable sum from iteration to iteration.
	   do i = 1, n
		  sum = sum + sig1(i+beg1-1) * sig2(i+beg2-1)
	   end do
	   cor = sum
	return
	end
