GSL does not seem to support non-symmetric or
non-Hermitian matrices at this time. If we want general
eigensystem solvers, LAPACK and its family might still be
the way to go. We can use GSL for the symmetric cases.
On a related note, it would be very useful if we could get
some sort of feedback on what functions and functionality
from PDL are being used.
Perhaps we could add some instrumentation to PDL to
collect various stats (module usage, function usage,...)
and periodically give a chance to email the results back
to our development team along with any suggestions.
--Chris
On July 09, 2006 1:37:26 PM, Craig DeForest wrote:
>
> 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 wrote:
> >
> > ... 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
> > ..._______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
