/*  -*- MaTX -*-
 *
 *  NAME
 *      dgramian() - Discrete-time controllability and observability gramians
 *
 *  SYNOPSIS
 *      G = dgramian(A, B)
 *         Matrix G;
 *         Matrix A, B;
 *
 *  DESCRIPTION
 *       dgramian(A,B) returns the controllability gramian.
 *       dgramian(A#,C#) returns the observability gramian.
 *
 *  SEE ALSO
 *      gramian
 */

Func Matrix dgramian(A,B)
	Matrix A,B;
{
	Matrix G,D,U,V,P;
	String msg;

	if (length((msg = abcdchk(A, B))) > 0) {
		error("dgramian(): " + msg);
	}

	{U,D,V} = svd(B);
	P = dlyap(U#*A*U, D*D#);
	G = U*P*U#;
	G = (G + G#)/2;
	return G;
}