0

I have a 3 dimensional array which I am seeking to shift, multiply (element wise) with the un-shifted array, and sum all the products. For a given set of the three integers defining the shift (i,j,k) the sum corresponding to this triple will be stored in another 3 dimensional array at said location (i,j,k).

I am attempting to do this for all possible values of i,j,k. Is there a way to do this without explicitly using for/while loops? I am highly interested in minimizing the computational time but IDL doesnt understand the current instructions I am providing:

FUNCTION ccs, vecx, vecy

TIC

l = 512
km = l/2
corvec = fltarr(km,km,km)
N = float(l)^3


corvec[*,*,*] = (total(vecx*shift(vecx,*,*,*))/N)+  (total(vecy*shift(vecy,*,*,*))/N)

return, corvec

TOC
end

IDL appears to only accept scalars or 1 element arrays as possible shifting arguments.

asynchronos
  • 583
  • 2
  • 14
Anthony Zealous
  • 25
  • 6

1 Answers1

0

Is the below what you are trying to do? It's not fast, but it will give something to test faster solutions against.

function ccs, vecx, vecy
  compile_opt strictarr

  tic

  l = 512
  km = l/2
  corvec = fltarr(km, km, km)
  n = float(l)^3


  for i = 0L, km - 1L do begin
    for j = 0L, km - 1L do begin
      for k = 0L, km - 1L do begin
        corvec[i, j, k] = total(vecx * shift(vecx, i, j, k)) / n +  total(vecy * shift(vecy, i, j, k)) / n
      endfor
    endfor
  endfor

  toc

  return, corvec
end

If this is what you are trying to do, I think it will be tough in IDL. Generally, vectorized solutions without loops for these types of problems need to create intermediate arrays — in this case, of all possible shifts of vecx and vecy, which would be a very large intermediate array. My suggestion would be to implement in C and call from IDL, but maybe I am missing something.

mgalloy
  • 2,356
  • 1
  • 12
  • 10
  • That's what I didn't want to do but have to unfortunately. I was told by my research advisor to figure out the fastest way to do it in IDL in order to help learn as Im new to it. I could easily just bang it out brute force in C or Fortran but was told not to. – Anthony Zealous Nov 06 '17 at 00:32