[SciPy-dev] Handling of float96 and complex192 in linalg

[Johannes Loehnert]

Hi,

I was looking through the code of the linalg module last week. I noticed that 
several of the functions use ``get_lapack_funcs`` from 
``Lib/linalg/lapack.py`` in order to find the appropriate lapack routine for 
a given data type.

However, ``get_lapack_funcs`` fails on arrays with dtype.char=='g' (float96) 
and dtype.char=='G' (complex 2*96). In both cases the double precision 
routine is given back. (I know there are no float96/complex192 routines in 
lapack).

While this somewhat acceptable for float96 (precision suffers), it is not very 
nice if complex192 is silently cast to float64 by simply throwing the 
imaginary part away. (See below for a "fatal" example.)

The "lazy" solution would be to let get_lapack_funcs raise an error for 'g' 
and 'G' typechar. 

Alternatively, complex192 could be treated by double precision complex 
routines ('z' prefix). For this approach, some kind of warning would be nice 
if float96 or complex192 arrays are cast to a lower precision.


Johannes Loehnert

Example:
In [1]: from numpy import *

In [2]: re = rand(10,10)

In [3]: im = rand(10,10)

In [4]: x = (re + 1j*im).astype('G')   # create random complex192 matrix

In [5]: from scipy.linalg import eig

In [6]: e1, dummy = eig(x)  # x is cast to float64, yielding wrong eigenvalues

In [7]: e2, dummy = eig(x.astype('D'))  # cast to complex128 -> correct values

In [8]: print e1
[ 4.99532769+0.j          1.33254924+0.j         -0.08283794+0.76930589j
 -0.08283794-0.76930589j -0.47488469+0.24342446j -0.47488469-0.24342446j
  0.33556027+0.15301716j  0.33556027-0.15301716j -0.02890708+0.31652241j
 -0.02890708-0.31652241j]

In [9]: print e2
[ 5.02636564+5.37963067j  1.50114472-0.01585185j  0.60033384-1.06642284j
 -0.46327973-0.84463702j -0.62357122-0.20919303j -0.65416101+0.66355244j
 -0.47412464+0.7726915j   0.37856415+0.8177398j   0.15827675-0.07332642j
  0.37618956+0.34544957j]