% Test routine for mtimesx, multi-dimensional speed and equality to MATLAB
%******************************************************************************
% 
%  MATLAB (R) is a trademark of The Mathworks (R) Corporation
% 
%  Function:    mtimesx_test_nd
%  Filename:    mtimesx_test_nd.m
%  Programmer:  James Tursa
%  Version:     1.40
%  Date:        October 4, 2010
%  Copyright:   (c) 2009,2010 by James Tursa, All Rights Reserved
%
%  This code uses the BSD License:
%
%  Redistribution and use in source and binary forms, with or without 
%  modification, are permitted provided that the following conditions are 
%  met:
%
%     * Redistributions of source code must retain the above copyright 
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright 
%       notice, this list of conditions and the following disclaimer in 
%       the documentation and/or other materials provided with the distribution
%      
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
%  POSSIBILITY OF SUCH DAMAGE.
%
%  Syntax:
%
%    A = mtimesx_test_nd    % default n=4 is used
%    A = mtimesx_test_nd(n)
%
%    where n = number of repetitions (should be 4 <= n <= 100)
%
%  Output:
%
%    Prints out speed and equality test results.
%    A = cell array with tabled results.
%
%  2010/Oct/04 --> 1.40, Added OpenMP support for custom code
%                        Expanded sparse * single and sparse * nD support
%
%--------------------------------------------------------------------------

function Cr = mtimesx_test_nd(n)
mtimesx; % load the mex routine into memory
if( nargin == 0 )
    n = 4;
else
    n = floor(n);
    if( ~(n >= 4 && n <= 100) )
        n = 4;
    end
end
cn = sprintf('%g',n);

disp(' ');
disp('MTIMESX multi-dimensional equality and speed tests');
disp('--------------------------------------------------');
disp(' ');
disp('(M x K) * ( K x N) equality tests, SPEED mode, M,K,N <= 4');
trans = 'NGTC';
cmpx = {'real ','cmpx '};
mtimesx('speed');
smallok = true;
for m=1:4
    for k=1:4
        for n=1:4
            for transa=1:4
                if( transa <= 2 )
                    ma = m;
                    ka = k;
                else
                    ma = k;
                    ka = m;
                end
                for transb=1:4
                    if( transb <= 2 )
                        kb = k;
                        nb = n;
                    else
                        kb = n;
                        nb = k;
                    end
                    for cmplxa=1:2
                        if( cmplxa == 1 )
                            A = floor(rand(ma,ka)*100+1);
                        else
                            A = floor(rand(ma,ka)*100+1) + floor(rand(ma,ka)*100+1)*1i;
                        end
                        for cmplxb=1:2
                            if( cmplxb == 1 )
                                B = floor(rand(kb,nb)*100+1);
                            else
                                B = floor(rand(kb,nb)*100+1) + floor(rand(kb,nb)*100+1)*1i;
                            end
                            Cm = mtimesx_sparse(A,trans(transa),B,trans(transb));
                            Cx = mtimesx(A,trans(transa),B,trans(transb));
                            if( isequal(Cm,Cx) )
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  EQUAL']);
                            else
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  NOT EQUAL']);
                                smallok = false;
                            end
                        end
                    end
                end
            end
        end
    end
end

if( mtimesx('openmp') )
disp(' ');
disp('(M x K) * ( K x N) equality tests, SPEEDOMP mode, M,K,N <= 4');
mtimesx('speedomp');
smallokomp = true;
for m=1:4
    for k=1:4
        for n=1:4
            for transa=1:4
                if( transa <= 2 )
                    ma = m;
                    ka = k;
                else
                    ma = k;
                    ka = m;
                end
                for transb=1:4
                    if( transb <= 2 )
                        kb = k;
                        nb = n;
                    else
                        kb = n;
                        nb = k;
                    end
                    for cmplxa=1:2
                        if( cmplxa == 1 )
                            A = floor(rand(ma,ka)*100+1);
                        else
                            A = floor(rand(ma,ka)*100+1) + floor(rand(ma,ka)*100+1)*1i;
                        end
                        A = reshape(repmat(A,1000,1),ma,ka,1000);
                        for cmplxb=1:2
                            if( cmplxb == 1 )
                                B = floor(rand(kb,nb)*100+1);
                            else
                                B = floor(rand(kb,nb)*100+1) + floor(rand(kb,nb)*100+1)*1i;
                            end
                            B = reshape(repmat(B,1000,1),kb,nb,1000);
                            Cm = mtimesx_sparse(A(:,:,1),trans(transa),B(:,:,1),trans(transb));
                            Cx = mtimesx(A,trans(transa),B,trans(transb));
                            if( isequal(Cm,Cx(:,:,1)) )
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  EQUAL']);
                            else
                                disp(['(' cmpx{cmplxa} num2str(m) ' x ' num2str(k) ')' trans(transa) ...
                                        ' * (' cmpx{cmplxb} num2str(k) ' x ' num2str(n) ')' trans(transb) '  NOT EQUAL']);
                                smallokomp = false;
                            end
                        end
                    end
                end
            end
        end
    end
end
end

