[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: QR updating changeset
From: |
Jaroslav Hajek |
Subject: |
Re: QR updating changeset |
Date: |
Thu, 6 Mar 2008 08:09:13 +0100 |
On Wed, Mar 5, 2008 at 8:30 PM, John W. Eaton <address@hidden> wrote:
>
> 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
>
Thanks for the verbose info, it helped me to locate the error quickly.
The attached changeset should fix it.
--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
qrincfix.diff
Description: Text document