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

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

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

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