/*  -*- MaTX -*-
 *
 *  NAME
 *      dbalreal() - Discrete-time balanced realization
 *
 *  SYNOPSIS
 *      {Ab,Bb,Cb,G,T} = dbalreal(A,B,C)
 *         Matrix Ab,Bb,Cb,G,T;
 *         Matrix A,B,C;
 *
 *  DESCRIPTION
 *       dbalreal() returns a balanced realization of the given system.
 *
 *       dbalreal() also returns a vector G containing the diagonal
 *       of the gramian, and matrix T for the similarity transformation.
 *
 *  SEE ALSO
 *      balreal and minreal
 */

Func List dbalreal(A, B, C)
	Matrix A, B, C;
{
	String msg;
	Matrix Ab, Bb, Cb, G, T;
	Matrix Gc, Go, R, U, D, V;

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

	Gc = dgramian(A, B);
	Go = dgramian(A#, C#);
	R = chol(Gc);

	{U, D, V} = svd(R * Go * R#);
	D = D .* sgn(U# * V);
	T = R# * V * vec2diag(diag2vec(D) .^ (-0.25));
	Ab = T \ A * T;
	Bb = T \ B;
	Cb = C * T;
	G = diag2vec(dgramian(Ab,Bb))#;

	return {Ab, Bb, Cb, G, T};
}