From 6fbb3c1a0c500c021844cf0a50005ac8e01ab65c Mon Sep 17 00:00:00 2001
From: Robert Smith <rssmith@hrl.com>
Date: Tue, 6 Sep 2022 19:29:35 -0700
Subject: [PATCH] wip

---
 src/compilers/approx.lisp  | 11 +++++++++
 src/matrix-operations.lisp | 50 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 61 insertions(+)

diff --git a/src/compilers/approx.lisp b/src/compilers/approx.lisp
index ea9a9470a..14d52e64f 100644
--- a/src/compilers/approx.lisp
+++ b/src/compilers/approx.lisp
@@ -234,6 +234,17 @@ Three self-explanatory restarts are offered: TRY-AGAIN, and GIVE-UP-COMPILATION.
       (give-up-compilation))))
 
 (defun orthogonal-decomposition (m)
+  (multiple-value-bind (diag-re diag-im left right)
+      (magicl:qz (.realpart m) (.imagpart m))
+    (declare (ignore right))
+    (let ((diag (.complex diag-re diag-im)))
+      (when (minusp (magicl:det left))
+        (dotimes (row (magicl:nrows left))
+          (setf (magicl:tref left row 0)
+                (- (magicl:tref left row 0)))))
+      (values left (magicl:@ diag diag) (magicl:transpose left)))))
+
+(defun old-orthogonal-decomposition (m)
   "Extracts from M a decomposition of E^* M E into A * D * B, where A and B are orthogonal and D is diagonal.  Returns the results as the VALUES triple (VALUES A D B)."
   (let* ((m (magicl:scale m (expt (magicl:det m) -1/4)))
          (a (diagonalizer-in-e-basis m))
diff --git a/src/matrix-operations.lisp b/src/matrix-operations.lisp
index 10a21c481..9426402aa 100644
--- a/src/matrix-operations.lisp
+++ b/src/matrix-operations.lisp
@@ -360,3 +360,53 @@ M1 and M2 are n x n."
              (magicl:nrows m2) (magicl:ncols m2)))
   (sqrt (- 1 (/ (abs (magicl:trace (magicl:@ m1 (magicl:dagger m2))))
                 (magicl:nrows m1)))))
+
+(defgeneric .realpart (m)
+  (:documentation "Compute the real part of all entries of the matrix M.")
+  (:method ((m magicl:matrix/single-float))
+    m)
+  (:method ((m magicl:matrix/double-float))
+    m)
+  (:method ((m magicl:matrix/complex-single-float))
+    (let ((re-m (magicl:zeros (magicl:shape m) :type 'single-float)))
+      (magicl::map-to #'realpart m re-m)
+      re-m))
+  (:method ((m magicl:matrix/complex-double-float))
+    (let ((re-m (magicl:zeros (magicl:shape m) :type 'double-float)))
+      (magicl::map-to #'realpart m re-m)
+      re-m)))
+
+(defgeneric .imagpart (m)
+  (:documentation "Compute the imaginary part of all entries of the matrix M.")
+  (:method ((m magicl:matrix/single-float))
+    (magicl:zeros (magicl:shape m) :type 'single-float))
+  (:method ((m magicl:matrix/double-float))
+    (magicl:zeros (magicl:shape m) :type 'double-float))
+  (:method ((m magicl:matrix/complex-single-float))
+    (let ((im-m (magicl:zeros (magicl:shape m) :type 'single-float)))
+      (magicl::map-to #'imagpart m im-m)
+      im-m))
+  (:method ((m magicl:matrix/complex-double-float))
+    (let ((im-m (magicl:zeros (magicl:shape m) :type 'double-float)))
+      (magicl::map-to #'imagpart m im-m)
+      im-m)))
+
+(defgeneric .complex (re im)
+  (:method ((re magicl:matrix/single-float) (im magicl:matrix/single-float))
+    (assert (equal (magicl:shape re) (magicl:shape im)))
+    (let ((z (magicl:zeros (magicl:shape re)
+                           :type '(complex single-float))))
+      (magicl::into! (lambda (i j)
+                       (complex (magicl:tref re i j)
+                                (magicl:tref im i j)))
+                     z)
+      z))
+  (:method ((re magicl:matrix/double-float) (im magicl:matrix/double-float))
+    (assert (equal (magicl:shape re) (magicl:shape im)))
+    (let ((z (magicl:zeros (magicl:shape re)
+                           :type '(complex double-float))))
+      (magicl::into! (lambda (i j)
+                       (complex (magicl:tref re i j)
+                                (magicl:tref im i j)))
+                     z)
+      z)))