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