--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/dlyap.m	Mon Aug 30 15:11:21 1993 +0000
@@ -0,0 +1,83 @@
+function x = dlyap (a, b)
+
+# Usage: x = dlyap (a, b)
+#
+# Solve a x a' - x + b = 0 (discrete Lyapunov equation) for square
+# matrices a and b.  If b is not square, then the function attempts 
+# to solve either
+#
+#  a x a' - x + b b' = 0
+#
+# or
+#
+#  a' x a - x + b' b = 0
+#
+# whichever is appropriate.  Uses Schur decomposition as in Kitagawa
+# (1977).
+
+# Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
+
+  if ((n = is_square (a)) == 0)
+    fprintf (stderr, "warning: dlyap: a must be square");
+  endif
+
+  if ((m = is_square (b)) == 0)
+    [n1, m] = size (b);
+    if (n1 == n)
+      b = b*b';
+      m = n1;
+    else
+      b = b'*b;
+      a = a';
+    endif
+  endif
+
+  if (n != m)
+    fprintf (stderr, "warning: dlyap: a,b not conformably dimensioned");
+  endif
+
+  # Solve the equation column by column.
+
+  [u, s] = schur (a);
+  b = u'*b*u;
+
+  j = n;
+  while (j > 0)
+    j1 = j;
+
+# Check for Schur block.
+
+    if (j == 1)
+      blksiz = 1;
+    elseif (s (j, j-1) != 0)
+      blksiz = 2;
+      j = j - 1;
+    else
+      blksiz = 1;
+    endif
+
+    Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n);
+
+    rhs = reshape (b (:, j:j1), blksiz*n, 1);
+
+    if (j1 < n)
+      rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)');
+      rhs = rhs + reshape (rhs2, blksiz*n, 1);
+    endif
+
+    v = - Ajj\rhs;
+    x (:, j) = v (1:n);
+
+    if(blksiz == 2)
+      x (:, j1) = v ((n+1):blksiz*n);
+    endif
+
+    j = j - 1;
+
+  endwhile
+
+# Back-transform to original coordinates.
+
+  x = u*x*u';
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/kron.m	Mon Aug 30 15:11:21 1993 +0000
@@ -0,0 +1,38 @@
+function x = kron (a, b)
+
+# Usage: x = kron (a, b)
+#
+# Form the Kronecker product of two matrices, defined block by block
+# as 
+#
+#   x = [a(i,j) b]
+
+# Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
+
+  if (nargin == 2)
+
+    [m, n] = size (b);
+    [ma, na] = size (a);
+
+# Do 1st column.
+
+    x = a (1, 1) * b;
+    for ii = 2:ma
+      x = [x; a (ii, 1)*b];
+    endfor
+
+# Do remaining columns.
+
+    for jj = 2:na
+      tmp = a (1, jj)*b;
+      for ii = 2:ma
+	tmp = [tmp; a (ii, jj)*b];
+      endfor
+      x = [x, tmp];
+    endfor
+
+  else
+    error ("usage: kron (a, b)");
+  endif
+
+endfunction