vector solve( matrix L, U, vector b)
/*
 * solve LUx=b for x, L is m-by-m, U is m-by-n
 */
{
    vector x, y = b;
    int i, j, q, m = rows(U), n = cols(U);

    q = min( [m, n]);

    for( i=2; i <= q; ++i)			/* forward elimination */
	y(i) = y(i) - L(i,1:i-1) * y(1:i-1);

    x = zeros(n,1);
    x(1:q) = y(1:q);

    if( U(q,q) == 0)
	x(q) = 0;
    else
	x(q) = x(q) / U(q,q);
						/* back substitution */
    for( i = q-1; i >= 1; --i)
	if( U(i,i) == 0)
	    x(i) = 0;
	else
	    x(i) = (x(i) - U(i,i+1:q) * x(i+1:q)) / U(i,i);

    return x;
}
