add MAGICL:CSD-BLOCKS, define MAGICL:CSD on it, and other lisp-native routines

.gitignore

@@ -1,2 +1,3 @@
 **/*.dylib
 **/*.so
+.DS_Store

magicl.asd

@@ -15,8 +15,7 @@
   :depends-on (#:magicl/core
                #:magicl/ext
                #:magicl/ext-blas
-               #:magicl/ext-lapack
-               #:magicl/ext-expokit))
+               #:magicl/ext-lapack))
 
 (asdf:defsystem #:magicl/core
   :license "BSD 3-Clause (See LICENSE.txt)"
@@ -59,7 +58,21 @@
                  (:file "types/specialized-vector")
                  (:file "constructors")
                  (:file "specialize-constructor")
-                 (:file "polynomial-solver")))
+                 (:file "polynomial-solver")
+                 (:module "matrix-functions"
+                  :serial t
+                  :components ((:file "row-column-matrices")
+                               (:file "mult-definition")
+                               (:file "mult-methods")
+                               (:file "givens")
+                               (:file "householder")
+                               (:file "qr")
+                               (:file "lu")
+                               (:file "tridiagonal-form")
+                               (:file "hermitian-eig")
+                               (:file "svd")
+                               (:file "random-hermitian")
+                               (:file "csd")))))
    (:file "magicl")))
 
 ;;; Extension common code

src/backend-function.lisp

@@ -61,30 +61,57 @@
 ;;; now the meat
 
 (define-condition no-applicable-implementation (simple-error)
-  ()
+  ((function-name :initarg :function-name :reader no-applicable-implementation-function-name)
+   (failed-callees :initarg :failed-callees :reader no-applicable-implementation-failed-callees
+                   :documentation "A list of NO-APPLICABLE-IMPLEMENTATION conditions that were a result of callee failures."))
   (:documentation "The error to signal when there's no applicable implementation, whether from a backend function or a CLOS generic function."))
 
