void lu( A, L, U)	/* Doolittle LU reduction by rows */
{
    int i, j, k, n = rows(A);

    if( narg != 3) {
	print("lu( A, L, U) takes 3 arguments");
	return;
    }

    if( n != cols(A) ) {
	print("lu( A, L, U) takes square matrix A");
	return;
    }

    U = A;
    L = eye( n, n);

    for( i = 2; i <= n; ++i) {
	L(i,1) = U(i,1)/U(1,1);
	for( j = 2; j <= i-1; ++j) {
	    L(i,j) = (U(i,j) - L(i,1:j-1)*U(1:j-1,j)) / U(j,j);
	}
	for( k = 1; k <= i-1; ++k)
	    U(i) = U(i) - L(i,k) * U(k);
    }
}
