title: Is there a way to optimize/vectorize these loops over elements in a 3D array, without requiring significantly more memory?
body: <p>I am looking for some help to speed up some Matlab code.</p>
  <p>A minimal example is shown below. The code is doing some calculations on a 3D matrix, defined on [x,y,z] coordinates. When using the profiler, the inner loop over <code>ind</code> is the time-consuming part, and so I am wondering if this loop can either be optimized, or removed/vectorized completely.</p>
  <pre class="lang-matlab prettyprint-override"><code>Nx = 8; % Number of grid points
  Ny = 6;
  Nz = 4;
  Ntot = Nx*Ny*Nz;
  
  xvals = rand(1,Nx); % Create grid vectors
  yvals = rand(1,Ny);
  zvals = rand(1,Nz);
  
  input_vec = rand(Ny,Nx,Nz); % Generate a dummy 3D matrix ( meshgrid convention, [y,x,z] )
  input_vec = reshape( permute(input_vec,[3,1,2]) , [Ntot 1]); % Unwrap to 1D, so z cycles fastest, then y, then x
  
  C1 = 5; % Loop counters
  C2 = 6;
  C3 = 7;
  
  output_vec = zeros(Ntot,1); % Preallocate
  temp_vec = zeros(Ntot,1);
  
  for cnt1 = 1:C1
      for cnt2 = 1:C2
          for cnt3 = 1:C3
              
              factor1 = xvals*cnt1; % Calculate some vectors which depend on cnt variables
              factor2 = yvals*cnt2;
              factor3 = zvals*cnt3;
              
              for ind = 1:Ntot % Loop over every grid point
                  j1 = floor( floor((ind-1) / Nz) / Ny) + 1; % +1 and -1 's account for Matlab's [1] indexing
                  j2 = mod( floor((ind-1)/Nz) , Ny ) + 1;
                  j3 = mod( (ind-1), Nz ) + 1;
                  temp_vec(ind) = input_vec(ind) * factor1(j1)*factor2(j2)*factor3(j3);
              end
              output_vec = output_vec + temp_vec;
          end
      end
  end
  </code></pre>
  <p>In my real application, the number of points is more like 1024x1024x512, and so I have tried to avoid using lots of <code>meshgrid</code> formatted variables (which contain a lot of repeated information) in order to keep the memory requirements down - this is the reason that the 3D array has been unwrapped to 1D in the code above. For example, one solution might be to precalculate all the <code>j1,j2,j3</code> values like so</p>
  <pre class="lang-matlab prettyprint-override"><code>j1 = 1:Nx;
  j2 = 1:Ny;
  j3 = 1:Nz;
  [J1,J2,J3] = meshgrid(j1,j2,j3);
  J1 = reshape( permute(J1,[3,1,2]) , [Ntot 1]); % Unwrap to 1D, so z cycles fastest, then y, then x 
  J2 = reshape( permute(J2,[3,1,2]) , [Ntot 1]); 
  J3 = reshape( permute(J3,[3,1,2]) , [Ntot 1]); 
  </code></pre>
  <p>but this requires much more RAM than calculating a single value of j each time depending on the value of ind.</p>
  <p>Can anyone help with a better/faster (but still memory efficient) way to do this? Thank you.</p>