Code covered by the BSD License

# frontal

### Felipe G. Nievinski (view profile)

15 Mar 2011 (Updated )

Do: C = frontal_mtimes(A, b); not: for k=1:size(A,3), C(:,:,k) = A(:,:,k) * b(:,:,k); end

frontal_mtimes (A, B)
```function C = frontal_mtimes (A, B)
if isscalar(A) || isscalar(B)
C = frontal_times(A, B);
return;
end
disable_lp = true;  % case ((mC*nC) < p).

[mA,nA,pA] = size(A);
[mB,nB,pB] = size(B);

if (nA ~= mB)
error('frontal:mtimes:innerdim', ...
'Inner matrix dimensions must agree.')
end
if ~( (pA == pB) || (pA == 1) || (pB == 1) )
error ('frontal:mtimes:thirddim', ...
'Third matrix dimensions must agree.')
end

mC = mA;
nC = nB;
p = max(pA, pB);

if (p == 1)
C = A * B;
elseif (pA ~= pB) && (pB == 1)
C = frontal_transpose(...
frontal_mtimes(...
frontal_transpose(B), ...
frontal_transpose(A)));
elseif (pA ~= pB) && (pA == 1) && (nB == 1)
C = permute(A * permute(B, [1,3,2]), [1,3,2]);
elseif (pA ~= pB) && (pA == 1) && ((mC*nC) < p)
C = zeros(mC,nC,p);
for i=1:mC, for j=1:nC
% following the definition of matrix product:
C(i,j,:) = sum(repmat(reshape(A(i,:),nA,1),[1,1,p]) .* B(:,j,:),1);
end, end
elseif (pA ~= pB) && (pA == 1) && ((mC*nC) >= p)
C = zeros(mC,nC,p);
for k=1:p
C(:,:,k) = A * B(:,:,k);
end
elseif ~disable_lp && ((mC*nC) < p)
C = zeros(mC,nC,p);
for i=1:mC, for j=1:nC
% following the definition of matrix product:
C(i,j,:) = sum(reshape(A(i,:,:),nA,1,p) .* B(:,j,:), 1);
%% for more detailed profiling:
%temp1 = A(i,:,:);
%temp2 = reshape(temp1,nA,1,p);
%temp3 = B(:,j,:);
%temp4 = temp2 .* temp3;
%C(i,j,:) = sum(temp4, 1);
end, end
%elseif ((mC*nC) >= p)
else
try
%try  N = maxNumCompThreads;  catch N = 1;  end
N = 1;
if (N == 1)
C = frontal_mtimes_helper(A, B);
return;
end
% if only cellfun were multi-threaded...
ppp = doit (p, N);
Ac = mat2cell(A, mA, nA, ppp);
Bc = mat2cell(B, mB, nB, ppp);
Cc = cellfun(@frontal_mtimes_helper, Ac, Bc, ...
'UniformOutput',false);
C = cell2mat(Cc);
catch s
%s = lasterror;
if ~any(strcmp(s.identifier, {...
'MATLAB:UndefinedFunction', ...  % _helper not compiled.
'frontal:mtimes:diffClasses', ...
'frontal:mtimes:nonFloat', ...
'frontal:mtimes:complex'}))
rethrow(s);
end
%whos A B  % DEBUG
C = zeros(mC,nC,p);
for k=1:p
C(:,:,k) = A(:,:,k) * B(:,:,k);
end
end
end
end

%!test
%! % frontal_mtimes
%! warning('off', 'test:noFuncCall');

%!test
%! lasterror('reset')

%!error
%! frontal_mtimes(zeros(1,1,3), zeros(1,1,4));

%!test
%! s = lasterror;
%! myassert(s.identifier, 'frontal:mtimes:thirddim')

%!test
%! lasterror('reset')

%!error
%! frontal_mtimes(zeros(3,2,3), zeros(3,2,3));

%!test
%! s = lasterror;
%! myassert(s.identifier, 'frontal:mtimes:innerdim')

%!test
%! warning('on', 'test:noFuncCall');

%!shared
%!  p = ceil(10*rand);
%! mA = ceil(10*rand);
%! nA = ceil(10*rand);
%! mB = nA;  % inner dimension must agree.
%! nB = ceil(10*rand);
%! A = rand(mA, nA, p);
%! B = rand(mB, nB, p);
%! A1 = A(:,:,1);
%! B1 = B(:,:,1);
%! mC = mA;
%! nC = nB;

%!test
%! % frontal matrices with only one page
%! % can be treated as 2d matrices.
%! C = frontal_mtimes(A1, B1);
%! C2 = A1 * B1;
%! myassert (C, C2);

%!test
%! % repeated frontal pages yield repeated frontal pages:
%! C = frontal_mtimes(repmat(A1, [1,1,p]), repmat(B1, [1,1,p]));
%! C2 = repmat(A1*B1, [1,1,p]);
%! %p, C, C2  % DEBUG
%! myassert (C, C2, -10*eps);

%!test
%! % matrix with only one frontal page can be multiplied with
%! % matrix with multiple frontal pages.
%! myassert(frontal_mtimes(A1, B), frontal_mtimes(repmat(A1, [1,1,p]), B), -10*eps);
%! myassert(frontal_mtimes(A, B1), frontal_mtimes(A, repmat(B1, [1,1,p])), -10*eps);

%!test
%! % empty input, empty output.
%! myassert(isempty(frontal_mtimes([], [])))
%! C = frontal_mtimes(zeros(0,0,1), zeros(0,3,3));
%! myassert(isempty(C));
%! myassert(C, repmat(zeros(0,0)*zeros(0,3), [1,1,3]))

%!test
%! % empty input, non-empty output.
%! myassert(isempty(frontal_mtimes([], [])))
%! C = frontal_mtimes(zeros(2,0,3), zeros(0,3,3));
%! myassert(~isempty(C));
%! myassert(C, repmat(zeros(2,0)*zeros(0,3), [1,1,3]))

%!test
%! % frontal_mtimes
%! lasterror('reset')

%!error
%! % empty input still have to have inner dimensions in agreement.
%! frontal_mtimes(zeros(2,0), zeros(2,2));

%!test
%! % frontal_mtimes
%! s = lasterror;
%! myassert(s.identifier, 'frontal:mtimes:innerdim')

%!test
%! % alternative implementation:
%! function C = frontal_mtimes2 (A, B)
%!     for i=1:mC, for j=1:nC
%!         % following the definition of matrix product:
%!         C(i,j,:) = sum(reshape(A(i,:,:),nA,1,p) .* B(:,j,:), 1);
%!         %% for more detailed profiling:
%!         %temp1 = A(i,:,:);
%!         %temp2 = reshape(temp1,nA,1,p);
%!         %temp3 = B(:,j,:);
%!         %temp4 = temp2 .* temp3;
%!         %C(i,j,:) = sum(temp4, 1);
%!     end, end
%!     %% caught and fixed a bug with alternative implementation:
%!     %C = frontal_mtimes(ones(1,1,2), 3*ones(1,2,2));
%!     %myassert(C, repmat([3,3], [1,1,2]))
%!     %C = frontal_mtimes(ones(2,1,2), 3*ones(1,1,2));
%!     %myassert(C, repmat([3;3], [1,1,2]))
%! end
%!
%! tic; C1 = frontal_mtimes (A,B); t1 = toc;
%! tic; C2 = frontal_mtimes2(A,B); t2 = toc;
%!
%! myassert (C1, C2, -10*eps);
%! %[t1, t2, cputime_tol]

%!test
%! % there is an alternative implementation for the for i=1:p loop.
%! function C = frontal_mtimes2 (A, B)
%!     C = cell2mat(cellfun(@mtimes, ...
%!         mat2cell(A, m, n, ones(p,1)), ...
%!         mat2cell(B, m, n, ones(p,1)), ...
%!         'UniformOutput',false));
%! end
%! m = 3;
%! n = 3;
%! p = 10000/2;
%! A = repmat(eye(m,n), [1,1,p]);
%! B = rand(m,n,p);
%! %[mC*nC, p]  % DEBUG
%!
%! tic; C1 = frontal_mtimes (A,B); t1 = toc;
%! tic; C2 = frontal_mtimes2(A,B); t2 = toc;
%!
%! myassert (C1, C2, -10*eps);
%! %[t1, t2, cputime_res]  % DEBUG

%!test
%! % complex-valued input:
%! A1 = complex(A1, A1);
%! B1 = complex(B1, B1);
%! C = frontal_mtimes(A1, B1);
%! C2 = A1 * B1;
%! myassert (C, C2);

%!test
%! % A is a single frontal page & B is a frontal collection of column vectors,
%! % then the multiplication is faster than replicating A:
%! pA = 1;
%! pB = ceil(10*rand);
%! pB = 100000 + pB;
%! mA = ceil(10*rand);
%! nA = ceil(10*rand);
%! mB = nA;  % inner dimension must agree.
%! nB = 1;
%! A = repmat(rand(mA, nA, pA), [1,1,pB]);
%! B = rand(mB, nB, pB);
%! tol = 2*cputime_res() + sqrt(eps);
%!
%! tic; C1 = frontal_mtimes (A,B); t1 = toc;
%! tic; C2 = frontal_mtimes (A(:,:,1),B); t2 = toc;
%!
%! myassert (C1, C2, -10*eps);
%! %[t1, t2, cputime_res]  % DEBUG
%! myassert (t2 <= t1 || abs(t2-t1) <= tol);

```