/* -*- MaTX -*-
 *
 *  NAME
 *      dlyap() - Solution of discrete-time Lyapunov equation
 *
 *  SYNOPSIS
 *      P = dlyap(A,Q)
 *         Matrix P;
 *         Matrix Q,Q;
 *
 *  DESCRIPTION
 *       dlyap() returns the solution of the discrete-time 
 *       Lyapunov equation:
 *
 *          A*P*A# + Q = P
 *
 *  SEE ALSO
 *      lyap
 */

Func Matrix dlyap(A, Q)
	Matrix A, Q;
{
	Integer n;
	Matrix A2, Q2, P;
	String msg;

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

	n = Rows(A);
	A2 = (A + I(n)) \ (A - I(n));
	Q2 = (I(n) - A2) * Q * (I(n) - A2#)/2;
	P = lyap(A2, Q2);
	return P;
}