Mercurial > matrix-functions
annotate matrixcomp/gqr.m @ 0:8f23314345f4 draft
Create local repository for matrix toolboxes. Step #0 done.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Wed, 06 May 2015 14:56:53 +0200 |
parents | |
children |
rev | line source |
---|---|
0
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
1 function [U, Q, L, S] = gqr(A, B, partial) |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
2 %GQR Generalized QR factorization. |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
3 % [U, Q, L, S] = GQR(A, B, partial) factorizes |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
4 % the m-by-n A and p-by-n B, where m >= n >= p, as |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
5 % A = U*L*Q^T, B = S*Q^T, with U and Q orthogonal |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
6 % and L = [0; L1], S = [S1 0], with L1 and S1 lower triangular. |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
7 % If a nonzero third argument is present then only a partial reduction |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
8 % of A is performed: the first p columns of A are not reduced to |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
9 % triangular form (which is sufficient for solving the LSE problem). |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
10 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
11 % Reference: |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
12 % N. J. Higham, Accuracy and Stability of Numerical Algorithms, |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
13 % Second edition, Society for Industrial and Applied Mathematics, |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
14 % Philadelphia, PA, 2002; sec. 20.9. |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
15 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
16 [m n] = size(A); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
17 [p n1] = size(B); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
18 if nargin < 3, partial = 0; end |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
19 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
20 if n ~= n1, error('A and B must have same number of columns!'), end |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
21 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
22 if partial |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
23 limit = p+1; |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
24 else |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
25 limit = 1; |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
26 end |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
27 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
28 [Q, S] = qr(B'); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
29 S = S'; |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
30 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
31 U = eye(m); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
32 A = A*Q; |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
33 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
34 % QL factorization of AQ. |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
35 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
36 for i = n:-1:limit |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
37 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
38 % Vector-reversal so that Householder eliminates leading |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
39 % rather than trailing elements. |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
40 temp = A(1:m-n+i,i); temp = temp(end:-1:1); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
41 [v, beta] = gallery('house',temp); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
42 v = v(end:-1:1); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
43 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
44 temp = A(1:m-n+i,1:i); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
45 A(1:m-n+i,1:i) = temp - beta*v*(v'*temp); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
46 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
47 % Put zeros where they're supposed to be! |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
48 A(1:m-n+i-1,i) = zeros(m-n+i-1,1); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
49 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
50 temp = U(:,1:m-n+i); |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
51 U(:,1:m-n+i) = temp - beta*temp*v*v'; |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
52 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
53 end |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
54 |
8f23314345f4
Create local repository for matrix toolboxes. Step #0 done.
Antonio Pino Robles <data.script93@gmail.com>
parents:
diff
changeset
|
55 L = A; |