% test SpMV using CSR, CSC and symmetric CSR
function err = SpMV_test(n,m,d)
%
% On input: 
%    create random n x m sparse matrix A with density d,
%       to test SpMV with A in CSR and CSC formats
%    create random symmetric n x n sparse matrix As with density d,
%       to test SpMV with As in symmetric CSR format
% On output: err = [
%    errCSR = error of CSR vs "dense" SpMV 
%    errCSC = error of CSC vs "dense" SpMV
%    errSCRC = error of symmetric CSR vs "dense" SpMV ]
%      in all cases, if err = O(1) then results agree to within roundoff
%   Also make spy plots of A, x and y=A*x of and As, x, and ys = As*x
%
A = full(sprand(n,m,d));
figure(1), clf
x = randn(m,1);
ydense = A*x; % get A*x using treating A as dense matrix
subplot(221), spy(A), title(['A: ',int2str(n),' by ',int2str(m)])
subplot(222), spy(x), title('x')
subplot(223), spy(ydense), title('ydense')
%
% create CSR version of A
nnzA = nnz(A);  % number of nonzeros in A
valACSR = zeros(nnzA,1);   % allocate space for nonzeros of A
colindexACSR = zeros(nnzA,1);  % allocate space for column indices of A
rowbeginACSR = zeros(n+1,1);  % allocate space for pointers to rows of A
nnzcnt = 0;
rowbeginACSR(1) = 1;
for i=1:n
    for j=1:m
        if (A(i,j) ~= 0),
            nnzcnt = nnzcnt + 1;
            valACSR(nnzcnt) = A(i,j);
            colindexACSR(nnzcnt) = j;
        end
    end
    rowbeginACSR(i+1) = nnzcnt+1;
end
%
% perform yCSR = A*x using A in CSR format
yCSR = SpMV_CSR(valACSR,colindexACSR,rowbeginACSR,x);
%
% compare CSR answer to dense code; 
% add realmin*eps to denominator to avoid 0/0
errCSR = norm( abs(yCSR-ydense)./(abs(A)*abs(x)+realmin*eps)/eps, 1);
%
% create CSC version of A
valACSC = zeros(nnzA,1);   % allocate space for nonzeros of A
rowindexACSC = zeros(nnzA,1);  % allocate space for row indices of A
colbeginACSC = zeros(m+1,1);  % allocate space for pointers to columns of A
nnzcnt = 0;
colbeginACSC(1) = 1;
for j=1:m
    for i=1:n
        if (A(i,j) ~= 0),
            nnzcnt = nnzcnt + 1;
            valACSC(nnzcnt) = A(i,j);
            rowindexACSC(nnzcnt) = i;
        end
    end
    colbeginACSC(j+1) = nnzcnt+1;
end
%
% perform yCSC = A*x using A in CSC format
yCSC = SpMV_CSC(valACSC,rowindexACSC,colbeginACSC,n,x);
%
% compare CSC answer to dense code; 
% add realmin*eps to denominator to avoid 0/0
errCSC = norm( abs(yCSC-ydense)./(abs(A)*abs(x)+realmin*eps)/eps, 1);
%
% Create random symmetric n x x sparse matrix in dense format
As = full(sprand(n,n,d/2)); As = As + As';
figure(2), clf
x = randn(n,1);
ysdense = As*x; % get As*x using treating As as dense matrix
subplot(221), spy(As), title(['As: ',int2str(n),' by ',int2str(n)])
subplot(222), spy(x), title('x')
subplot(223), spy(ysdense), title('ysdense')
%
% Create CSR version of As, storing 
%    strict upper half, and
%    half value along diagonal, ie As(i,i)/2 if As(i,i) nonzero
% Calling this upper triangular matrix Asu (u for "upper"), we see
% that As = Asu + Asu', so we can just use the above two algorithms
% for multiplying by A in CSR format (for Asu) in CSC format (for Asu')
% without handling the diagonal differently. 
% Alternatively, we could store As = U' + D + U where U is strictly
% upper triangular in CSR format, and D is diagonal, to save
% a number of flops equal to the number of nonzero diagonal entries
%
nnzAs = nnz(triu(A,0));  % number of nonzeros in upper half of As
valAsCSR = zeros(nnzAs,1);   % allocate space for (half) nonzeros of As
colindexAsCSR = zeros(nnzAs,1);  % allocate space for column indices of As
rowbeginAsCSR = zeros(n+1,1);  % allocate space for pointers to rows of As
nnzcnt = 0;
rowbeginAsCSR(1) = 1;
for i=1:n
    for j=i:n
        if (As(i,j) ~= 0),
            nnzcnt = nnzcnt + 1;
            valAsCSR(nnzcnt) = As(i,j);
            if (i==j) valAsCSR(nnzcnt) = As(i,i)/2; end
            colindexAsCSR(nnzcnt) = j;
        end
    end
    rowbeginAsCSR(i+1) = nnzcnt+1;
end
%
% Perform ysCSR = As*x with As in symmetric CSR format
ysCSR = SpMV_CSR(valAsCSR,colindexAsCSR,rowbeginAsCSR,x) + ...
        SpMV_CSC(valAsCSR,colindexAsCSR,rowbeginAsCSR,n,x);
%
% compare symmetric CSR answer to dense code; 
% add realmin*eps to denominator to avoid 0/0
errSCSR = norm( abs(ysCSR-ysdense)./(abs(As)*abs(x)+realmin*eps)/eps, 1);
err = [errCSR, errCSC, errSCSR];