home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / la / qzhess.m < prev    next >
Text File  |  1999-04-29  |  2KB  |  82 lines
  1. ## Copyright (C) 1996, 1997 John W. Eaton
  2. ##
  3. ## This file is part of Octave.
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it
  6. ## under the terms of the GNU General Public License as published by
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License for more details.
  14. ##
  15. ## You should have received a copy of the GNU General Public License
  16. ## along with Octave; see the file COPYING.  If not, write to the Free
  17. ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  18. ## 02111-1307, USA.
  19.  
  20. ## Usage: [aa, bb, q, z] = qzhess (a, b)
  21. ##
  22. ## Compute the qz decomposition of the matrix pencil (a - lambda b)
  23. ##
  24. ## result: (for Matlab compatibility):
  25. ##
  26. ##   aa = q*a*z and bb = q*b*z, with q, z orthogonal, and
  27. ##   v = matrix of generalized eigenvectors.
  28. ##
  29. ## This ought to be done in a compiled program
  30. ##
  31. ## Algorithm taken from Golub and Van Loan, Matrix Computations, 2nd ed.
  32.  
  33. ## Author: A. S. Hodel <scotte@eng.auburn.edu>
  34. ## Created: August 1993
  35. ## Adapted-By: jwe
  36.  
  37. function [aa, bb, q, z] = qzhess (a, b)
  38.  
  39.   if (nargin != 2)
  40.     error ("usage: [aa, bb, q, z] = qzhess (a, b)");
  41.   endif
  42.  
  43.   [na, ma] = size (a);
  44.   [nb, mb] = size (b);
  45.   if (na != ma || na != nb || nb != mb)
  46.     error ("qzhess: incompatible dimensions");
  47.   endif
  48.  
  49.   ## Reduce to hessenberg-triangular form.
  50.  
  51.   [q, bb] = qr (b);
  52.   aa = q' * a;
  53.   q = q';
  54.   z = eye (na);
  55.   for j = 1:(na-2)
  56.     for i = na:-1:(j+2)
  57.  
  58.       ## disp (["zero out aa(", num2str(i), ",", num2str(j), ")"])
  59.  
  60.       rot = givens (aa (i-1, j), aa (i, j));
  61.       aa ((i-1):i, :) = rot *aa ((i-1):i, :);
  62.       bb ((i-1):i, :) = rot *bb ((i-1):i, :);
  63.       q  ((i-1):i, :) = rot *q  ((i-1):i, :);
  64.  
  65.       ## disp (["now zero out bb(", num2str(i), ",", num2str(i-1), ")"])
  66.  
  67.       rot = givens (bb (i, i), bb (i, i-1))';
  68.       bb (:, (i-1):i) = bb (:, (i-1):i) * rot';
  69.       aa (:, (i-1):i) = aa (:, (i-1):i) * rot';
  70.       z  (:, (i-1):i) = z  (:, (i-1):i) * rot';
  71.  
  72.     endfor
  73.   endfor
  74.  
  75.   bb (2, 1) = 0.0;
  76.   for i = 3:na
  77.     bb (i, 1:(i-1)) = zeros (1, i-1);
  78.     aa (i, 1:(i-2)) = zeros (1, i-2);
  79.   endfor
  80.  
  81. endfunction
  82.