Mercurial > octave

comparison scripts/linear-algebra/qzhess.m @ 2303:5cffc4b8de57

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-06-24 09:15:24 by jwe]
author jwe
date Mon, 24 Jun 1996 09:15:24 +0000
parents 5d29638dd524
children 2b5788792cad
comparison
equal deleted inserted replaced
2302:470c856bf55a 2303:5cffc4b8de57
1 # Copyright (C) 1996 John W. Eaton 1 ### Copyright (C) 1996 John W. Eaton
2 # 2 ###
3 # This file is part of Octave. 3 ### This file is part of Octave.
4 # 4 ###
5 # Octave is free software; you can redistribute it and/or modify it 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 the 6 ### under the terms of the GNU General Public License as published by
7 # Free Software Foundation; either version 2, or (at your option) any 7 ### the Free Software Foundation; either version 2, or (at your option)
8 # later version. 8 ### any later version.
9 # 9 ###
10 # Octave is distributed in the hope that it will be useful, but WITHOUT 10 ### Octave is distributed in the hope that it will be useful, but
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 11 ### WITHOUT ANY WARRANTY; without even the implied warranty of
12 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 12 ### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 # for more details. 13 ### General Public License for more details.
14 # 14 ###
15 # You should have received a copy of the GNU General Public License 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 16 ### along with Octave; see the file COPYING. If not, write to the Free
17 # Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 17 ### Software Foundation, 59 Temple Place - Suite 330, Boston, MA
18 ### 02111-1307, USA.
18 19
19 function [aa, bb, q, z] = qzhess (a, b) 20 function [aa, bb, q, z] = qzhess (a, b)
20 21
21 # Usage: [aa, bb, q, z] = qzhess (a, b) 22 ## Usage: [aa, bb, q, z] = qzhess (a, b)
22 # 23 ##
23 # Compute the qz decomposition of the matrix pencil (a - lambda b) 24 ## Compute the qz decomposition of the matrix pencil (a - lambda b)
24 # 25 ##
25 # result: (for Matlab compatibility): 26 ## result: (for Matlab compatibility):
26 # 27 ##
27 # aa = q*a*z and bb = q*b*z, with q, z orthogonal, and 28 ## aa = q*a*z and bb = q*b*z, with q, z orthogonal, and
28 # v = matrix of generalized eigenvectors. 29 ## v = matrix of generalized eigenvectors.
29 # 30 ##
30 # This ought to be done in a compiled program 31 ## This ought to be done in a compiled program
31 # 32 ##
32 # Algorithm taken from Golub and Van Loan, Matrix Computations, 2nd ed. 33 ## Algorithm taken from Golub and Van Loan, Matrix Computations, 2nd ed.
33 34
34 # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993. 35 ## Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
35 36
36 if (nargin != 2) 37 if (nargin != 2)
37 error ("usage: [aa, bb, q, z] = qzhess (a, b)"); 38 error ("usage: [aa, bb, q, z] = qzhess (a, b)");
38 endif 39 endif
39 40
41 [nb, mb] = size (b); 42 [nb, mb] = size (b);
42 if (na != ma || na != nb || nb != mb) 43 if (na != ma || na != nb || nb != mb)
43 error ("qzhess: incompatible dimensions"); 44 error ("qzhess: incompatible dimensions");
44 endif 45 endif
45 46
46 # Reduce to hessenberg-triangular form. 47 ## Reduce to hessenberg-triangular form.
47 48
48 [q, bb] = qr (b); 49 [q, bb] = qr (b);
49 aa = q' * a; 50 aa = q' * a;
50 q = q'; 51 q = q';
51 z = eye (na); 52 z = eye (na);
52 for j = 1:(na-2) 53 for j = 1:(na-2)
53 for i = na:-1:(j+2) 54 for i = na:-1:(j+2)
54 55
55 # disp (["zero out aa(", num2str(i), ",", num2str(j), ")"]) 56 ## disp (["zero out aa(", num2str(i), ",", num2str(j), ")"])
56 57
57 rot = givens (aa (i-1, j), aa (i, j)); 58 rot = givens (aa (i-1, j), aa (i, j));
58 aa ((i-1):i, :) = rot *aa ((i-1):i, :); 59 aa ((i-1):i, :) = rot *aa ((i-1):i, :);
59 bb ((i-1):i, :) = rot *bb ((i-1):i, :); 60 bb ((i-1):i, :) = rot *bb ((i-1):i, :);
60 q ((i-1):i, :) = rot *q ((i-1):i, :); 61 q ((i-1):i, :) = rot *q ((i-1):i, :);
61 62
62 # disp (["now zero out bb(", num2str(i), ",", num2str(i-1), ")"]) 63 ## disp (["now zero out bb(", num2str(i), ",", num2str(i-1), ")"])
63 64
64 rot = givens (bb (i, i), bb (i, i-1))'; 65 rot = givens (bb (i, i), bb (i, i-1))';
65 bb (:, (i-1):i) = bb (:, (i-1):i) * rot'; 66 bb (:, (i-1):i) = bb (:, (i-1):i) * rot';
66 aa (:, (i-1):i) = aa (:, (i-1):i) * rot'; 67 aa (:, (i-1):i) = aa (:, (i-1):i) * rot';
67 z (:, (i-1):i) = z (:, (i-1):i) * rot'; 68 z (:, (i-1):i) = z (:, (i-1):i) * rot';