-
Notifications
You must be signed in to change notification settings - Fork 1
/
mp_myqr.m
98 lines (91 loc) · 2.42 KB
/
mp_myqr.m
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
function [Q,R]=mp_myqr(A,flag)
%QR Computes QR decomposition of A.
% It uses the Gram-Schmidt algorithm.
% The inner product is given by <x,y>=x^ty.
% The command is [Q,R]=myqr(A,flag).
% If flag=0, the "economy size" answer is provided
% based upon GRAM.M, Copyright 1993 Terry Lawson
% Terry Lawson, Math Department, Tulane University, 11/93
% Modified by CLV, 20061112
% Nows deals correctly with rank deficient matrix A
% Modified by CLV, 20041229
% "flag" is added as a parameter, to specify the "economy size" decomposition or not
if (nargin < 2), flag = 1;end %request full size decomposition
[m,n]=size(A);
zero=A(1)*0; %a zero of the same data type as matrix A
switch class(A)
case 'double'
EPS=eps;
otherwise
EPS=eps(zero+1);
end
tol = max([m,n])*EPS*norm(A,'inf');
v=A;w=A;
%R=zeros(zero,m,n);
R=zero*zeros(m,n);
k=1;
w(:,1)=v(:,1);
R(1,1)=sqrt(w(:,1)'*w(:,1));
Q(:,1)=w(:,1)/R(1,1);nn(1)=R(1,1);
for j= 2:n
wtemp = v(:,j);
for i= 1:k
t(i,j)=v(:,j)'*w(:,i)/(w(:,i)'*w(:,i));
wtemp = wtemp-t(i,j)*w(:,i);
R(i,j)=t(i,j)*nn(i);
end
if sqrt(wtemp'*wtemp) > tol
k=k+1;
w(:,k)=wtemp;
nn(k)=sqrt(wtemp'*wtemp);
R(k,j)=nn(k);
Q(:,k)=w(:,k)/nn(k);
else
k=k+1;
[Q,nn,w,k]=fill_with_something(Q,nn,w,m,n,k,tol,zero);
nn(k)=0;
end
end
%exit if flag is explicitly set to zero
if flag==0;
R=R(1:k,:);
return
end
%if m>n, we should invent something for the rest of the columns
k=n+1;
while k<=m
[Q,nn,w,k]=fill_with_something(Q,nn,w,m,n,k,tol,zero);
k=k+1;
% dummy=rand(m,1)+zero;
% for i=1:(k-1)
% rik=Q(:,i)'*dummy;
% dummy=dummy-rik*Q(:,i);
% end
% qk=dummy;
% rkk=sqrt(qk'*qk);
% if rkk<tol
% %you are really unlucky! the random vector is collinear with all the previous
% %ones. But, do not worry! We can save your life doing nothing!
% else
% Q(:,k)=qk/rkk;
% k=k+1;
% end
end
%keyboard
R(m,n)=0;
return
function [Q,nn,w,k]=fill_with_something(Q,nn,w,m,n,k,tol,zero);
rkk=0;
while rkk<tol
dummy=rand(m,1)+zero;
for i=1:(k-1)
rik=Q(:,i)'*dummy;
dummy=dummy-rik*Q(:,i);
end
qk=dummy;
rkk=sqrt(qk'*qk);
end
Q(:,k)=qk/rkk;
nn(k)=rkk;
w(:,k)=qk;
%k=k+1;