vector lsolve( matrix L, vector b)
/*
 * solve Lx=b for x, L is r-by-r non-singular lower-triangular
 */
{
    vector x = b;
    int i, r = rows(L);

    x(1) = x(1) / L(1,1);

    for( i = 2; i <= r;  ++i)
	x(i) = (x(i) - L(i,1:i-1) * x(1:i-1)) / L(i,i);

    return x;
}
