3D convolution

[Jayashree_N_]

I am using MKL's VSL multi-dimensional convolution routines to perform a 3D convolution. As a test, I am using the input function as a 3D gaussian as follows:

f(x,y,z) = exp(- (x-2)**2) exp(-(y+2)**2) exp(- (z**2))

and the convolution kernel as

g(x,y,z) = exp(-(x**2 + y**2 + z**2))

The resulting function is formally C(x,y,z), also a gaussian. I have some trouble understanding the access of the exact indices of the ouput 3D array. Given below is an output for the non-zero values of the convoluted function C, in the last column, calculated with the following code (forgive the lengthy information).

  26  26  26   26    3    1     0.00000018
  27  27  27   27    3    1     0.00000496
  28  28  28   28    3    1     0.00008922
  29  29  29   29    3    1     0.00102816
  30  30  30   30    3    1     0.00759755
  31  31  31   31    3    1     0.03599289
  32  32  32   32    3    1     0.10934334
  33  33  33   33    3    1     0.21295892
  34  34  34   34    3    1     0.26596953
  35  35  35   35    3    1     0.21295892
  36  36  36   36    3    1     0.10934334
  37  37  37   37    3    1     0.03599289
  38  38  38   38    3    1     0.00759755
  39  39  39   39    3    1     0.00102816
  40  40  40   40    3    1     0.00008922
  41  41  41   41    3    1     0.00000496
  42  42  42   42    3    1     0.00000018

      !n is number of points on the 3D grid.
      n=30   
      rank=3
      cdim=n+n-1
      r1shape(1)=n;r1shape(2)=n;r1shape(3)=n;r2shape=r1shape;
      convshape(1)=cdim;convshape(2)=cdim;convshape(3)=cdim
      r1stride=1; r2stride=1; convstride=1
      allocate(conv_mkl(cdim,cdim,cdim));conv_mkl=0.d0
      ok=.true.
      precis=VSL_CONV_PRECISION_SINGLE
      stat=vslconvsetmode(task,mode)
      stat=vslconvsetinternalprecision(task,precis)
      stat=vsldconvnewtask(task,mode,rank,r1shape,r2shape,convshape)
      stat=vsldconvexec(task,r1,r1stride,r2,r2stride,conv_mkl,convstride)
      ier=vslconvdeletetask(task)
      !normalise this
      conv_mkl=conv_mkl/sum(conv_mkl)
      if(stat/=VSL_STATUS_OK) then
       print*,'Error: Bad status : ',stat
       ok=.false.
      end if
      if(ier/=0) then
       print*,'Failed to destroy task'
      end if
      write(*,*)'Successful? ',ok
      do i=1,cdim
       do j=1,cdim
        do k=1,cdim
         write(23,'(3I4,3I5,F15.8)')i,i-j+3,i-j-k+4,i,j,k,conv_mkl(i,j,k)
        end do
       end do
      end do

I find that the value of j and k are always 3 and 1 respectively, even when I change parameters in function f.  What is so special about 3 and 1?

Thanks.

[Zhang_Z_Intel]

Please take a look at the "Data Allocation" section in the MKL reference manual that describes in great details about multi-dimensional input and output data vectors in convolution routines: http://software.intel.com/en-us/node/470762