gpu arrayfun matrix multipltication, but without using concatenation or varargin varargout

Hello,

I was wondering how i can apply the same matrix multiplication, as encountered in linear unmixing of data sets, upon the elements of multiple arrays, using gpu enabled arrayfun? the problem i am encountering is that only elementary operations are permitted, so i cannot concatenate the scalars that are passed to the function inside arrayfun to construct the vector to be multiplied by a matrix. Furthermore, i cannot use variable inputs/outputs to handle unknown number of datasets.

a working piece of code is as follows with 2 inputs/datasets {A,B} hard coded:

 A = rand(2048,2048,11);
 B = rand(2048,2048,11);
 K = rand(2,2);
 invKmatrix = inv(K);

 function [Aout,Bout] = gpuApplyInvKmatrix(invKmatrix,A,B)
 %GPUAPPLYINVKMATRIX 
   function [output1,output2] = DoStuff(A,B)
       output1 = invKmatrix(1)*A + invKmatrix(3)*B; 
       output2 = invKmatrix(2)*A + invKmatrix(4)*B;
    end
 [Aout,Bout] = arrayfun(@DoStuff,A,B);
 end

the above fulfills the gpu requirement of elementary operations only, but i don't know how to generalize to multiple datasets. for example if the user wants to input 5 datasets...so K matrix of rand(5,5); with dataset inputs of A,B,C,D,E... do i just hard code each case?

stackoverflow question that is similar suggested using a multiple right hand side call, but there seems to be a caveat that this approach allocates too much memory?

Answer 1

[Frederick Chang]

For reference, i wrote an in-place memory version of the above code, going the route of the suggested stackoverflow approach - but for general n datasets. seems like the right solution!

 function cellDataIn = gpuApplyInvKmatrix(kMatrix,cellDataIn)
 %GPUAPPLYINVKMATRIX will take kmatrix and solve the system given varargin
 %datasets.  should be gpu friendly.  datasets can be ndimensional

 sizeADataset    = size(cellDataIn{1});

 cellDataIn = cellfun(@(x) x(:),cellDataIn,'UniformOutput',false);
 cellDataIn = [cellDataIn{:}]';
 cellDataIn = kMatrix \ cellDataIn;
 cellDataIn = num2cell(cellDataIn,2);
 cellDataIn = cellfunNonUniformOutput(@(x) reshape(x,sizeADataset),cellDataIn);