% this is just an example, of how one would use the two functions
% at the end (add_to and action), to hande a sparse matrix A.
    % We'll only store the non-zero elements. To that end we have
    % for k in [1,c] : A(i,k,1) = colum index j, A(i,k,2) = value of Aij

% n - dimension of the matrix, c - max number of non-zero elements per row.
function error_norm = sparse_example(n,c)
    % create an empty matrix. 
    A = zeros(n,c,2);
    % we want to test our code, so we'll keep the matrix in full form too
    Afull = zeros(n,n);
    
    % now let us randomly fill our matrices:
    for i=1:n
        for k=1:c
            j = ceil(rand(1)*n);
            value = randn(1);
            Afull(i,j) = Afull(i,j) + value;
            A = add_to(A, i,j, value);
        end
    end
    
    
    % supposedly now both A and Afull represent the same matrix
    % this is what they look:
    Afull
    A
    
    % lets test them against some random vector.
    x = rand(n,1);
    
    b1 = action(A,x);
    b2 = Afull*x;
    
    b1
    b2
    b1-b2
    % the vectors b1 and b2 differ by a vector with norm:
    error_norm = norm(b1-b2);
end

%The following two functions are ALL that we need to handle
% a sparse matrix A, without storing all the zeros in our 2D 
% finite element code... (think about how to (easily!) impose 
% Dirichlet boundary conditions)

%The main function! A(i,j) = A(i,j) + d;
function newA = add_to(A,i,j,d)
    col = size(A,2);
    for k = 1:col
        % entry for the j-th column found add 'd' to it.
        if (A(i,k,1)==j)
            A(i,k,2) = A(i,k,2)+d;
            newA = A; return;
        end
        % empty entry found - create one for the j-th column.
        if (A(i,k,1)==0)
            A(i,k,1) = j; A(i,k,2) = A(i,k,2)+d;
            newA = A; return;
        end
    end
    
    % we'll only get this far, if all entries were filled, and neither
    %  was for the j-th column. means the matrix A is not as sparse as we
    %  thought !
end

% the action function b=Ax
function b = action(A,x)
    m = size(A,1);
    b = zeros(m,1);
    for i=1:m
        for j=1:size(A,2)
            if (A(i,j,1)~=0)
                b(i) = b(i) + A(i,j,2)*x(A(i,j,1));
            end
        end
    end
end