[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: QR updating changeset
From: |
John W. Eaton |
Subject: |
Re: QR updating changeset |
Date: |
Wed, 05 Mar 2008 14:30:57 -0500 |
On 5-Mar-2008, Jaroslav Hajek wrote:
| On Wed, Mar 5, 2008 at 3:52 AM, John W. Eaton <address@hidden> wrote:
| >
| > Finally, I'm seeing this error now (even without my changes, so I
| > don't think it is a problem that I introduced):
| >
| > ====>>>>> processing /tmp/jwe/octave/src/DLD-FUNCTIONS/qrinsert.cc
| > ***** test
| > A = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
| > 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
| > 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
| > 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
| > 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
| >
| > x = [0.20351 + 0.05401i;
| > 0.13141 + 0.43708i;
| > 0.29808 + 0.08789i;
| > 0.69821 + 0.38844i;
| > 0.74871 + 0.25821i ];
| >
| > [Q,R] = qr(A);
| > [Q,R] = qrinsert(Q,R,3,x);
| > assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
| > assert(norm(vec(triu(R)-R),Inf) == 0)
| > assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps)
| >
| > !!!!! test failed
| > error: assert (norm (vec (Q * R - [A(:, 1:2), x, A(:, 3)]), Inf) < norm
(A) * 1e1 * eps) failed
| >
| > Your patch and my changes are all in my public archive now, so will
| > you please update and take a look?
| >
| > Thanks,
| >
| > jwe
| >
| >
| >
|
| John,
| I can't reproduce this failure, all tests are passing when I compile.
| Would you please check again with the newest changesets? (Submitted in
| the more recent thread).
| In case it persists, please run the attached script (the failing test
| but more verbosely) so that I can see whether this is a numerical
| issue. The difference 1e1*eps is an ad hoc value though it seems more
| than enough on my machine.
Here is what I see:
octave:1> more off
octave:2> A = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
> 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
> 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
> 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
> 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
octave:3>
octave:3> x = [0.20351 + 0.05401i;
> 0.13141 + 0.43708i;
> 0.29808 + 0.08789i;
> 0.69821 + 0.38844i;
> 0.74871 + 0.25821i ];
octave:4>
octave:4> [Q,R] = qr(A)
Q =
Columns 1 through 3:
-0.3126136 - 0.4821956i 0.1427138 - 0.1060751i -0.3019292 - 0.3902357i
-0.2968279 - 0.3317875i -0.2376419 + 0.1795628i 0.0067884 + 0.1550643i
-0.0467395 - 0.1741870i 0.7862817 + 0.1492957i 0.0081376 + 0.5375256i
-0.4595937 - 0.3633142i -0.2964661 + 0.2030497i 0.4531435 + 0.1829374i
-0.0568631 - 0.3042824i 0.3316416 - 0.0053703i -0.0494593 - 0.4496805i
Columns 4 and 5:
-0.1359284 + 0.1721766i -0.5264857 + 0.2634149i
-0.6359470 - 0.4479185i 0.2889609 + 0.0210917i
0.0639545 - 0.1589105i -0.0900039 + 0.0224230i
0.5001559 + 0.1207159i -0.0700878 - 0.1385980i
0.1954472 + 0.1206172i 0.7312090 - 0.0457203i
R =
-1.98457 + 0.00000i -0.54697 + 0.32204i -1.46847 + 0.04324i
0.00000 + 0.00000i 1.07725 + 0.00000i 0.88459 + 0.08638i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.70444 + 0.00000i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
octave:5> [Q,R] = qrinsert(Q,R,3,x)
Q =
Columns 1 through 3:
-0.312614 - 0.482196i 0.142714 - 0.106075i -0.474593 + 0.072920i
-0.296828 - 0.331788i -0.237642 + 0.179563i -0.084010 + 0.293504i
-0.046739 - 0.174187i 0.786282 + 0.149296i -0.116315 + 0.258789i
-0.459594 - 0.363314i -0.296466 + 0.203050i 0.323550 - 0.191302i
-0.056863 - 0.304282i 0.331642 - 0.005370i 0.568194 - 0.362790i
Columns 4 and 5:
-0.463765 + 0.224506i -0.376576 - 0.017169i
-0.026834 - 0.056892i 0.707904 - 0.338831i
0.414950 - 0.242240i 0.053897 + 0.112589i
0.208323 - 0.427587i -0.291584 + 0.274088i
-0.135596 + 0.505057i 0.241764 + 0.088233i
R =
-1.98457 + 0.00000i -0.54697 + 0.32204i -1.46847 + 0.04324i
0.00000 + 0.00000i 1.07725 + 0.00000i 0.88459 + 0.08638i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.70444 + 0.00000i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
octave:5> [Q,R] = qrinsert(Q,R,3,x)
Q =
Columns 1 through 3:
-0.312614 - 0.482196i 0.142714 - 0.106075i -0.474593 + 0.072920i
-0.296828 - 0.331788i -0.237642 + 0.179563i -0.084010 + 0.293504i
-0.046739 - 0.174187i 0.786282 + 0.149296i -0.116315 + 0.258789i
-0.459594 - 0.363314i -0.296466 + 0.203050i 0.323550 - 0.191302i
-0.056863 - 0.304282i 0.331642 - 0.005370i 0.568194 - 0.362790i
Columns 4 and 5:
-0.463765 + 0.224506i -0.376576 - 0.017169i
-0.026834 - 0.056892i 0.707904 - 0.338831i
0.414950 - 0.242240i 0.053897 + 0.112589i
0.208323 - 0.427587i -0.291584 + 0.274088i
-0.135596 + 0.505057i 0.241764 + 0.088233i
R =
Columns 1 through 3:
-1.98457 + 0.00000i -0.54697 + 0.32204i -1.05169 - 0.11184i
0.00000 + 0.00000i 1.07725 + 0.00000i 0.55728 - 0.02241i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.49600 + 0.47446i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
Column 4:
-1.46847 + 0.04324i
0.88459 + 0.08638i
0.38363 + 0.00000i
-0.20263 + 0.55498i
0.00000 + 0.00000i
octave:6> norm(vec(Q'*Q - eye(5)),Inf)
ans = 2.2348e-16
octave:7> norm(vec(triu(R)-R),Inf)
ans = 0
octave:8> norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf)
ans = 0.27774
octave:9> disp(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
1
octave:10> disp(norm(vec(triu(R)-R),Inf) == 0)
1
octave:11> disp(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps)
0
octave:12> Q*R
ans =
0.620405 + 0.956953i 0.480013 + 0.048806i 0.082001 + 0.290764i
0.402627 + 0.338171i
0.589077 + 0.658457i 0.013205 + 0.279323i -0.034270 + 0.593248i
0.229284 + 0.721929i
0.092758 + 0.345687i 0.928679 + 0.241052i 0.290722 + 0.327169i
0.764536 + 0.832406i
0.912098 + 0.721024i 0.049018 + 0.269452i 0.533303 + 0.611921i
0.730029 + 0.796517i
0.112849 + 0.603871i 0.486352 + 0.142337i 0.664426 + 0.405587i
0.355646 + 0.151496i
jwe