int pivot( matrix A, vector s, int j)
/*
 * find largest element for pivot.  s is the row scale factors.
 */
{
    int i, k = j, m = rows(A);
    double v, u = abs(A(j,j))/s(j);

    for( i = j+1; i <= m; ++i)
	if( (v = abs(A(i,j))/s(i)) > u) {
	    k = i;
	    u = v;
	}
    return k;
}

void lup( matrix A, L, U, vector p)
/*
 * LU reduction using Gaussian Elimination with scaled partial pivoting
 * L*U = A(p)
 */
{
    int i, j, k, q, m = rows(A), n = cols(A);
    double t;
    vector s;

    if( narg != 4) {
	print("lup( A, L, U, p) takes four arguments");
	return;
    }

    U = A;
    L = eye( m, m);
    p = [1:m];
    s = max(abs(A));
    s = s + (s==0);

    q = min( [n, m-1]);

    for( j = 1; j <= q; ++j) {	/* go across the columns */

	if( (k = pivot(U,s,j)) != j) {	/* swap indices and rows */
	    p([j k]) = p([k j]);
	    U([j k]) = U([k j]);
	    if( j > 1)
		L([j k],1:j-1) = L([k j],1:j-1);
	}

	if( U(j,j) == 0)		/* already zeros in col j */
	    continue;

	for( i = j+1; i <= m; ++i) {	/* put zeros below diagonal in col j */

	    t = U(i,j)/U(j,j);
	    U(i,j:n) = U(i,j:n) - t * U(j,j:n);
	    L(i,j) = t;
	}

	U(j+1:m,j) = 0;
    }
}
