[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: polyvalm problems for defective matrix input
|
From: |
Rolf Fabian |
|
Subject: |
Re: polyvalm problems for defective matrix input |
|
Date: |
Mon, 4 Feb 2008 08:30:42 -0800 (PST) |
Many many thanks to Ben Abbott who transfered the 'polyvalmh_demo.m'
test script to MatLab and tested and ran it also under ML ( his results
are labelled 'MATLAB RESULT' in the overview below.
In summary all the matlab results for the examples given by Ben
confirm the results returned by proposed 'polyvalmh.m' and discard
all the results from current Octave V3.0.0 'polyvalm.m'.
The most outstanding difference can be seen in EXA #3,
where MatLab and 'polyvalmh.m' return a value of '5'
for the upper right corner element, whereas Octave's
polyvalm returns 1.8014e+016 for that element, which is
more than only a bit different.
EXA #6 shows that a non-trivial, non-nilpotent, regular but 'defective' x
may also lead to subtle large differences between Octave's
'polyvam.m' result and both 'polyvalmh.m' and MatLab polyvalm.
Again MatLab confirms the result of 'polyvalmh'.
The MatLab results strongly support my suggestion to replace
Octave's current 'polyvalm.m' by a renamed 'polyvalmh.m'.
I already pointed out that you can do this PROVIDED that my
name appears as author in the copyright notice.
You might wonder whether I do insist so much on this point ?
The answer is very simple. Actually it has been me who has written most
of the code of the current Octave 'polyvalm.m' function.
The original version by Tony Richardson (dated 1994) failed horribly without
warning
for even the most simple matrices. It was essentially a '3 line thing'
defining only a
unitary diagonalization algorithm to do polyvalm. Thus it was actually only
applicable
to the special case of hermitian (more precisely 'normal) input matrices and
persistently
failed for the more common general case, which requires a general (not only
unitary)
similarity transformation algorithm. Thus I decided on Aug 19, 1997 to send
a bug report
<http://velveeta.che.wisc.edu/octave/lists/archive//bug-octave.1997/msg00362.html>
where I suggested a complete revision of 'polyvalm' for Octave 2.0.9 . The
reply
by jwe dated Aug. 28, 2997
<http://velveeta.che.wisc.edu/octave/lists/archive//bug-octave.1997/msg00388.html
accepts the proposed function.
If you have a look at the 1st URL you'll easily see that the current Octave
V3.0.0
'polyvalm.m' corresponds nearly exactly character-by-character to my
proposed
function from Aug 19, 1997. Besides some changes in the help part
(texinfo-related)
The only essential difference is that you simply removed my autorship notice
which I've written aside to the authorship notice of Tony (check my
authorship
notice in the URL).
This treatment wasn't really polite and even though I don't know much about
licenses I
wonder whether such a policy can be conform to the terms of the GNU public
license.
I'm pretty much shure that such a policy concerning authorship on free
software is
discouraging other contributors also, not only me.
Now that I've also solved the even more sophisticated problem of handling
defective matrix input correctly the function doesn't have to do anything
anymore
with the '3-lines-code' of Tony.
My offer : You can have the code, but I insist on getting single authorship.
Rolf Fabian
=======================================================
OVERVIEW TEST DEMOS
=======================================================
octave-3.0.0.exe:1> [y_Octave,y_polyvalmh]=polyvalmh_demo(1)
EXA #1
the 'classical', most simple nilpotent + defective matrix x [note: x^ (k >
2) = 0]
x =
0 1
0 0
p =
1 2 3
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 5e-293.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
3 0
0 3
y_polyvalmh =
3 2
0 3
MATLAB RESULT
:> polyvalm(p,x)
ans =
3 2
0 3
=======================================================
octave-3.0.0.exe:2> [y_Octave,y_polyvalmh]=polyvalmh_demo(2)
EXA #2
similar to #1, but complex input matrix
x =
0 + 0i 1 + 1i
0 + 0i 0 + 0i
p =
1 2 3
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 5e-293.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
3 0
0 3
y_polyvalmh =
3 + 0i 2 + 2i
0 + 0i 3 + 0i
MATLAB RESULT
:> polyvalm(p,x)
ans =
3.0000 2.0000 + 2.0000i
0 3.0000
=======================================================
octave-3.0.0.exe:3> [y_Octave,y_polyvalmh]=polyvalmh_demo(3)
EXA #3
simple defective regular matrix.
p =
1 2 3
x =
1 1 1
0 1 1
0 0 1
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 2.5e-032.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
6.0000e+000 0.0000e+000 1.8014e+016
0.0000e+000 6.0000e+000 0.0000e+000
0.0000e+000 0.0000e+000 6.0000e+000
y_polyvalmh =
6 4 5
0 6 4
0 0 6
MATLAB RESULT
:> polyvalm(p,x)
ans =
6 4 5
0 6 4
0 0 6
=======================================================
octave-3.0.0.exe:4> [y_Octave,y_polyvalmh]=polyvalmh_demo(4)
EXA #4
same matrix as #3, scaled to very small norm.
p =
1 2 3
scalefac = 1.0000e-060
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 2.5e-032.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
3.0000e+000 0.0000e+000 9.0072e+015
0.0000e+000 3.0000e+000 0.0000e+000
0.0000e+000 0.0000e+000 3.0000e+000
y_polyvalmh =
3.00000 0.00000 0.00000
0.00000 3.00000 0.00000
0.00000 0.00000 3.00000
MATLAB RESULT
:> polyvalm(p,x)
ans =
3.0000 0.0000 0.0000
0 3.0000 0.0000
0 0 3.0000
=======================================================
octave-3.0.0.exe:5> [y_Octave,y_polyvalmh]=polyvalmh_demo(6)
EXA #6
a non-trivial, non-nilpotent, regular (det(x)=1) but 'defective' x
x =
0 1 0 0
1 0 2 3
0 0 0 1
0 0 1 0
p =
2 3 4 5 6 7 8 9 10 11
{message from Octave for 'polyvalm(p,x)'}
warning: matrix singular to machine precision, rcond = 2.22045e-017
warning: attempting to find minimum norm solution
warning: dgelsd: rank deficient 4x4 matrix, rank = 3
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 3.7e-017.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
3.5000e+001 3.2527e+001 5.2000e+001 -2.5000e+001
3.0000e+001 3.5355e+001 5.6000e+001 -2.8000e+001
-1.9984e-015 -3.2658e-015 5.3975e+000 -2.9996e+000
-8.8818e-016 -7.5364e-016 -2.7014e+000 1.5013e+000
y_polyvalmh =
35.00000 30.00000 220.00000 230.00000
30.00000 35.00000 290.00000 310.00000
0.00000 0.00000 35.00000 30.00000
0.00000 0.00000 30.00000 35.00000
MATLAB RESULT
:> polyvalm(p,x)
ans =
35 30 220 230
30 35 290 310
0 0 35 30
0 0 30 35
=======================================================
octave-3.0.0.exe:6> [y_Octave,y_polyvalmh]=polyvalmh_demo(8)
EXA #6
matrix from example #6 scaled to very small norm.
x =
0 1 0 0
1 0 2 3
0 0 0 1
0 0 1 0
p =
2 3 4 5 6 7 8 9 10 11
scalefac = 1.0000e-040
{message from Octave for 'polyvalm(p,x)'}
warning: matrix singular to machine precision, rcond = 2.22045e-017
warning: attempting to find minimum norm solution
warning: dgelsd: rank deficient 4x4 matrix, rank = 3
warning: polyvalmh: numeric underflow probable. Result should be checked.
warning: polyvalmh: found 'defective' input: rcond (eigensystem) = 3.7e-017.
polyvalmh: 'Horner algorithm' used for 'defective' input.
y_Octave =
1.1000e+001 2.4749e+000 1.2000e+001 -1.7500e+000
-4.4409e-016 8.1317e+000 4.0000e+000 1.0000e+000
-1.4655e-015 -2.8262e-015 8.9933e+000 -4.1586e-001
9.7700e-016 2.1981e-015 -8.4660e+000 3.9148e-001
y_polyvalmh =
11.00000 0.00000 0.00000 0.00000
0.00000 11.00000 0.00000 0.00000
0.00000 0.00000 11.00000 0.00000
0.00000 0.00000 0.00000 11.00000
MATLAB RESULT
:> polyvalm(p,x)
ans =
11.0000 0.0000 0.0000 0.0000
0.0000 11.0000 0.0000 0.0000
0 0 11.0000 0.0000
0 0 0.0000 11.0000
=======================================================
-----
Rolf Fabian
<r dot fabian at jacobs-university dot de>
--
View this message in context:
http://www.nabble.com/polyvalm-problems-for-defective-matrix-input-tp15183200p15269047.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.