1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
|
function [H,R] = myqr (A)
%MYQR QR factorization using Householder reflections
% uses function [v,beta,xnorm] = hmake1 (x)
% and function hx = happly (v, beta, x)
%
% Example
% [H,R] = myqr (A)
% See also: testall
% CXSparse, Copyright (c) 2006-2022, Timothy A. Davis. All Rights Reserved.
% SPDX-License-Identifier: LGPL-2.1+
[m n] = size (A) ;
H = zeros (m,n) ;
R = zeros (m,n) ;
for k = 1:n
% apply prior H's
% fprintf ('\n-----------------init %d\n', k) ;
x = A (:,k) ;
for i = 1:k-1
v = H(((i+1):m),i) ;
v = [1 ; v] ; %#ok
beta = H (i,i) ;
% n1 = norm (x (i:m)) ;
x (i:m) = happly (v, beta, x (i:m)) ;
% n2 = norm (x (i:m)) ;
% fprintf ('=============== i %d %g %g\n', i, n1, n2) ;
% beta
% v'
% X = x'
% pause
% i
% x
end
% k
% x
% make Hk
% fprintf ('x(k:m) = ') ; x (k:m)
[v,beta,xnorm] = hmake1 (x (k:m)) ;
H (k,k) = beta ;
H (k+1:m, k) = v (2:end) ;
R (1:(k-1),k) = x (1:(k-1)) ;
R (k,k) = xnorm ;
% full (R)
% pause
end
% s2 = svd (full (R)) ;
% [s1 s2 s1-s2]
|
