Mercurial > octave-antonio

annotate scripts/linear-algebra/qzhess.m @ 29:69497d3b3b75

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-10 21:42:45 by jwe] Initial revision
author jwe
date Tue, 10 Aug 1993 21:42:45 +0000
parents
children 98eb51c870b2
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
29
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
1 function [aa,bb,q,z] = qzhess(a,b)
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
2
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
3 % compute the qz decomposition of the matrix pencil (a - lambda b)
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
4 % result: (for MATLAB compatibility):
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
5 % q*a*z = aa, q*b*z = bb
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
6 % q, z orthogonal,
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
7 % v = matrix of generalized eigenvectors
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
8 %
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
9 % This ought to be done in a compiled program
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
10 % Algorithm taken from Golub and Van Loan, MATRIX COMPUTATIONS, 2nd ed.
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
11
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
12 [na,ma] = size(a);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
13 [nb,mb] = size(b);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
14 if( (na ~= ma) | (na ~= nb) | (nb ~= mb) )
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
15 disp('qz: incompatible dimensions');
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
16 return;
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
17 endif
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
18
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
19 % reduce to hessenberg-triangular form
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
20 [q,bb] = qr(b);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
21 aa = q'*a;
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
22 q = q';
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
23 z = eye(na);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
24 for j=1:(na-2)
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
25 for i=na:-1:(j+2)
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
26 %disp(['zero out aa(',num2str(i),',',num2str(j),')'])
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
27 rot= givens(aa(i-1,j),aa(i,j));
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
28 aa((i-1):i,:) = rot*aa((i-1):i,:);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
29 bb((i-1):i,:) = rot*bb((i-1):i,:);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
30 q( (i-1):i,:) = rot*q( (i-1):i,:);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
31
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
32 %disp(['now zero out bb(',num2str(i),',',num2str(i-1),')'])
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
33 rot = givens(bb(i,i),bb(i,i-1))';
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
34 bb(:,(i-1):i) = bb(:,(i-1):i)*rot';
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
35 aa(:,(i-1):i) = aa(:,(i-1):i)*rot';
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
36 z(:, (i-1):i) = z(:, (i-1):i)*rot';
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
37 endfor
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
38 endfor
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
39
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
40 bb(2,1) = 0;
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
41 for i=3:na
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
42 bb(i,1:(i-1)) = zeros(1,i-1);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
43 aa(i,1:(i-2)) = zeros(1,i-2);
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
44 endfor
69497d3b3b75 [project @ 1993-08-10 21:42:45 by jwe]
jwe
parents:
diff changeset
45 endfunction