fix lapack's EIG on real matrices

Fixes quil-lang#177.

Thanks to @jcguu95 for the investigation.

Author: stylewarning

examples.lisp

@@ -140,15 +140,26 @@
     (format t "e^iH~%~a~%" expih)
     (format t "det(e^iH) = ~D~%" d)))
 
+(defun ensure-complex (x)
+  (etypecase x
+    (complex x)
+    (real (complex x))))
+
 (defun eig-printing (m)
-  (multiple-value-bind (vals vects)
-      (magicl:eig m)
-    (let ((val-diag (funcall #'magicl:from-diag vals)))
-      (format t "M~%~a~%" m)
-      (format t "Eigenvalues LAMBDA~%~a~%" val-diag)
-      (format t "Eigenvectors V~%~a~%" vects)
-      (format t "M*V~%~a~%" (magicl:@ m vects))
-      (format t "V*LAMBDA~%~a~%" (magicl:@ vects val-diag)))))
+  (let ((m-complex (magicl:.complex m (magicl:zeros (magicl:shape m)
+                                                    :type (magicl:element-type m)))))
+    (multiple-value-bind (vals vects)
+        (magicl:eig m)
+      (let ((val-diag (magicl:from-diag (mapcar #'ensure-complex vals))))
+        (format t "M~%~a~%" m)
+        (format t "Eigenvalues LAMBDA~%~a~%" val-diag)
+        (format t "Eigenvectors V~%~a~%" vects)
+        (format t "M*V~%~a~%" (magicl:@ m-complex vects))
+        (format t "V*LAMBDA~%~a~%" (magicl:@ vects val-diag))
+        (format t "M*V = V*LAMBDA (within 10^-10)? ~:[no~;yes~]~%"
+                (magicl:= (magicl:@ m-complex vects)
+                          (magicl:@ vects val-diag)
+                          1d-10))))))
 
 (defun eig-example ()
   (let ((m (magicl:from-list (list -2 -1 1 -2 1 1 -9 -3 4) '(3 3) :type 'double-float)))

src/extensions/lapack/lapack-templates.lisp

@@ -241,8 +241,6 @@
                    (from-array smat (list u-cols vt-rows) :input-layout :column-major)
                    (from-array vt (list vt-rows cols) :input-layout :column-major)))))))
 
-;; TODO: This returns only the real parts when with non-complex
-;; numbers. Should do something different?
 (defun generate-lapack-eig-for-type (class type eig-function &optional real-type)
   ` (defmethod lapack-eig ((m ,class))
       (policy-cond:with-expectations (> speed safety)
@@ -273,7 +271,55 @@
               ;; run it again with optimal workspace size
               (,eig-function jobvl jobvr rows a rows ,@(if real-type `(w) `(wr wi))
                              vl 1 vr rows work lwork ,@(when real-type `(rwork)) info)
-              (values (coerce ,@(if real-type `(w) `(wr)) 'list) (from-array vr (list rows cols) :input-layout :column-major))))))))
+              ,(if real-type
+                   `(values (coerce w 'list)
+                            (from-array vr (list rows cols) :input-layout :column-major))
+                   `(values (cl:map 'list (lambda (a b)
+                                            (if (zerop b)
+                                                a
+                                                (complex a b)))
+                                    wr wi)
+                            (let* ((evecs (magicl:zeros (list rows cols) :type '(complex ,type)))
+                                   (storage (magicl::storage evecs)))
+                              ;; square matrix
+                              (loop :with col-lapack := 0
+                                    :with col-result := 0
+                                    :with skip := nil
+                                    :for zr :across wr
+                                    :for zi :across wi
+                                    :do (cond
+                                          ;; real eigenvalue
+                                          ((zerop zi)
+                                           (when skip
+                                             (warn "SKIP is T when we reached a real eigenvalue."))
+                                           (dotimes (r rows)
+                                             ;; column-major
+                                             (setf (aref storage (+ r (* col-result rows)))
+                                                   (complex (aref vr (+ r (* col-lapack rows))))))
+                                           (incf col-result)
+                                           (incf col-lapack))
+                                          ;; complex eigenvalue with conjugate
+                                          (skip
+                                           (unless (cl:= skip (- zi))
+                                             (error "Reached a non-conjugate eigenvalue"))
+                                           (setf skip nil))
+                                          ;; New complex eigenvalue
+                                          (t
+                                           ;; expect a conjugate in the next iteration
+                                           (setf skip zi)
+                                           (dotimes (r rows)
+                                             ;; column-major
+                                             (setf (aref storage (+ r (* col-result rows)))
+                                                   (complex
+                                                    (aref vr (+ r (* col-lapack rows)))
+                                                    (aref vr (+ r (* (1+ col-lapack) rows)))))
+                                             (setf (aref storage (+ r (* (1+ col-result) rows)))
+                                                   (complex
+                                                    (aref vr (+ r (* col-lapack rows)))
+                                                    (- (aref vr (+ r (* (1+ col-lapack) rows)))))))
+                                           (incf col-result 2)
+                                           (incf col-lapack 2)))
+                                    :finally (return evecs)))))))))))
 
 (defun generate-lapack-hermitian-eig-for-type (class type eig-function real-type)
   `(defmethod lapack-hermitian-eig ((m ,class))