Numeric Python, variable types and my sanity

[leeg]

I'm writing a simulation of a problem in quantum mechanics, which boils down
to populating a huge matrix, then finding its eigenvalues and eigenvectors
using MLab.eig() (or LinearAlgebra.eigenvectors()). One of the key facts
about matrices related to quantum operators is that they are Hermitian, and
therefore have definitely real eigenvalues.

In fact, the matrices I create are real and symmetric; in my code I raise an
exception should this not be so and this never happens unless I
deliberately break something to test the exception handler :-). Now a real
symmetric matrix *really* ought to have both real eigenvalues and real
eigenvectors, but there are certain values of the parameters in my program
(unfortunately, many of which I want to investigate) in which this is not
the case.

I know exactly what's happening here; precision error. An eigenvalue that
*ought* to be zero is coming out as something along the lines of
xe-15+ye-15j. The problem I have is that such a squiffy eigenvalue has a
near-arbitrary argument, and this *really* screws up the associated
eigenvector, so much so that it's useless. So far I've just been throwing
out any eigenvectors associated with complex artefacts, and converting the
remainder to pure reals on the fly. However I'm conscious that I could be
throwing away signal here.

Is there a way to force MLab to deal solely in real numbers? I'd like to
get it to not even try to compute the imaginary part of any results, so
that these artefacts do not occur and the eigenvalues retrieved by the
program are at least usable, even if they still have the precision error in
the real part.

Cheers,

leeg.
--
I am leeg, for we are many.
http://users.ox.ac.uk/~wadh1342

[David M. Cooke]

At some point, leeg <le**@teaching.physics.ox.ac.uk.valid> wrote:

In fact, the matrices I create are real and symmetric; in my code I raise an
exception should this not be so and this never happens unless I
deliberately break something to test the exception handler :-). Now a real
symmetric matrix *really* ought to have both real eigenvalues and real
eigenvectors, but there are certain values of the parameters in my program
(unfortunately, many of which I want to investigate) in which this is not
the case.

I know exactly what's happening here; precision error. An eigenvalue that
*ought* to be zero is coming out as something along the lines of
xe-15+ye-15j. The problem I have is that such a squiffy eigenvalue has a
near-arbitrary argument, and this *really* screws up the associated
eigenvector, so much so that it's useless. So far I've just been throwing
out any eigenvectors associated with complex artefacts, and converting the
remainder to pure reals on the fly. However I'm conscious that I could be
throwing away signal here.

I have a similiar use case, and ended up wrapping the LAPACK routine
DSYEVD myself (which does real, symmetric matrices), and using that.
Much better to use this than trying to coerce complex to real.

With f2py that should be pretty painless to do (I handwrote a C wrapper).

--
|>|\/|<
/--------------------------------------------------------------------------\
|David M. Cooke
|cookedm(at)physics(dot)mcmaster(dot)ca

[leeg]

David M. Cooke wrote:

I have a similiar use case, and ended up wrapping the LAPACK routine
DSYEVD myself (which does real, symmetric matrices), and using that.
Much better to use this than trying to coerce complex to real.

With f2py that should be pretty painless to do (I handwrote a C wrapper).

I don't have too much Fxx knowledge, and it's mainly F90. But I'll look
into that method. Thanks :-)
--
I am leeg, for we are many.
http://users.ox.ac.uk/~wadh1342