| Hmmm...
I had a look at eigens, and it appears to be an algorithmic problem, not a simple coding error in the wrapper. For example, the 2x2 matrix pdl([4,2],[2,0]) yields correct answers, as does the 3x3 matrix ([4,2,0],[2,0,0],[0,0,0]). The related matrix ([4,2,0],[2,0,0],[0,0,1]) yields correct results as well (it is really a decoupled 2x2 and 1x1). Adding cross terms, as in ([4,2,0],[2,0,0.5],[0,0.3,1]), breaks all the eigenvectors.
Sorry to have let folks down here -- SSL seemed to be a pretty stable package and I inserted it into the code without sufficient additional testing.
I just spent a few minutes browsing around the newer GSL documentation, and it seems to be the way to go. For our overhaul, perhaps we should link into the GSL more extensively -- it has become the de facto industry standard for tested, fast code...
Cheers, Craig
On Jul 8, 2006, at 12:50 PM, Bill Coffman wrote: Funny you should mention ...
The eigens for assymmetric matrices seems to be completely broken. It seems to only produce correct results for 2x2 matrices. I implemented code that detects a symmetric matrix, and passes to the eigens_sym program, and that seems to work pretty well. But more extensive testing shows that results for asymmetric matrixes are consistently wrong. I also found that the PDL::Slatec::eigsys function is also wrong. It produces correct eigenvalues (for symmetric matrices at least), but the eigenvectors are wrong.
Here's a simple test:
my ($v,$w) = eigens $A;
my ($d) = $v->dims;
for my $j (0 .. $d-1) {
my $prod1 = $z x $v->index($j)->transpose;
my $prod2 = $w->index($j) * $v->index($j)->transpose;
printf " dim %d error %.5f\n", $j, sum(abs($prod1-$prod2));
}
The good news is that it seems to work great for symmetric matrices (not the slatec one).
If I check in my code, all tests will pass, but it's still broken.
-Bill Coffman
On 7/8/06, Chris Marshall <[EMAIL PROTECTED]> wrote: I'm in the process of re-running the full PDL
install process with cygwin/XP to verify operation
for the latest CVS code and to collect cygwin
specific notes for the upcoming bug-fix release.
Starting with a default build from the current
cvs code, all tests pass except for test #20 in
t/matrixops.t which fails with the following
message:
>> Failure in hqr2 function. Do not trust the given eigenvectors and -values
>> eigensum for the 8x8: nan
>> not ok 20
>> # Failed test 20 in t/matrixops.t at line 118
>> # t/matrixops.t line 118 is: ok($esum == 61.308);
I'm ignoring the problem as I continue my run-through of
the install process but it looks like the root cause of this
problem can not be solved with a simple fix to eigens.
The eigens() routine does not handle complex results
nor the general matrix case where left- and right-
eigenvectors may be needed (returning NaNs is not
sufficient).
I suggest the following triage:
1. Fix CVS eigens() and documentation to handle symmetric
matrix inputs only (as in 2.4.2).
2. Remove broken support for non-symmetric matrix eigen
calculations.
3. Put general matrix eigen- support on the list for the
next major release.
One way forward would be to use the corresponding functions
from the CLAPACK library with possible addition of ATLAS
support to PDL as in the following links:
http://www.netlib.org/lapack/lug/lapack_lug.html
http://math-atlas.sourceforge.net/faq.html
--Chris
_______________________________________________
PDL-porters mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/pdl-porters
--
Bill Coffman
~~~_\\\|/
~~~-(@ @)
=oOO=(_)==OOo===_______________________________________________ Perldl mailing list |
