EIG on real matrices is buggy and wrong

[stylewarning]

The EIG call to LAPACK is buggy and wrong. It lops off imaginary parts and just returns real components. The vectors are also all real.

CL-USER> *m
#<MAGICL:MATRIX/DOUBLE-FLOAT (3x3):
  -0.683    -2.002    -2.352
   0.810     0.292    -0.054
   1.921     0.673     2.364>

CL-USER> (magicl::eig *m) ; lapack
(0.5488087694607198d0 0.5488087694607198d0 0.8753843383425303d0)
#<MAGICL:MATRIX/DOUBLE-FLOAT (3x3):
  -0.758     0.000    -0.366
  -0.052     0.321    -0.575
   0.441     0.355     0.732>

CL-USER> (magicl::eig-lisp *m) ; lispified lapack
(#C(0.5488087694607201d0 1.9472471755824976d0)
 #C(0.5488087694607201d0 -1.9472471755824976d0) 0.8753843383425299d0)
#<MAGICL:MATRIX/COMPLEX-DOUBLE-FLOAT (3x3):
   0.632 + 0.074j     0.632 - 0.074j     0.374 + 0.000j
   0.000 + 0.531j     0.000 - 0.531j     0.417 + 0.000j
   0.404 - 0.679j     0.404 + 0.679j     0.614 + 0.000j>

I won't claim with 100% certainty eig-lisp is correct, but you can manually verify with Ax = lx that it is correct enough.

[jcguu95]

The eigenvalues you posted are correct (lisp backend), according to Wolframalpha. I can reproduce this bug with

(defparameter *m (from-list '(-0.683 -2.002 -2.352 0.810 0.292 -0.054 1.921 0.673 2.364) '(3 3)))
(defparameter *m2 (from-list '(#c(-0.683 0) -2.002 -2.352 0.810 0.292 -0.054 1.921 0.673 2.364) '(3 3)))
(funcall (magicl.backends:backend-implementation 'eig :lapack) *m)

The bug does not happen to *m2.

The bug is in the last form in the definition of #'generate-lapack-eig-for-type:

(values (coerce ,@(if real-type `(w) `(wr)) 'list) (from-array vr (list rows cols) :input-layout :column-major))))))))

The (wr) should have been something like ((.complex-lisp wr wi)). However, when I made such change

              (values (coerce ,@(if real-type
                                    `(w)
                                    `((let ((magicl.backends::*backends* `(:lisp)))
                                        (magicl::.complex-lisp wr wi))))
                              'list)
                      (from-array vr (list rows cols) :input-layout :column-major))

it errors out because #'complex-lispdoes not care my dynamic binding for *backends* and keep tries to use lapack for .complex. The definition of these functions are all done by a huge macro.. and I didn't know how to make it work. One quick around is to reimplement #'.complex-lisp there without depending on *backends*, but that is not satisfactory.

Please let me know about your thought.

[stylewarning]

@jcguu95 Thanks a lot for looking into this! I'll try to look soon. What you say sounds promising.

[jcguu95]

Thanks for building this too. By the way I forgot to say, that wr denotes the real part and wi denotes the imaginary part. Changing (wr) to ((wr wi)) does return the imaginary part correctly - I've tested this. So the only thing left is indeed to glue them together with #'.complex.