Mercurial > matrix-functions

comparison matrixcomp/gqr.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [U, Q, L, S] = gqr(A, B, partial)
2 %GQR Generalized QR factorization.
3 % [U, Q, L, S] = GQR(A, B, partial) factorizes
4 % the m-by-n A and p-by-n B, where m >= n >= p, as
5 % A = U*L*Q^T, B = S*Q^T, with U and Q orthogonal
6 % and L = [0; L1], S = [S1 0], with L1 and S1 lower triangular.
7 % If a nonzero third argument is present then only a partial reduction
8 % of A is performed: the first p columns of A are not reduced to
9 % triangular form (which is sufficient for solving the LSE problem).
10
11 % Reference:
12 % N. J. Higham, Accuracy and Stability of Numerical Algorithms,
13 % Second edition, Society for Industrial and Applied Mathematics,
14 % Philadelphia, PA, 2002; sec. 20.9.
15
16 [m n] = size(A);
17 [p n1] = size(B);
18 if nargin < 3, partial = 0; end
19
20 if n ~= n1, error('A and B must have same number of columns!'), end
21
22 if partial
23 limit = p+1;
24 else
25 limit = 1;
26 end
27
28 [Q, S] = qr(B');
29 S = S';
30
31 U = eye(m);
32 A = A*Q;
33
34 % QL factorization of AQ.
35
36 for i = n:-1:limit
37
38 % Vector-reversal so that Householder eliminates leading
39 % rather than trailing elements.
40 temp = A(1:m-n+i,i); temp = temp(end:-1:1);
41 [v, beta] = gallery('house',temp);
42 v = v(end:-1:1);
43
44 temp = A(1:m-n+i,1:i);
45 A(1:m-n+i,1:i) = temp - beta*v*(v'*temp);
46
47 % Put zeros where they're supposed to be!
48 A(1:m-n+i-1,i) = zeros(m-n+i-1,1);
49
50 temp = U(:,1:m-n+i);
51 U(:,1:m-n+i) = temp - beta*temp*v*v';
52
53 end
54
55 L = A;