function [ y ] = Lsolve( L,b )
%Lsolve: solve Ly = b (for y) where L is lower triangular and b is a vector
% note: it is not checked that L is actually lower triangular
    [bRows,bCols] = size(b);
    [LRows,LCols] = size(L);
 
    if LRows ~= LCols || bRows ~= LRows || bCols ~= 1
        error('Either L or b is the wrong size.')
    end
    if any( diag(L) == 0 )
        error('There are zeros on the diagonal of L')
    end

    y = b;
    for k = 1:LRows
      for j = 1:k-1
        y(k) = y(k) - L(k,j)*y(j);
      end
      y(k) = y(k)/L(k,k);
    end
end