-(defun no-applicable-implementation (name)
+(defun describe-callees (s conditions &optional (indent 1))
+  ;; used to format NO-APPLICABLE-IMPLEMENTATION errors
+  (etypecase conditions
+    ((member :|<unavailable>|)
+     (format s "~vT* <unavailable>~%" (* 4 indent)))
+    (no-applicable-implementation
+     (format s "~vT* ~S~%" (* 4 indent)
+             (no-applicable-implementation-function-name conditions))
+     (describe-callees s (no-applicable-implementation-failed-callees conditions) (1+ indent)))
+    (list
+     (dolist (c conditions)
+       (describe-callees s c indent)))))
+
+(defun nai-fmt (s arguments &optional colon-modifier at-modifier)
+  ;; used to format NO-APPLICABLE-IMPLEMENTATION errors
+  (declare (ignore colon-modifier at-modifier))
+  (destructuring-bind (name callees) arguments
+    (let ((*package* (find-package "KEYWORD")))
+      (format s "No applicable implementations found for the (supposed) ~
+                 backend function ~S.~2%" name)
+      (describe-backend-function name s)
+      (terpri s)
+      (format s "There was no implementation of ~S that was applicable to the given ~
+                 arguments. This may be because:~2%" name)
+      (format s "    * An implementation exists, but a function this implementation ~
+                       itself uses doesn't have an implementation. (See below.)~2%")
+      (format s "    * No implementation exists in a currently active backend,~2%")
+      (format s "    * A generic function implements the backend, but the generic ~
+                       function does not specialize on the argument, or~2%")
+      (format s "    * An implementation exists, but signaled that it wasn't ~
+                       applicable for the given arguments.~2%")
+      (unless (null callees)
+        (format s "The callees that failed to find an implementation were:~2%")
+        (describe-callees s callees))
+      (terpri s)
+      (format s "Ensure the proper backends are activated, and the implementations ~
+                 in those backends can handle the argument types appropriately. If ~
+                 an implementation does not exist, consider writing one!"))))
+
+(defun no-applicable-implementation (name &key (failed-callees ':|<unavailable>|))
   "Call this function to signal an error indicating the caller (or any callers above it) are not applicable to the current backend function being invoked."
   (error 'no-applicable-implementation
-         :format-control
-         (let ((*package* (find-package "KEYWORD")))
-           (with-output-to-string (s)
-             (format s "No applicable implementations found for the supposed ~
-                        backend function ~S.~2%" name)
-             (describe-backend-function name s)
-             (terpri s)
-             (format s "There was no implementation of ~S that was applicable to the given ~
-                        arguments. This may be because:~2%" name)
-             (format s "    * No implementation exists in a currently active backend,~2%")
-             (format s "    * A generic function implements the backend, but the generic ~
-                              function does not specialize on the argument, or~2%")
-             (format s "    * An implementation exists, but signaled that it wasn't ~
-                              applicable for the given arguments.~2%")
-             (format s "Ensure the proper backends are activated, and the implementations ~
-                        in those backends can handle the argument types appropriately. If ~
-                        an implementation does not exist, consider writing one!")))
-         :format-arguments nil))
+         :function-name name
+         :failed-callees failed-callees
+         :format-control "~/magicl.backends::nai-fmt/"
+         :format-arguments (list (list name failed-callees))))
 
 (defmacro define-compatible-no-applicable-method-behavior (&rest generic-function-names)
   "Ensure the (unquoted symbol) generic function names GENERIC-FUNCTION-NAMES will behave correctly when used as implementations for backend functions.
@@ -99,7 +126,7 @@ Without using this, a backend function may error if no method is found."
          :collect
          `(defmethod cl:no-applicable-method ((gf (eql #',name)) &rest args)
             (declare (ignore args))
-            (no-applicable-implementation ',name)))))
+            (no-applicable-implementation ',name :failed-callees nil)))))
 
 ;;; Backend Names
 
@@ -115,7 +142,7 @@ Without using this, a backend function may error if no method is found."
     "A list in priority order the backends that MAGICL should use for functionality.
 
 It is preferable to use WITH-BACKENDS instead of this.")
-  
+
   (defun backend-name-p (x)
     (and (symbolp x)
          (find x *known-backends*)
@@ -224,16 +251,16 @@ Backend implementations are defined with DEFINE-BACKEND-IMPLEMENTATION."
        (initialize-function-for-implementations ',name)
        (defun ,name ,lambda-list
          ,@(and doc (list doc))
-         (dolist (,backend *backend*)
-           (handler-case
-               (let ((,func (backend-implementation ',name ,backend)))
-                 (unless (null ,func)
-                   (return-from ,name
-                     ,(interface:calling-form func lambda-list))))
-             (no-applicable-implementation (c)
-               (declare (ignore c))
-               nil)))
-         (no-applicable-implementation ',name)))))
+         (let ((callees nil))
+           (dolist (,backend *backend*)
+             (handler-case
+                 (let ((,func (backend-implementation ',name ,backend)))
+                   (unless (null ,func)
+                     (return-from ,name
+                       ,(interface:calling-form func lambda-list))))
+               (no-applicable-implementation (c)
+                 (push c callees))))
+           (no-applicable-implementation ',name :failed-callees callees))))))
 
 (defmacro define-backend-implementation (name backend funcallable-expression)
   "Define the implementation of the backend function named NAME, for the backend BACKEND, as the function FUNCALLABLE-EXPRESSION (evaluated).
@@ -289,7 +316,7 @@ NOTE: If your implementation is a generic function, then the generic function's
            (terpri stream)
            (cond
              ((null inactive)
-              (format stream "There are no inactive backends.~%"))
+              (format stream "There are no loaded backends that are inactive.~%"))
              (t
               (format stream "Inactive Backends:~%")
               (dolist (b inactive)

src/extensions/lapack/lapack-csd.lisp

@@ -89,7 +89,7 @@
 
       (values mu1 mu2 mv1h mv2h (list (atan s c))))))
 
-(defmethod lapack-csd ((x matrix/complex-double-float) p q)
+(defmethod csd-blocks-extension ((x matrix/complex-double-float) p q)
   (let* ((m (nrows x))
          (layout (layout x))
          (xcopy (deep-copy-tensor x)))
@@ -173,43 +173,3 @@
                   (from-array v2t (list (- m q) (- m q)) :input-layout layout)
                   (coerce theta 'list)))))))
 
-(defmethod csd-extension ((matrix matrix/complex-double-float) p q)
-  (labels ((csd-from-blocks (u1 u2 v1t v2t theta)
-             "Calculates the matrices U, SIGMA, and VT of the CSD of a matrix from its intermediate representation, as calculated from LAPACK-CSD."
-             (let ((p (nrows u1))
-                   (q (nrows v1t))
-                   (m (+ (nrows u1) (nrows u2)))
-                   (r (length theta)))
-               (let ((u (direct-sum u1 u2))
-                     (sigma (const 0 (list m m) :type (element-type matrix)))
-                     (vt (direct-sum v1t v2t)))
-                 (let ((diag11 (min p q))
-                       (diag12 (min p (- m q)))
-                       (diag21 (min (- m p) q))
-                       (diag22 (min (- m p) (- m q))))
-                   (let ((iden11 (- diag11 r))
-                         (iden12 (- diag12 r))
-                         (iden21 (- diag21 r))
-                         (iden22 (- diag22 r)))
-                     ;; Construct sigma from theta
-                     (loop :for i :from 0 :to (1- iden11)
-                           do (setf (tref sigma i i) 1))
-                     (loop :for i :from iden11 :to (1- diag11)
-                           do (setf (tref sigma i i) (cos (nth (- i iden11) theta))))
-                     (loop :for i :from 0 :to (1- iden12)
-                           do (setf (tref sigma (- p 1 i) (- m 1 i)) -1))
-                     (loop :for i :from iden12 :to (1- diag12)
-                           do (setf (tref sigma (- p 1 i) (- m 1 i))
-                                    (- (sin (nth (- r 1 (- i iden12)) theta)))))
-                     (loop :for i :from 0 :to (1- iden21)
-                           do (setf (tref sigma (- m 1 i) (- q 1 i)) 1))
-                     (loop :for i :from iden21 :to (1- diag21)
-                           do (setf (tref sigma (- m 1 i) (- q 1 i))
-                                    (sin (nth (- r 1 (- i iden21)) theta))))
-                     (loop :for i :from 0 :to (1- iden22)
-                           do (setf (tref sigma (+ p i) (+ q i)) 1))
-                     (loop :for i :from iden22 :to (1- diag22)
-                           do (setf (tref sigma (+ p i) (+ q i)) (cos (nth (- i iden22) theta))))))
-                 (values u sigma vt)))))
-    (multiple-value-bind (u1 u2 v1t v2t theta) (lapack-csd matrix p q)
-      (csd-from-blocks u1 u2 v1t v2t theta))))

src/extensions/lapack/lapack-generics.lisp

@@ -19,9 +19,7 @@
 
 (magicl:define-extensible-function (magicl:inv lapack-inv :lapack) (matrix))
 
-(magicl:define-extensible-function (magicl:csd csd-extension :lapack) (matrix p q))
-
-(defgeneric lapack-csd (matrix p q))
+(magicl:define-extensible-function (magicl:csd-blocks csd-blocks-extension :lapack) (matrix p q))
 
 (magicl:define-extensible-function (magicl:svd lapack-svd :lapack) (matrix &key reduced))
 

src/extensions/lapack/lapack-templates.lisp

@@ -15,24 +15,24 @@
        (policy-cond:with-expectations (> speed safety)
            ((type (member nil :n :t :c) transa)
             (type (member nil :n :t :c) transb))
-         (let* ((m (if (eq :n transa) (nrows a) (ncols a)))
-                (k (if (eq :n transa) (ncols a) (nrows a)))
-                (n (if (eq :n transb) (ncols b) (nrows b)))
-                (brows (if (eq :n transb) (nrows b) (ncols b))))
+         (let* ((m (if (eq ':n transa) (nrows a) (ncols a)))
+                (k (if (eq ':n transa) (ncols a) (nrows a)))
+                (n (if (eq ':n transb) (ncols b) (nrows b)))
+                (brows (if (eq ':n transb) (nrows b) (ncols b))))
            (policy-cond:with-expectations (> speed safety)
                ((assertion (cl:= k brows))
                 (assertion (or (not target) (equal (shape target) (list m n)))))
              (let ((ta
-                     (if (eql :row-major (layout a))
+                     (if (eq ':row-major (layout a))
                          (case transa
-                           (:n :t)
-                           (:t :n)
+                           (:n ':t)
+                           (:t ':n)
                            (:c (error "Specifying TRANSA to be :C is not allowed if A is ROW-MAJOR")))
                          transa))
-                   (tb (if (eql :row-major (layout b))
+                   (tb (if (eq ':row-major (layout b))
                            (case transb
-                             (:n :t)
-                             (:t :n)
+                             (:n ':t)
+                             (:t ':n)
                              (:c (error "Specifying TRANSB to be :C is not allowed if B is ROW-MAJOR")))
                            transb))
                    (target (or target
@@ -53,9 +53,9 @@
                 k
                 alpha
                 (magicl::storage a)
-                (if (eql :n ta) m k)
+                (if (eq ':n ta) m k)
                 (magicl::storage b)
-                (if (eql :n tb) k n)
+                (if (eq ':n tb) k n)
                 beta
                 (magicl::storage target)
                 m)
@@ -64,16 +64,16 @@
        (policy-cond:with-expectations (> speed safety)
            ((type (member nil :n :t :c) transa)
             (assertion (null transb)))
-         (let* ((m-op (if (eq :n transa) (nrows a) (ncols a)))
-                (n-op (if (eq :n transa) (ncols a) (nrows a))))
+         (let* ((m-op (if (eq ':n transa) (nrows a) (ncols a)))
+                (n-op (if (eq ':n transa) (ncols a) (nrows a))))
            (policy-cond:with-expectations (> speed safety)
                ((assertion (cl:= n-op (size x)))
                 (assertion (or (not target) (equal (shape target) (list m-op)))))
              (let ((ta
-                     (if (eql :row-major (layout a))
+                     (if (eq ':row-major (layout a))
                          (case transa
-                           (:n :t)
-                           (:t :n)
+                           (:n ':t)
+                           (:t ':n)
                            (:c (error "Specifying TRANS to be :C is not allowed if A is ROW-MAJOR")))
                          transa))
                    (target (or target
@@ -85,8 +85,8 @@
                   (:t "T")
                   (:c "C")
                   (:n "N"))
-                (if (eql :column-major (layout a)) (nrows a) (ncols a))
-                (if (eql :column-major (layout a)) (ncols a) (nrows a))
+                (if (eq ':column-major (layout a)) (nrows a) (ncols a))
+                (if (eq ':column-major (layout a)) (ncols a) (nrows a))
                 alpha
                 (magicl::storage a)
                 (if (eql :column-major (layout a)) (nrows a) (ncols a))

src/extensions/lapack/package.lisp

@@ -26,5 +26,4 @@
    #:lapack-ql-q
    #:lapack-qr-q
    #:lapack-rq-q
-   #:lapack-lq-q
-   #:lapack-csd))
+   #:lapack-lq-q))

src/high-level/abstract-tensor.lisp

@@ -265,8 +265,8 @@ If TARGET is not specified then a new tensor is created with the same element ty
 
 (define-backend-implementation ./ :lisp
   (lambda (source1 source2 &optional target)
-    ;; This won't do the right thing if / is not closed in the set of
-    ;; choice, like integers.
+    ;; FIXME: This won't do the right thing if / is not closed in the
+    ;; set of choice, like integers.
     (binary-operator #'/ source1 source2 target)))
 
 (define-backend-function .^ (source1 source2 &optional target)
@@ -275,8 +275,6 @@ If TARGET is not specified then a new tensor is created with the same element ty
 
 (define-backend-implementation .^ :lisp
   (lambda (source1 source2 &optional target)
-    ;; This won't do the right thing if / is not closed in the set of
-    ;; choice, like integers.
     (binary-operator #'expt source1 source2 target)))
 
 (define-extensible-function (.max max-lisp) (source1 source2 &optional target)
@@ -320,7 +318,6 @@ If TARGET is not specified then a new tensor is created with the same element ty
   (:method ((source abstract-tensor) &optional target)
     (unary-operator #'log source target)))
 
-
 (define-extensible-function (= =-lisp) (source1 source2 &optional epsilon)
   (:documentation "Check the equality of tensors with an optional EPSILON")
   (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional epsilon)