# HG changeset patch
# User jwe
# Date 745104643 0
# Node ID 98eb51c870b27906f99200641e7edf50a9ddf38c
# Parent  4565ad8b4697431f117bf0211f67973ad8fb9999
[project @ 1993-08-11 21:29:50 by jwe]

diff -r 4565ad8b4697 -r 98eb51c870b2 scripts/general/is_square.m
--- a/scripts/general/is_square.m	Wed Aug 11 21:21:34 1993 +0000
+++ b/scripts/general/is_square.m	Wed Aug 11 21:30:43 1993 +0000
@@ -9,7 +9,7 @@
 
   if (nargin == 1)
     [nr, nc] = size (x);
-    if( nr == nc) 
+    if (nr == nc) 
       retval = nr;
     else
       retval = 0;
diff -r 4565ad8b4697 -r 98eb51c870b2 scripts/linear-algebra/qzhess.m
--- a/scripts/linear-algebra/qzhess.m	Wed Aug 11 21:21:34 1993 +0000
+++ b/scripts/linear-algebra/qzhess.m	Wed Aug 11 21:30:43 1993 +0000
@@ -1,45 +1,56 @@
-function [aa,bb,q,z] = qzhess(a,b)
+function [aa, bb, q, z] = qzhess (a, b)
 
-% compute the qz decomposition of the matrix pencil (a - lambda b)
-% result: (for MATLAB compatibility):
-%	q*a*z = aa, q*b*z = bb
-%	q, z orthogonal, 
-%	v = matrix of generalized eigenvectors
-%
-% This ought to be done in a compiled program
-% Algorithm taken from Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed.
+# Compute the qz decomposition of the matrix pencil (a - lambda b)
+#
+# result: (for Matlab compatibility):
+#
+#   aa = q*a*z and bb = q*b*z, with q, z orthogonal, and
+#   v = matrix of generalized eigenvectors.
+#
+# This ought to be done in a compiled program
+#
+# Algorithm taken from Golub and Van Loan, Matrix Computations, 2nd ed.
 
-[na,ma] = size(a);
-[nb,mb] = size(b);
-if( (na ~= ma) | (na ~= nb) | (nb ~= mb) ) 
-  disp('qz: incompatible dimensions');
-  return;
-endif
+  if (nargin != 2)
+    error ("usage: [aa, bb, q, z] = qzhess (a, b)");
+  endif
+
+  [na, ma] = size (a);
+  [nb, mb] = size (b);
+  if (na != ma || na != nb || nb != mb)
+    error ("qzhess: incompatible dimensions");
+  endif
+
+# Reduce to hessenberg-triangular form.
 
-% reduce to hessenberg-triangular form
-[q,bb] = qr(b);
-aa = q'*a;
-q = q';
-z = eye(na);
-for j=1:(na-2)
-  for i=na:-1:(j+2)
-    %disp(['zero out aa(',num2str(i),',',num2str(j),')'])
-    rot= givens(aa(i-1,j),aa(i,j));
-    aa((i-1):i,:) = rot*aa((i-1):i,:);
-    bb((i-1):i,:) = rot*bb((i-1):i,:);
-    q( (i-1):i,:) = rot*q( (i-1):i,:);
-    
-    %disp(['now zero out bb(',num2str(i),',',num2str(i-1),')'])
-    rot = givens(bb(i,i),bb(i,i-1))';
-    bb(:,(i-1):i) = bb(:,(i-1):i)*rot';
-    aa(:,(i-1):i) = aa(:,(i-1):i)*rot';
-    z(:, (i-1):i) = z(:, (i-1):i)*rot';
+  [q, bb] = qr (b);
+  aa = q' * a;
+  q = q';
+  z = eye (na);
+  for j = 1:(na-2)
+    for i = na:-1:(j+2)
+
+# disp (["zero out aa(", num2str(i), ",", num2str(j), ")"])
+
+      rot = givens (aa (i-1, j), aa (i, j));
+      aa ((i-1):i, :) = rot *aa ((i-1):i, :);
+      bb ((i-1):i, :) = rot *bb ((i-1):i, :);
+      q  ((i-1):i, :) = rot *q  ((i-1):i, :);
+
+# disp (["now zero out bb(", num2str(i), ",", num2str(i-1), ")"])
+
+      rot = givens (bb (i, i), bb (i, i-1))';
+      bb (:, (i-1):i) = bb (:, (i-1):i) * rot';
+      aa (:, (i-1):i) = aa (:, (i-1):i) * rot';
+      z  (:, (i-1):i) = z  (:, (i-1):i) * rot';
+
+    endfor
   endfor
-endfor
 
-bb(2,1) = 0;
-for i=3:na
-  bb(i,1:(i-1)) = zeros(1,i-1);
-  aa(i,1:(i-2)) = zeros(1,i-2);
-endfor
+  bb (2, 1) = 0.0;
+  for i = 3:na
+    bb (i, 1:(i-1)) = zeros (1, i-1);
+    aa (i, 1:(i-2)) = zeros (1, i-2);
+  endfor
+
 endfunction