Mercurial > octave-nkf

comparison scripts/linear-algebra/qzhess.m @ 54:98eb51c870b2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-11 21:29:50 by jwe]
author jwe
date Wed, 11 Aug 1993 21:30:43 +0000
parents 69497d3b3b75
children 2d480148756b
comparison
equal deleted inserted replaced
53:4565ad8b4697 54:98eb51c870b2
1 function [aa,bb,q,z] = qzhess(a,b) 1 function [aa, bb, q, z] = qzhess (a, b)
2 2
3 % compute the qz decomposition of the matrix pencil (a - lambda b) 3 # Compute the qz decomposition of the matrix pencil (a - lambda b)
4 % result: (for MATLAB compatibility): 4 #
5 % q*a*z = aa, q*b*z = bb 5 # result: (for Matlab compatibility):
6 % q, z orthogonal, 6 #
7 % v = matrix of generalized eigenvectors 7 # aa = q*a*z and bb = q*b*z, with q, z orthogonal, and
8 % 8 # v = matrix of generalized eigenvectors.
9 % This ought to be done in a compiled program 9 #
10 % Algorithm taken from Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed. 10 # This ought to be done in a compiled program
11 #
12 # Algorithm taken from Golub and Van Loan, Matrix Computations, 2nd ed.
11 13
12 [na,ma] = size(a); 14 if (nargin != 2)
13 [nb,mb] = size(b); 15 error ("usage: [aa, bb, q, z] = qzhess (a, b)");
14 if( (na ~= ma) | (na ~= nb) | (nb ~= mb) ) 16 endif
15 disp('qz: incompatible dimensions');
16 return;
17 endif
18 17
19 % reduce to hessenberg-triangular form 18 [na, ma] = size (a);
20 [q,bb] = qr(b); 19 [nb, mb] = size (b);
21 aa = q'*a; 20 if (na != ma || na != nb || nb != mb)
22 q = q'; 21 error ("qzhess: incompatible dimensions");
23 z = eye(na); 22 endif
24 for j=1:(na-2) 23
25 for i=na:-1:(j+2) 24 # Reduce to hessenberg-triangular form.
26 %disp(['zero out aa(',num2str(i),',',num2str(j),')']) 25
27 rot= givens(aa(i-1,j),aa(i,j)); 26 [q, bb] = qr (b);
28 aa((i-1):i,:) = rot*aa((i-1):i,:); 27 aa = q' * a;
29 bb((i-1):i,:) = rot*bb((i-1):i,:); 28 q = q';
30 q( (i-1):i,:) = rot*q( (i-1):i,:); 29 z = eye (na);
31 30 for j = 1:(na-2)
32 %disp(['now zero out bb(',num2str(i),',',num2str(i-1),')']) 31 for i = na:-1:(j+2)
33 rot = givens(bb(i,i),bb(i,i-1))'; 32
34 bb(:,(i-1):i) = bb(:,(i-1):i)*rot'; 33 # disp (["zero out aa(", num2str(i), ",", num2str(j), ")"])
35 aa(:,(i-1):i) = aa(:,(i-1):i)*rot'; 34
36 z(:, (i-1):i) = z(:, (i-1):i)*rot'; 35 rot = givens (aa (i-1, j), aa (i, j));
36 aa ((i-1):i, :) = rot *aa ((i-1):i, :);
37 bb ((i-1):i, :) = rot *bb ((i-1):i, :);
38 q ((i-1):i, :) = rot *q ((i-1):i, :);
39
40 # disp (["now zero out bb(", num2str(i), ",", num2str(i-1), ")"])
41
42 rot = givens (bb (i, i), bb (i, i-1))';
43 bb (:, (i-1):i) = bb (:, (i-1):i) * rot';
44 aa (:, (i-1):i) = aa (:, (i-1):i) * rot';
45 z (:, (i-1):i) = z (:, (i-1):i) * rot';
46
47 endfor
37 endfor 48 endfor
38 endfor
39 49
40 bb(2,1) = 0; 50 bb (2, 1) = 0.0;
41 for i=3:na 51 for i = 3:na
42 bb(i,1:(i-1)) = zeros(1,i-1); 52 bb (i, 1:(i-1)) = zeros (1, i-1);
43 aa(i,1:(i-2)) = zeros(1,i-2); 53 aa (i, 1:(i-2)) = zeros (1, i-2);
44 endfor 54 endfor
55
45 endfunction 56 endfunction