disp(' ');
if( smallok )
    disp('All small matrix multiplies are OK in SPEED mode');
else
    disp('ERROR --> One or more of the small matrix multiplies was not equal in SPEED mode');
end
if( mtimesx('openmp') )
if( smallokomp )
    disp('All small matrix multiplies are OK in SPEEDOMP mode');
else
    disp('ERROR --> One or more of the small matrix multiplies was not equal in SPEEDOMP mode');
end
end

disp(' ');
disp(['mtimesx multi-dimensional test routine using ' cn ' repetitions']);

if( mtimesx('OPENMP') )
    topm = 6;
else
    topm = 4;
end
Cr = cell(6,topm+1);
Cr{1,1} = 'All operands real';

for m=2:topm+1
if( m == 2 )
    mtimesx('BLAS');
elseif( m == 3 )
    mtimesx('LOOPS');
elseif( m == 4 )
    mtimesx('MATLAB');
elseif( m == 5 )
    mtimesx('SPEED');
elseif( m == 6 )
    mtimesx('LOOPSOMP');
else
    mtimesx('SPEEDOMP');
end
Cr{1,m} = mtimesx;

disp(' ');
disp('--------------------------------------------------------------');
disp('--------------------------------------------------------------');
disp(' ');
disp(['MTIMESX mode: ' mtimesx]);
disp(' ');
disp('(real 3x5x1x4x3x2x1x8) * (real 5x7x3x1x3x2x5) example');
Cr{2,1} = '(3x5xND) *(5x7xND)';
A = rand(3,5,1,4,3,2,1,8);
B = rand(5,7,3,1,3,2,5);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(3,7,3,4,3,2,5,8);
for k1=1:3
    for k2=1:4
        for k3=1:3
            for k4=1:2
                for k5=1:5
                    for k6=1:8
                        Cm(:,:,k1,k2,k3,k4,k5,k6) = A(:,:,1,k2,k3,k4,1,k6) * B(:,:,k1,1,k3,k4,k5);
                    end
                end
            end
        end
    end
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{2,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 3x3x1000000) * (real 3x3x1000000) example');
Cr{3,1} = '(3x3xN) *(3x3xN)';
A = rand(3,3,1000000);
B = rand(3,3,1000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(3,3,1000000);
for k1=1:1000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{3,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 2x2x2000000) * (real 2x2x2000000) example');
Cr{4,1} = '(2x2xN) *(2x2xN)';
A = rand(2,2,2000000);
B = rand(2,2,2000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(2,2,2000000);
for k1=1:2000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{4,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

disp(' ');
disp('--------------------------------------------------------------');
disp('(real 2x2x2000000) * (real 1x1x2000000) example');
Cr{5,1} = '(2x2xN) *(1x1xN)';
A = rand(2,2,2000000);
B = rand(1,1,2000000);
% mtimes
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = zeros(2,2,2000000);
for k1=1:2000000
    Cm(:,:,k1) = A(:,:,k1) * B(:,:,k1);
end
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{5,m} = faster;
disp(' ');
disp(['mtimes  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB mtimes' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches mtimes:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match mtimes:  NOT EQUAL , max diff = ' num2str(dx)])
end

try
bsxfun(@times,1,1);
Cr{6,1} = 'above vs bsxfun';
A = rand(2,2,2000000);
B = rand(1,1,2000000);
% bsxfun
tm = zeros(1,n);
for k=1:n
clear Cm
A(1) = 2*A(1);
B(1) = 2*B(1);
tic
Cm = bsxfun(@times,A,B);
tm(k) = toc;
end
% mtimesx
tx = zeros(1,n);
for k=1:n
clear Cx
tic
Cx = mtimesx(A,B);
tx(k) = toc;
end
% results
tm = median(tm);
tx = median(tx);
if( tx < tm )
    faster = sprintf('%7.1f',100*(tm)/tx-100);
    slower = '';
else
    faster = sprintf('%7.1f',-(100*(tx)/tm-100));
    slower = ' (i.e., slower)';
end
Cr{6,m} = faster;
disp(' ');
disp(['bsxfun  Elapsed time ' num2str(tm) ' seconds.']);
disp(['MTIMESX Elapsed time ' num2str(tx) ' seconds.']);
disp(['MTIMESX ' mtimesx ' mode is ' faster '% faster than MATLAB bsxfun with @times' slower])
if( isequal(Cx,Cm) )
    disp(['MTIMESX ' mtimesx ' mode result matches bsxfun with @times:  EQUAL'])
else
    dx = max(abs(Cx(:)-Cm(:)));
    disp(['MTIMESX ' mtimesx ' mode result does not match bsxfun with @times:  NOT EQUAL , max diff = ' num2str(dx)])
end
catch
    disp('Could not perform comparison with bsxfun, possibly because your version of');
    disp('MATLAB does not have it. You can download a substitute for bsxfun from the');
    disp('FEX here: http://www.mathworks.com/matlabcentral/fileexchange/23005-bsxfun-substitute');
end

end

disp(' ');
disp('Percent Faster Results Table');
disp(' ');
disp(Cr);

disp(' ');
disp('Done');
disp(' ');

end