double norm1( matrix a)		/* 1-norm is maximum absolute column sum */
{
    return max(+/abs(a'));
}

double norminf( matrix a)	/* infinite-norm is maximum absolute row sum */
{
    return max(+/abs(a));
}

matrix pinv( matrix a)		/* pseudo-inverse is  V * pinv(S) * U' */
{
    matrix u, s, v;
    int rank;

    rank = svd( a, u, s, v);

    if( rank == 0)
	return zeros( cols(a), rows(a));
    else
	return v(:,1:rank) * diag( 1/diag(s(1:rank,1:rank))) * u(:,1:rank)';
}

matrix N( matrix a)		/* Null space */
{
    matrix v;
    int rank;

    rank = svd( a, :, :, v);		/* don't need u or s */

    if( rank == cols(a))		/* there is no null space! */
	return [];
    else
	return v(:,rank+1:cols(v));
}
