gpuArray related memory issue

Hello, i am trying to do some calculations on the GPU but unfortunately possibly due to size of matrix (1010x1010x1601) I get GPU memory run out error (4gb Nvidia), and without GPU this calculation takes days. Is there any way to modify gpuArray for separate it subArrays to overcome this memory problem?

Thank you for your assistance.

function proj2d = projection(data3d,param, iview)  
angle_rad = param.deg(iview)/360*2*pi;
proj2d = (zeros(param.nu,param.nv,'single'));

[uu,vv] = meshgrid(param.us,param.vs);
[xx,yy] = meshgrid(param.xs,param.ys);

if param.gpu == 1
    data3d = gpuArray(single(data3d));
    rx = gpuArray(((xx.*cos(angle_rad) - yy.*sin(angle_rad)) - xx(1,1))/param.dx + 1);
    ry = gpuArray(((xx.*sin(angle_rad) + yy.*cos(angle_rad)) - yy(1,1))/param.dy + 1);
else
    rx = (((xx.*cos(angle_rad) - yy.*sin(angle_rad)) - xx(1,1))/param.dx + 1);
    ry = (((xx.*sin(angle_rad) + yy.*cos(angle_rad)) - yy(1,1))/param.dy + 1);
end

for iz = 1:param.nz   

      data3d(:,:,iz) = interp2(data3d(:,:,iz),rx,ry, param.interptype);

end

data3d(isnan(data3d))=0;

data3d = permute(data3d,[1 3 2]);

[xx,zz] = meshgrid(param.xs,param.zs);

for iy = 1:param.ny

      Ratio = (param.ys(iy)+param.DSO)/(param.DSD);

      pu = uu*Ratio;
      pv = vv*Ratio;    

      pu = (pu - xx(1,1))/(param.dx)+1; 
      pv = (pv - zz(1,1))/(param.dz)+1; 

      if param.gpu == 1
          tmp = gather(interp2(gpuArray(single(data3d(:,:,iy))),gpuArray(single(pv)),gpuArray(single(pu)),param.interptype));
      else
          tmp = (interp2((single(data3d(:,:,iy))),(single(pv)),(single(pu)),param.interptype));
      end

      tmp(isnan(tmp))=0;

      proj2d = proj2d + tmp';
  end

dist = sqrt((param.DSD)^2 + uu.^2 + vv.^2)./(param.DSD)*param.dy;

proj2d = proj2d .* dist';

[Joss Knight]

It looks to me like a single call to interpn is what you need, not multiple calls to interp2 with a permute in the middle. And your second loop is summing over the y-axis, so why not use sum?

Answer 1

[Joss Knight]

Your first chunk of operations are element-wise, so you can divide the arrays up however you like (perhaps along the z-axis?) and process it in chunks, gathering each result to free up GPU memory.

Your two loops are over the 3rd dimension of your data, so you could move each slice of data to the GPU, process it, and gather it back.

[Emre Topal]

Hello Mr.Knight, thank you for your fast response. I actually tried to divide arrays along z axis but then the further parts of algorithm started to give mismatch on matrix dimension error, in theory your suggestion should work but i really don't know why it didn't work for my case. Do you have any further suggestion on how can successfully divide arrays and process it in chunks? I know, I am asking too much but really need some help because just one reconstruction takes almost two days.

[Joss Knight]

Well, I mean, just don't put data3d on the GPU and then go

   for iz = 1:param.nz   
         slice = gpuArray(data3d(:,:,iz));
         data3d(:,:,iz) = gather(interp2(slice,rx,ry, param.interptype));
   end

But really you should be using interpn.

[Emre Topal]

Many thanks for your comment and assistance, this is a simple solution, and I think, I was trying more complicated and useless approaches to divide arrays than your suggestion. I will also try to use interpn within today, but I think I need to use interpn2 because the input argument is a double array.

[Joss Knight]

You're permuting your data just to move the z-axis to dim 2 so that you can call meshgrid and interp2. With sensible use of ndgrid and interpn, or meshgrid and interp3, your code could be essentially identical but without the permute which is very slow.