% compare the time required to compute matrix vector products while 
% accessing the contents of the matrix and vector five different ways
% Issues:
%       1. Row based access of the matrix
%       2. Column based access of the matrix
%       3. Row oriented vector operations
%       4. Column oriented vector operations
%       5. Direct matrix vector product through MATLAB's call to BLAS

mvtimes=zeros(4,5);
% 1 access A by rows
for N=1:4
    M=10^N; % matrix sizes of 10 x 10, 100 x 100, 1000 x 1000, 10^4 x 10^4
    s(N)=M; % store size for plotting
    A=rand(M); % create a random matrix 
    b=rand(M,1); % create a random vector
    tic % set clock
    for i=1:M
        c(i)=0;
        for j=1:M
            c(i)=c(i)+A(i,j)*b(j);
        end
    end
    mvtimes(N,1)=toc;
end
% 1 access A by columns
for N=1:4
    M=10^N;
    A=rand(M);
    b=rand(M,1);
    c=zeros(M,1);
    tic
    for j=1:M
        for i=1:M
            c(i)=c(i)+A(i,j)*b(j);
        end
    end
    mvtimes(N,2)=toc;
end
% 1 access A by columns, but vectorially

for N=1:4
    M=10^N;
    A=rand(M);
    b=rand(M,1);
    c=zeros(M,1);
    tic
    for j=1:M
        c(:)=c(:)+A(:,j)*b(j);
    end
    mvtimes(N,3)=toc;
end
% 1 access A by rows, but vectorially

for N=1:4
    M=10^N;
    A=rand(M);
    b=rand(M,1);
    c=zeros(M,1);
    tic
    for i=1:M
        c(i)=c(i)+dot(A(i,:),b(:)); % note dot product row times column
    end
    mvtimes(N,4)=toc;
end
% do it by a single statement (using Matlab internal call to BLAS)
for N=1:4
    M=10^N;
    A=rand(M);
    b=rand(M,1);
    c=zeros(M,1);
    tic
    c=A*b;
    mvtimes(N,5)=toc;
end

mvtimes