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

Johannes Loehnert a.u.r.e.l.i.a.n at gmx.net
Mon Jun 12 10:01:31 CDT 2006

```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]

```