diff --git a/README.md b/README.md index ce362881..d6a50674 100755 --- a/README.md +++ b/README.md @@ -39,6 +39,10 @@ You can run the MAGICL tests from your Lisp REPL with: (asdf:test-system :magicl) ``` +## High-level Interface + +See [high-level doc](doc/high-level.md). + ## Showing Available Functions Some distributions of a library don't actually provide all of the functions of the reference BLAS and LAPACK. One can look at a summary of available and unavailable functions with the function `magicl:print-availability-report`. By default, it will show all functions and their availability. There are three arguments to fine-tune this behavior: diff --git a/VERSION.txt b/VERSION.txt index 55e30fcc..482f2571 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -"0.6.5" +"0.7.0" diff --git a/doc/high-level.md b/doc/high-level.md new file mode 100644 index 00000000..ea6c1ca4 --- /dev/null +++ b/doc/high-level.md @@ -0,0 +1,83 @@ +# High Level Bindings + +The purpose of the high-level magicl bindings is to allow for MATLAB-like multidimensional arrays in lisp. + +## Constructors + +The construction of tensors can be done with any of the given constructors. The constructors take a shape and arguments for method of construction. + +Tensors are specialized on both the shape and the element type. The class of a tensor will be of the form `$CLASS/$TYPE` (e.g. `MATRIX/DOUBLE-FLOAT`). + +| Number of dimensions | Tensor Class | +|----------------------|-----------------| +| 1 | `VECTOR` | +| 2 | `MATRIX` | +| * | `TENSOR` | + +| Element Type | Class Suffix | +|--------------------------|------------------------| +| `SINGLE-FLOAT` | `SINGLE-FLOAT` | +| `DOUBLE-FLOAT` | `DOUBLE-FLOAT` | +| `(COMPLEX SINGLE-FLOAT)` | `COMPLEX-SINGLE-FLOAT` | +| `(COMPLEX DOUBLE-FLOAT)` | `COMPLEX-DOUBLE-FLOAT` | +| `(SIGNED-BYTE 32)` | `INT32` | + +The type of the elements of the tensor can be specified with the `:type` keyword, or the constructor will attempt to find an appropriate type from the given arguments. The default element type for a tensor is `DOUBLE-FLOAT`. + +The layout of the tensor (column-major or row-major) can be specified with the `:layout` keyword. This affects the underlying storage of the tensor and will affect how it carries out operations with LAPACK. + + +## Other Library Equivalents + +This table was adapted largely from the [NumPy Equivalents Table](https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html#linear-algebra-equivalents). + +### Basic Accessors + +| MAGICL | MATLAB | NumPy | Description | +|----------------|------------|-------------------------|---------------------------------------------------------------| +| `(order a)` | `ndims(a)` | `ndim(a)` or `a.ndim` | Get the number of dimensions of the array. | +| `(size a)` | `numel(a)` | `size(a)` or `a.size` | Get the number of elements of the array. | +| `(shape a)` | `size(a)` | `shape(a)` or `a.shape` | Get the shape of the array. | +| `(tref a 1 4)` | `a(2,5)` | `a[1, 4]` | Get the element in the second row, fifth column of the array. | + +### Constructors + +| MAGICL | MATLAB | NumPy | Description | +|-------------------------------------------------|--------------------|-----------------------------------|--------------------------------------------------------------------------------------| +| `(from-list '(1d0 2d0 3d0 4d0 5d0 6d0) '(2 3))` | `[ 1 2 3; 4 5 6 ]` | `array([[1.,2.,3.], [4.,5.,6.]])` | Create a 2x3 matrix from given elements. | +| `(zeros '(2 3 4))` or `(const 0d0 '(2 3 4))` | `zeros(2,3,4)` | `zeros((2,3,4))` | Create a 2x3x4 dimensional array of zeroes of double-float element type. | +| `(ones '(3 4)` or `(const 1d0 '(3 4))` | `ones(3,4)` | `ones((3,4))` | Create a 3x4 dimensional array of ones of double-float element type. | +| `(eye 1d0 '(3 3))` | `eye(3)` | `eye(3)` | Create a 3x3 identity array of double-float element type. | +| `(from-diag a)` | `diag(a)` | `diag(a)` | Create a square matrix from the diagonal entries in `a` with zeroes everywhere else. | +| `(rand '(3 4))` | `rand(3,4)` | `random.rand(3,4)` | Create a random 3x4 array. | + +### Basic Operations + +| MAGICL | MATLAB | NumPy | Description | +|------------|----------|------------------|-----------------------------| +| `(@ a b)` | `a * b` | `a @ b` | Matrix multiplication | +| `(.+ a b)` | `a + b` | `a + b` | Element-wise add | +| `(.- a b)` | `a - b` | `a - b` | Element-wise subtract | +| `(.* a b)` | `a .* b` | `a * b` | Element-wise multiply | +| `(./ a b)` | `a./b` | `a/b` | Element-wise divide | +| `(.^ a b)` | `a.^b` | `np.power(a,b)` | Element-wise exponentiation | + +### Linear Algebra + +| MAGICL | MATLAB | NumPy | Description | +|---------------------------------------------------------|-------------------|----------------------------------------------------------------|--------------------------------------------| +| `(det a)` | `det(a)` | `linalg.det(a)` | Determinant of matrix | +| `(trace a)` | `trace(a)` | `trace(a)` | Trace (sum of diagonal elements) of matrix | +| `(upper-triangular a)` | `triu(a)` | `triu(a)` | Upper triangular part of matrix | +| `(lower-triangular a)` | `tril(a)` | `tril(a)` | Lower triangular part of matrix | +| `(transpose a)` | `a.'` | `a.transpose()` or `a.T` | Transpose of matrix | +| `(conjugate-transpose a)` or `(dagger a)` | `a'` | `a.conj().transpose()` or `a.H` | Conjugate transpose of matrix | +| `(inv a)` | `inv(a)` | `linalg.inv(a)` | Inverse of matrix | +| `(svd a)` (Returns `(VALUES U SIGMA Vt)`) | `[U,S,V]=svd(a)` | `U, S, Vh = linalg.svd(a), V = Vh.T` | Singular value decomposition of matrix | +| `(eig a)` (Returns `(VALUES EIGENVALUES EIGENVECTORS)`) | `[V,D]=eig(a)` | `D,V = linalg.eig(a)` | Eigenvalues and eigenvectors of matrix | +| `(qr a)` (Returns `(VALUES Q R)`) | `[Q,R,P]=qr(a,0)` | `Q,R = scipy.linalg.qr(a)` | QR factorization of matrix | +| `(ql a)` (Returns `(VALUES Q L)`) | | | QL factorization of matrix | +| `(rq a)` (Returns `(VALUES R Q)`) | | `R,Q = scipy.linalg.rq(a)` | RQ factorization of matrix | +| `(lq a)` (Returns `(VALUES L Q)`) | | | LQ factorization of matrix | +| `(lu a)` (Returns `(VALUES LU IPIV)`) | `[L,U,P]=lu(a)` | `L,U = scipy.linalg.lu(a)` or `LU,P=scipy.linalg.lu_factor(a)` | LU decomposition of matrix | +| `(csd a)` (Returns `(VALUES U SIGMA VT)`) | | | Cosine-sine decomposition of matrix | diff --git a/examples.lisp b/examples.lisp index 054c6460..70290c08 100644 --- a/examples.lisp +++ b/examples.lisp @@ -1,5 +1,5 @@ (defpackage #:magicl-examples - (:use #:common-lisp #:magicl #:magicl-transcendental) + (:use #:common-lisp #:magicl-transcendental) #+package-local-nicknames (:local-nicknames (:blas :magicl.blas-cffi) (:lapack :magicl.lapack-cffi) @@ -16,238 +16,148 @@ (in-package #:magicl-examples) +(defconstant +double-comparison-threshold-loose+ (* 256 double-float-epsilon)) + (defun print-matrix (m) (princ m) (terpri) (terpri)) -(defun dot (u v) - (assert (= (length u) (length v)) (u v)) - (let ((n (length u))) - (with-blapack - (let ((cx (magicl::copy-matrix-storage u)) - (cy (magicl::copy-matrix-storage v))) - (format t "x: ~A~%y: ~A~%" cx cy) - (blas:%zdotc - n - cx - 1 - cy - 1))))) - (defun dot-example () - (let ((a (magicl::make-C-storage 4)) - (b (magicl::make-C-storage 4))) - (dotimes (i 4) - (setf (aref a i) #C(1.0f0 0.0f0) - (aref b i) #C(2.0f0 0.0f0))) - (format t "a^t = ~A~%b^t = ~A~%a^t b = ~A~%~%" - a b (blas:%cdotu 4 a 1 b 1)))) + (let ((a (magicl:const #C(1d0 0d0) '(4))) + (b (magicl:const #C(2d0 0d0) '(4)))) + (format t "a = ~A~%b = ~A~% = ~A~%~%" + a b (magicl:dot a b)))) (defun eigenvalue-example () - ;; Set the traps - (with-blapack - - ;; An eigenvalue example. Note that we have no matrix abstraction a - ;; this point. We pretend 4-vectors are 2-by-2 matrices. - - ;; BLAS/LAPACK expects column major order, we are creating the - ;; (matlab notation) matrix M = [1 2; 2 3]. - (let ((M (magicl::make-D-storage 4))) - (setf (aref M 0) 1.0d0 - (aref M 1) 2.0d0 - (aref M 2) 2.0d0 - (aref M 3) 3.0d0) - - (let ((V (magicl::make-D-storage 4)) - (D (magicl::make-D-storage 2)) - (lwork 4096) - (liwork 4096) - (info 0) - (eigs-found 0)) - - (lapack:%dsyevr "V" "A" "U" 2 (magicl::copy-matrix-storage M) 2 0.0d0 0.0d0 - 0 0 -1.0d0 eigs-found D V 2 (magicl::make-int32-storage 4) - (magicl::make-D-storage lwork) lwork - (magicl::make-int32-storage liwork) liwork - info) - (format t "M = ~A~%V=~A~%D=~A~%~%" M V D) - - ;; Construct a "matlab-style D" --- is there a better way? - (let ((Df (magicl::make-D-storage 4))) - (setf (aref Df 0) (aref D 0) - (aref Df 3) (aref D 1)) - ;; Reconstruct M as V*Df*V'; - (let ((Mri (magicl::make-D-storage 4)) - (Mr (magicl::make-D-storage 4))) - (blas:%dgemm "N" "N" 2 2 2 1.0d0 V 2 Df 2 0.0d0 Mri 2) - (blas:%dgemm "N" "T" 2 2 2 1.0d0 Mri 2 V 2 0.0d0 Mr 2) - (format t "Reconstructed M = ~A~%" Mr))))))) + (let ((a (magicl:from-list '(1 2 3 4) '(2 2) :type 'double-float))) + (multiple-value-bind (val vect) + (magicl:eig a) + (format t "Eigenvalues = ~A~%Eigenvectors=~A~%~%" val vect)))) (defun qr-example () - (let ((a (make-complex-matrix 3 2 '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2) 4 #C(0 -2.9d0))))) + (let ((a (magicl:from-list '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2)) '(2 2) :type '(complex double-float)))) (qr-printing a))) (defun ql-example () - (let ((a (make-complex-matrix 3 2 '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2) 4 #C(0 -2.9d0))))) + (let ((a (magicl:from-list '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2)) '(2 2) :type '(complex double-float)))) (ql-printing a))) (defun rq-example () - (let ((a (make-complex-matrix 2 3 '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2) 4 #C(0 -2.9d0))))) + (let ((a (magicl:from-list '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2)) '(2 2) :type '(complex double-float)))) (rq-printing a))) (defun lq-example () - (let ((a (make-complex-matrix 2 3 '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2) 4 #C(0 -2.9d0))))) + (let ((a (magicl:from-list '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2)) '(2 2) :type '(complex double-float)))) (lq-printing a))) (defun svd-example () - (let ((a (make-complex-matrix 3 2 '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2) 4 #C(0 -2.9d0))))) + (let ((a (magicl:from-list '(#C(1 2) #C(-4 3) #C(-3 -3) #C(9 2)) '(2 2) :type '(complex double-float)))) (svd-printing a))) (defun qr-printing (a) (multiple-value-bind (q r) - (qr a) - (let ((a-reconst (multiply-complex-matrices q r))) - (format t "A~%") - (print-matrix a) - (format t "Q~%") - (print-matrix q) - (format t "R~%") - (print-matrix r) - (format t "Reconstructed A~%") - (print-matrix a-reconst)))) + (magicl:qr a) + (format t "A~%~a~%" a) + (format t "Q~%~a~%" q) + (format t "R~%~a~%" r) + (let ((a-reconst (magicl:@ q r))) + (format t "Reconstructed A~%~a~%" a-reconst) + (assert (magicl:= a a-reconst +double-comparison-threshold-loose+))))) (defun ql-printing (a) (multiple-value-bind (q l) - (ql a) - (let ((a-reconst (multiply-complex-matrices q l))) - (format t "A~%") - (print-matrix a) - (format t "Q~%") - (print-matrix q) - (format t "L~%") - (print-matrix l) - (format t "Reconstructed A~%") - (print-matrix a-reconst)))) + (magicl:ql a) + (format t "A~%~a~%" a) + (format t "Q~%~a~%" q) + (format t "L~%~a~%" l) + (let ((a-reconst (magicl:@ q l))) + (format t "Reconstructed A~%~a~%" a-reconst) + (assert (magicl:= a a-reconst +double-comparison-threshold-loose+))))) (defun rq-printing (a) (multiple-value-bind (q r) - (rq a) - (let ((a-reconst (multiply-complex-matrices r q))) - (format t "A~%") - (print-matrix a) - (format t "R~%") - (print-matrix r) - (format t "Q~%") - (print-matrix q) - (format t "Reconstructed A~%") - (print-matrix a-reconst)))) + (magicl:rq a) + (format t "A~%~a~%" a) + (format t "Q~%~a~%" q) + (format t "R~%~a~%" r) + (let ((a-reconst (magicl:@ r q))) + (format t "Reconstructed A~%~a~%" a-reconst) + (assert (magicl:= a a-reconst +double-comparison-threshold-loose+))))) (defun lq-printing (a) (multiple-value-bind (q l) - (lq a) - (let ((a-reconst (multiply-complex-matrices l q))) - (format t "A~%") - (print-matrix a) - (format t "L~%") - (print-matrix l) - (format t "Q~%") - (print-matrix q) - (format t "Reconstructed A~%") - (print-matrix a-reconst)))) + (magicl:lq a) + (format t "A~%~a~%" a) + (format t "Q~%~a~%" q) + (format t "L~%~a~%" l) + (let ((a-reconst (magicl:@ l q))) + (format t "Reconstructed A~%~a~%" a-reconst) + (assert (magicl:= a a-reconst +double-comparison-threshold-loose+))))) (defun svd-printing (a) - (multiple-value-bind (u sigma vt) (svd a) - (let* ((rows (matrix-rows a)) - (cols (matrix-cols a)) - (complex-sigma (magicl::make-Z-storage (* rows cols)))) - (dotimes (j cols) - (dotimes (i rows) - ;; TODO: COLUMN-MAJOR-INDEX - (setf (aref complex-sigma (+ (* rows j) i)) - (coerce (ref sigma i j) '(complex double-float))))) - (let ((a-reconst (multiply-complex-matrices - (multiply-complex-matrices u - (make-matrix :rows rows - :cols cols - :data complex-sigma)) - vt))) - (format t "A~%") - (print-matrix a) - (format t "U~%") - (print-matrix u) - (format t "SIGMA~%") - (print-matrix sigma) - (format t "VT~%") - (print-matrix vt) - (format t "Reconstructed A~%") - (print-matrix a-reconst))))) + (multiple-value-bind (u sigma vt) (magicl:svd a) + (let ((complex-sigma (magicl::coerce-type sigma '(complex double-float)))) + (let ((a-reconst (magicl:@ u + complex-sigma + vt))) + (format t "A~%~a~%" a) + (format t "U~%~a~%" u) + (format t "SIGMA~%~a~%" sigma) + (format t "VT~%~a~%" vt) + (format t "Reconstructed A~%~a~%" a-reconst))))) (defun csd-example () - (let ((x (make-complex-matrix 2 2 (list -0.894288 #C(-0.372336 -0.248224) #C(0.372336 -0.248224) -0.894288)))) + (let ((x (magicl:from-list '(-0.894288 #C(0.372336 -0.248224) #C(-0.372336 -0.248224) -0.894288) + '(2 2) + :type '(complex double-float)))) (csd-printing x 1 1))) (defun csd-printing (x p q) (multiple-value-bind (u sigma vt) - (csd x p q) - (let ((x-reconst (multiply-complex-matrices u (multiply-complex-matrices sigma vt)))) - (format t "X~%") - (print-matrix x) - (format t "U~%") - (print-matrix u) - (format t "SIGMA~%") - (print-matrix sigma) - (format t "VT~%") - (print-matrix vt) - (format t "Reconstructed A~%") - (print-matrix x-reconst)))) + (magicl:csd x p q) + (let ((x-reconst (magicl:@ u sigma vt))) + (format t "X~%~a:~a~%" (magicl::layout x) x) + (format t "U~%~a~%" u) + (format t "SIGMA~%~a~%" sigma) + (format t "VT~%~a~%" vt) + (format t "Reconstructed A~%~a~%" x-reconst)))) (defun det-example () - (let* ((x (make-complex-matrix 3 3 (list 2 3 5 #C(1.5 -2) -3 1.5 #C(0 2) 0 #C(0 -3)))) - (d (det x))) + (let* ((x (magicl:from-list (list 2 3 5 #C(1.5 -2) -3 1.5 #C(0 2) 0 #C(0 -3)) '(3 3) :type '(complex double-float))) + (d (magicl:det x))) (format t "X~%") (print-matrix x) (format t "det(X) = ~D~%" d))) (defun inv-example () - (let* ((x (make-complex-matrix 3 3 (list 2 3 5 #C(1.5 -2) -3 1.5 #C(0 2) 0 #C(0 -3)))) - (inv-x (inv x)) - (id (multiply-complex-matrices x inv-x))) - (format t "X~%") - (print-matrix x) - (format t "X^-1~%") - (print-matrix inv-x) - (format t "X*X^-1~%") - (print-matrix id))) + (let* ((x (magicl:from-list (list 2 3 5 #C(1.5 -2) -3 1.5 #C(0 2) 0 #C(0 -3)) '(3 3) :type '(complex double-float))) + (inv-x (magicl:inv x)) + (id (magicl:@ x inv-x))) + (format t "X~%~a~%" x) + (format t "X^-1~%~a~%" inv-x) + (format t "X*X^-1~%~a~%" id))) (defun expm-example () - (let* ((x (make-complex-matrix 4 4 (list 0 0 0 0 0 0 1.5 0 0 -1.5 0 0.5 0 0 -0.5 0))) - (expx (expm x)) - (d (det expx))) - (format t "X~%") - (print-matrix x) - (format t "e^X~%") - (print-matrix expx) + (let* ((x (magicl:from-list (list 0 0 0 0 0 0 1.5 0 0 -1.5 0 0.5 0 0 -0.5 0) '(4 4) :type 'double-float)) + (expx (magicl-transcendental:expm x)) + (d (magicl:det expx))) + (format t "X~%~a~%" x) + (format t "e^X~%~a~%" expx) (format t "det(X) = ~D~%" d))) (defun eig-printing (m) (multiple-value-bind (vals vects) - (eig m) - (let* ((rows (matrix-rows m)) - (val-diag (funcall #'diag rows rows vals))) - (format t "M~%") - (print-matrix m) - (format t "Eigenvalues LAMBDA~%") - (print-matrix val-diag) - (format t "Eigenvectors V~%") - (print-matrix vects) - (format t "M*V~%") - (print-matrix (multiply-complex-matrices m vects)) - (format t "V*LAMBDA~%") - (print-matrix (multiply-complex-matrices vects val-diag))))) + (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))))) (defun eig-example () - (let ((m (make-complex-matrix 3 3 (list -2 -1 1 -2 1 1 -9 -3 4)))) + (let ((m (magicl:from-list (list -2 -1 1 -2 1 1 -9 -3 4) '(3 3) :type 'double-float))) (eig-printing m))) (defun all-examples () diff --git a/magicl-tests.asd b/magicl-tests.asd index f0ed7f68..caab2b0f 100644 --- a/magicl-tests.asd +++ b/magicl-tests.asd @@ -6,7 +6,8 @@ :license "BSD 3-Clause (See LICENSE.txt)" :description "Regression tests for MAGICL." :author "Rigetti Computing" - :depends-on (#:uiop + :depends-on (#:alexandria + #:uiop #:magicl #:magicl-transcendental #:magicl-examples @@ -18,4 +19,9 @@ :serial t :components ((:file "packages") (:file "suite") + (:file "constants") + (:file "util-tests") + (:file "abstract-tensor-tests") + (:file "specialization-tests") + (:file "matrix-tests") (:file "high-level-tests"))) diff --git a/magicl.asd b/magicl.asd index 20e1324a..372f5ec7 100755 --- a/magicl.asd +++ b/magicl.asd @@ -10,8 +10,13 @@ :version (:read-file-form "VERSION.txt") :depends-on (#:alexandria #:cffi - #:cffi-libffi) + #:cffi-libffi + #:abstract-classes + #:policy-cond) :in-order-to ((asdf:test-op (asdf:test-op #:magicl-tests))) + :around-compile (lambda (compile) + (let (#+sbcl (sb-ext:*derive-function-types* t)) + (funcall compile))) :serial t :pathname "src/" :components @@ -20,18 +25,33 @@ (:file "with-array-pointers") (:file "cffi-types") (:module "bindings" - :components - ((:file "lapack00-cffi") - (:file "lapack01-cffi") - (:file "lapack02-cffi") - (:file "lapack03-cffi") - (:file "lapack04-cffi") - (:file "lapack05-cffi") - (:file "lapack06-cffi") - (:file "lapack07-cffi") - (:file "blas-cffi"))) - (:file "high-level/high-level") - (:file "high-level/random") - (:file "high-level/einsum") - (:file "high-level/polynomial-solver") + :components ((:file "lapack00-cffi") + (:file "lapack01-cffi") + (:file "lapack02-cffi") + (:file "lapack03-cffi") + (:file "lapack04-cffi") + (:file "lapack05-cffi") + (:file "lapack06-cffi") + (:file "lapack07-cffi") + (:file "blas-cffi"))) + (:module "high-level" + :serial t + :components ((:file "util") + (:file "shape") + (:file "abstract-tensor") + (:file "specialize-tensor") + (:file "tensor") + (:file "matrix") + (:file "vector") + (:file "lapack-generics") + (:file "types/single-float") + (:file "types/double-float") + (:file "types/complex-single-float") + (:file "types/complex-double-float") + (:file "types/int32") + (:file "lapack-templates") + (:file "lapack-bindings") + (:file "constructors") + (:file "specialize-constructor") + (:file "polynomial-solver"))) (:file "magicl"))) diff --git a/src/high-level/abstract-tensor.lisp b/src/high-level/abstract-tensor.lisp new file mode 100644 index 00000000..9f7a7e4f --- /dev/null +++ b/src/high-level/abstract-tensor.lisp @@ -0,0 +1,265 @@ +;;;; abstract-tensor.lisp +;;;; +;;;; Author: Cole Scott +;;;; +;;;; Collaborators: Robert Smith + +(in-package #:magicl) + +(defstruct (abstract-tensor + (:constructor nil) + (:copier nil)) + "Abstract tensor class. Superclass for implementing the abstract-tensor protocol") + +;;; abstract-tensor protocol +;;; These methods do not have generic definitions and must be implemented by subclasses + +(defgeneric shape (tensor) + (:documentation "The shape (dimensions) of the tensor. eg. '(2 3) for a 2x3 tensor")) + +(defgeneric tref (tensor &rest pos) + (:documentation "Get a reference to element at position")) + +(defgeneric (setf tref) (new-value tensor &rest pos) + (:documentation "Set the value of element at position")) + +(defgeneric copy-tensor (tensor &rest args) + (:documentation "Create a new tensor with the same properties as the given tensor, creating new storage without initializing contents")) + +(defgeneric deep-copy-tensor (tensor &rest args) + (:documentation "Create a new tensor with the same properties as the given tensor, copying the contents of the storage")) + +;;; abstract-tensor generic methods +;;; These methods have default generic implementations but can be optimized by subclasses + +(defgeneric order (tensor) + (:documentation "Order (number of dimensions) of the tensor") + (:method ((tensor abstract-tensor)) + (length (shape tensor)))) + +(defgeneric size (tensor) + (:documentation "Total number of elements in the tensor") + (:method ((tensor abstract-tensor)) + (reduce #'* (shape tensor)))) + +(defgeneric element-type (tensor) + (:documentation "The type of the elements in the tensor") + (:method ((tensor abstract-tensor)) + (declare (ignore tensor)) + t)) + +(defgeneric lisp-array (tensor &optional target) + (:documentation "Return a lisp array containing the data from the tensor + +If TARGET is specified then the contents of the tensor are copied into the array. +In the event TARGET is not specified, the result may return an array sharing memory with the input tensor.") + (:method ((tensor abstract-tensor) &optional target) + (policy-cond:with-expectations (> speed safety) + ((assertion (or (null target) + (and (= (order tensor) (array-rank target)) + (equal (shape tensor) (array-dimensions target)))))) + (let ((arr (or target (make-array (shape tensor) :element-type (element-type tensor))))) + (map-indexes + (shape tensor) + (lambda (&rest pos) + (let ((val (apply #'tref tensor pos))) + (apply #'(setf aref) val arr pos)))) + arr)))) + +(defgeneric map! (function tensor) + (:documentation "Map elements of TENSOR by replacing the value with the output of FUNCTION on the element") + (:method ((function function) (tensor abstract-tensor)) + (map-indexes + (shape tensor) + (lambda (&rest dims) + (apply #'(setf tref) (funcall function (apply #'tref tensor dims)) tensor dims))) + tensor)) + +(defgeneric into! (function tensor) + (:documentation "Map indices of TENSOR by replacing the value with the output of FUNCTION on the index at the element.") + (:method ((function function) (tensor abstract-tensor)) + (map-indexes + (shape tensor) + (lambda (&rest dims) + (apply #'(setf tref) (apply function dims) tensor dims))) + tensor)) + +(defgeneric foreach (function tensor) + (:documentation "Call FUNCTION with each element of TENSOR") + (:method ((function function) (tensor abstract-tensor)) + (map-indexes + (shape tensor) + (lambda (&rest dims) + (funcall function (apply #'tref tensor dims)))) + tensor)) + +(defgeneric map-to (function source target) + (:documentation "Map elements of SOURCE by replacing the corresponding element of TARGET the output of FUNCTION on the source element") + (:method ((function function) (source abstract-tensor) (target abstract-tensor)) + (policy-cond:with-expectations (> speed safety) + ((assertion (equalp (shape source) (shape target)))) + (map-indexes + (shape source) + (lambda (&rest dims) + (apply #'(setf tref) (funcall function (apply #'tref source dims)) target dims)))))) + +(defgeneric map (function tensor) + (:documentation "Map elements of TENSOR, storing the output of FUNCTION on the element into the corresponding element of a new tensor") + (:method ((function function) (tensor abstract-tensor)) + (map! function (deep-copy-tensor tensor)))) + +(defgeneric into (function tensor) + (:documentation "Map indices of TENSOR, storing the output of FUNCTION on the index into the corresponding element of a new tensor + +If LAYOUT is specified then traverse TENSOR in the specified order (column major or row major).") + (:method ((function function) (tensor abstract-tensor)) + (into! function (deep-copy-tensor tensor)))) + +(defgeneric sum (tensor) + (:documentation "Get the sum of the elements of TENSOR") + (:method ((tensor abstract-tensor)) + (let ((sum 0)) + (foreach (lambda (x) (incf sum x)) tensor) + sum))) + +(defgeneric scale (tensor factor) + (:documentation "Scale TENSOR by FACTOR, returning a new tensor of the same type as TENSOR") + (:method ((tensor abstract-tensor) (factor number)) + (map! (lambda (x) (* x factor)) (deep-copy-tensor tensor)))) + +(defgeneric scale! (tensor factor) + (:documentation "Scale TENSOR by FACTOR, storing back into the tensor") + (:method ((tensor abstract-tensor) (factor number)) + (map! (lambda (x) (* x factor)) tensor))) + +(defgeneric slice (tensor from to) + (:documentation "Slice a tensor from FROM to TO, returning a new tensor with the contained elements") + (:method ((tensor abstract-tensor) from to) + (declare (type sequence from to)) + (policy-cond:with-expectations (> speed safety) + ((assertion (and (valid-index-p from (shape tensor)) + (cl:every #'< from (shape tensor)))) + (assertion (and (cl:= (order tensor) (length to)) + (valid-shape-p to) + (cl:every #'<= to (shape tensor)))) + (assertion (cl:every #'<= from to))) + (let* ((dims (mapcar #'- to from)) + (target (empty dims + :layout (layout tensor) + :type (element-type tensor)))) + (map-indexes dims + (lambda (&rest dims) + (setf (apply #'tref target dims) + (apply #'tref tensor (mapcar #'+ dims from))))) + target)))) + +(defgeneric binary-operator (function source1 source2 &optional target) + (:documentation "Perform a binary operator on tensors elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((function function) (source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (policy-cond:with-expectations (> speed safety) + ((assertion (equalp (shape source1) (shape source2)))) + (let ((target (or target (copy-tensor source1)))) + (map-indexes + (shape source1) + (lambda (&rest dims) + (apply #'(setf tref) + (funcall function + (apply #'tref source1 dims) + (apply #'tref source2 dims)) + target dims))) + target)))) + +(defgeneric .+ (source1 source2 &optional target) + (:documentation "Add tensors elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (binary-operator #'+ source1 source2 target))) + +(defgeneric .- (source1 source2 &optional target) + (:documentation "Subtract tensors elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (binary-operator #'- source1 source2 target))) + +(defgeneric .* (source1 source2 &optional target) + (:documentation "Multiply tensors elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (binary-operator #'* source1 source2 target))) + +(defgeneric ./ (source1 source2 &optional target) + (:documentation "Add tensors elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (binary-operator #'/ source1 source2 target))) + +(defgeneric .^ (source1 source2 &optional target) + (:documentation "Exponentiate SOURCE1 by SOURCE2 elementwise, optionally storing the result in TARGET. +If TARGET is not specified then a new tensor is created with the same element type as the first source tensor") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional target) + (binary-operator #'expt source2 source1 target))) + +(defgeneric = (source1 source2 &optional epsilon) + (:documentation "Check the equality of tensors with an optional EPSILON") + (:method ((source1 abstract-tensor) (source2 abstract-tensor) &optional epsilon) + (unless (equal (shape source1) (shape source2)) + (return-from = nil)) + (map-indexes + (shape source1) + (if (null epsilon) + (lambda (&rest pos) + (unless (equal (apply #'tref source1 pos) + (apply #'tref source2 pos)) + (return-from = nil))) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref source1 pos) + (apply #'tref source2 pos))) + epsilon) + (return-from = nil))))) + t)) + +(defgeneric every (predicate tensor &rest more-tensors) + (:documentation "Check that every element in TENSOR meets a given condition") + (:method ((predicate function) (tensor abstract-tensor) &rest more-tensors) + (map-indexes + (shape tensor) + (lambda (&rest dims) + (unless (apply predicate + (apply #'tref tensor dims) + (mapcar (apply #'alexandria:rcurry #'magicl:tref dims) + more-tensors)) + (return-from every nil)))) + t)) + +(defgeneric some (predicate tensor &rest more-tensors) + (:documentation "Check that some elements in TENSOR meet a given condition") + (:method ((predicate function) (tensor abstract-tensor) &rest more-tensors) + (map-indexes + (shape tensor) + (lambda (&rest dims) + (when (apply predicate + (apply #'tref tensor dims) + (mapcar (apply #'alexandria:rcurry #'magicl:tref dims) + more-tensors)) + (return-from some t)))) + nil)) + +(defgeneric notevery (predicate tensor &rest more-tensors) + (:documentation "Check that not every element in TENSOR meets a given condition") + (:method ((predicate function) (tensor abstract-tensor) &rest more-tensors) + (not (apply #'every predicate tensor more-tensors)))) + +(defgeneric notany (predicate tensor &rest more-tensors) + (:documentation "Check that not any element in TENSOR meets a given condition") + (:method ((predicate function) (tensor abstract-tensor) &rest more-tensors) + (not (apply #'some predicate tensor more-tensors)))) + +(defgeneric coerce-type (tensor type) + (:documentation "Coerce element type of TENSOR to TYPE by creating a new tensor") + (:method ((tensor abstract-tensor) type) + (let ((target (empty (shape tensor) :layout (layout tensor) :type type))) + (map-to (lambda (x) (coerce x type)) + tensor + target) + target))) diff --git a/src/high-level/constructors.lisp b/src/high-level/constructors.lisp new file mode 100644 index 00000000..21b1efda --- /dev/null +++ b/src/high-level/constructors.lisp @@ -0,0 +1,177 @@ +;;;; constructors.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(defvar *default-tensor-type* 'double-float + "The default element type to be used for constructing tensors when TYPE is not specified.") + +(defgeneric make-tensor (class shape &key initial-element layout storage) + (:documentation "Make a dense tensor with elements of the specified type")) + +(defun empty (shape &key (type *default-tensor-type*) layout) + "Create a tensor without intializing the contents of the storage + +If TYPE is not specified then *DEFAULT-TENSOR-TYPE* is used. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor specialized on the specified SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type shape shape)) + (let ((tensor-type (infer-tensor-type type shape nil))) + (make-tensor tensor-type shape :layout layout)))) + +(defun const (const shape &key type layout) + "Create a tensor with the specified SHAPE with each element being set to CONST + +If TYPE is not specified then it is inferred from the type of CONST. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type shape shape)) + (let ((tensor-class (infer-tensor-type type shape const))) + (make-tensor tensor-class shape :layout layout :initial-element const)))) + +(defun rand (shape &key (type *default-tensor-type*) layout distribution) + "Create tensor with random elements from DISTRIBUTION + +DISTRIBUTION is a function with no arguments which returns a value for the element. +If DISTRIBUTION is not specified then a uniform distribution on [0,1] (or [0,1] + [0,1]i for complex types) is used. +If TYPE is not specified then *DEFAULT-TENSOR-TYPE* is used. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type shape shape)) + (let* ((tensor-class (infer-tensor-type type shape nil)) + (rand-function + (or distribution + (cond + ((subtypep type 'complex) + (lambda () + (complex + (random 1d0) + (random 1d0)))) + (t + (lambda () + (random 1d0)))))) + (f (lambda (&rest rest) + (declare (ignore rest)) + (coerce (funcall rand-function) type)))) + (into! f (make-tensor tensor-class shape :layout layout))))) + +(defun eye (shape &key value (offset 0) (type *default-tensor-type*) layout) + "Create an identity tensor + +SHAPE can either be a list of dimensions or a fixnum defining the length of the side a square matrix. + +If VALUE is not specified then 1 is used. +If TYPE is not specified then it is inferred from the type of VALUE, defaulting to *DEFAULT-TENSOR-TYPE*. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type (or shape fixnum) shape) + (type fixnum offset)) + (let* ((shape (if (integerp shape) + (fixnum-to-shape shape) + shape)) + (tensor-class (infer-tensor-type (if value nil type) shape value)) + (tensor (make-tensor tensor-class shape :layout layout)) + (shape-length (length shape)) + (fill-value (coerce (or value 1) (element-type tensor)))) + (loop :for i :below (reduce #'min shape) + :do (setf (apply #'tref tensor (make-list shape-length :initial-element i)) fill-value)) + tensor))) + +(defun arange (range &key (type *default-tensor-type*)) + "Create a 1-dimensional tensor of elements from 0 up to but not including the RANGE + +If TYPE is not specified then it is inferred from the type of RANGE. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on TYPE with shape (floor(RANGE))." + (let ((tensor-class (infer-tensor-type type (list (floor range)) range))) + (let ((tensor (make-tensor tensor-class (list (floor range)))) + (f (lambda (index) + (coerce index type)))) + (into! f tensor)))) + +;; NOTE: This does not accout for changing array layout. If the input +;; is in row-major but the constructor specifies column-major +;; then the tensor will not return the correct values. Not a +;; huge deal but something that be improved (similarly to how +;; FROM-LIST does it.) +(defun from-array (array shape &key type (layout :row-major)) + "Create a tensor from ARRAY, calling ADJUST-ARRAY on ARRAY to flatten to a 1-dimensional array of length equal to the product of the elements in SHAPE + +If TYPE is not specified then it is inferred from the element type of ARRAY. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type shape shape)) + (let* (;; TODO: This gets the type of an array from the array, not + ;; by inspecting the elements. This can cause issues when + ;; the element-type of the array is T + (element-type + (if (null type) + (array-element-type array) + type)) + (tensor-class (infer-tensor-type element-type shape nil))) + (adjust-array array (list (reduce #'* shape)) :element-type element-type) + (make-tensor tensor-class shape + :storage array + :layout layout)))) + +(defun from-list (list shape &key type layout (input-layout :row-major)) + "Create a tensor with the elements of LIST, placing in layout INPUT-LAYOUT + +If INPUT-LAYOUT is not specified then row-major is assumed. + +If TYPE is not specified then it is inferred from the type of the first element of LIST. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (policy-cond:with-expectations (> speed safety) + ((type shape shape) + (assertion (cl:= (length list) (reduce #'* shape)))) + (let* ((tensor-class (infer-tensor-type type shape (first list))) + (tensor (make-tensor tensor-class shape :layout layout))) + (into! + (lambda (&rest pos) + (nth + (if (eql input-layout :row-major) + (row-major-index pos shape) + (column-major-index pos shape)) + list)) + tensor)))) + +(defun from-diag (list &key (order 2) type layout) + "Create a tensor from a list, placing along the diagonal + +If ORDER is specified then the tensor will be of that order, otherwise 2 is assumed. +If TYPE is not specified then it is inferred from the type of the first element of LIST. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor is specialized on SHAPE and TYPE." + (let* ((length (length list)) + (shape (fixnum-to-shape length order)) + (tensor-class (infer-tensor-type type shape (first list))) + (tensor (make-tensor tensor-class shape :layout layout))) + (loop :for i :below length + :do (setf (apply #'tref tensor (make-list order :initial-element i)) + (pop list))) + tensor)) + +;;; Constructors for convenience + +(defun zeros (shape &key (type *default-tensor-type*) layout) + "Create a tensor with the specified SHAPE of zeros + +If TYPE is not specified then *DEFAULT-TENSOR-TYPE* is used. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor specialized on the specified SHAPE and TYPE." + (const 0 shape :type type :layout layout)) + +(defun ones (shape &key (type *default-tensor-type*) layout) + "Create a tensor with the specified SHAPE of ones + +If TYPE is not specified then *DEFAULT-TENSOR-TYPE* is used. +LAYOUT specifies the internal storage representation ordering of the returned tensor. +The tensor specialized on the specified SHAPE and TYPE." + (const 1 shape :type type :layout layout)) diff --git a/src/high-level/high-level.lisp b/src/high-level/high-level.lisp deleted file mode 100644 index fb0d9b9c..00000000 --- a/src/high-level/high-level.lisp +++ /dev/null @@ -1,1031 +0,0 @@ -;;;; high-level.lisp -;;;; -;;;; Author: Joseph Lin -;;;; Robert Smith - -(in-package #:magicl) - -(deftype matrix-storage () - "Representation of the smallest supertype of matrix storage." - `(simple-array * (*))) - -(deftype matrix-dimension () - "Representation of a valid dimension of a matrix." - `(integer 1 ,array-total-size-limit)) - -(deftype matrix-index () - "Representation of a valid matrix index." - `(integer 0 (,array-total-size-limit))) - -(defun make-int32-storage (n) - (make-array n :element-type '(signed-byte 32) :initial-element 0)) - -(defun make-S-storage (n) - (make-array n :element-type 'single-float :initial-element 0.0f0)) - -(defun make-D-storage (n) - (make-array n :element-type 'double-float :initial-element 0.0d0)) - -(defun make-C-storage (n) - (make-array n :element-type '(complex single-float) :initial-element #C(0.0f0 0.0f0))) - -(defun make-Z-storage (n) - (make-array n :element-type '(complex double-float) :initial-element #C(0.0d0 0.0d0))) - -(deftype lapack-data-type () - "The data type symbols used in LAPACK." - `(member S D C Z)) - -(deftype lapack-storage-type () - "The storage type of the data for LAPACK." - `(member - ;; general - ge gb gt gg - ;; symmetric - sy sb sp st - ;; Hermitian - he hb hp - ;; SPD / HPD - po pb pp pt - ;; triangular - tr tb tp tg - ;; upper Hessenberg - hs hg - ;; trapezoidal - tz - ;; orthogonal - or op - ;; unitary - un up - ;; diagonal - di - ;; bidiagonal - bd)) - -(defstruct (matrix (:constructor %make-matrix (rows cols data-type data))) - "Representation of a dense matrix." - (rows (error "Required argument") - :type matrix-dimension - :read-only t) - (cols (error "Required argument") - :type matrix-dimension - :read-only t) - (data-type (error "Required argument") - :type lapack-data-type - :read-only t) - (storage-type 'ge - :type lapack-storage-type - :read-only t) - (data (error "Required argument") - :type matrix-storage - :read-only t)) - -;; Allow MATRIX objects to be dumped. -(defmethod make-load-form ((object matrix) &optional environment) - (make-load-form-saving-slots object :slot-names '(rows cols data-type storage-type data) - :environment environment)) - -(defun matrix-storage-size (v) - "Compute the size (i.e., number of elements) of the matrix storage vector V." - (declare (type matrix-storage v)) - (length v)) - -(defun make-matrix (&key rows cols data) - (declare (type matrix-dimension rows cols) - (type matrix-storage data)) - (assert (= (matrix-storage-size data) - (* rows cols)) - (rows cols) - "The number of rows=~D and cols=~D does not match the storage size=~D" - rows - cols - (matrix-storage-size data)) - (let ((data-type (etypecase data - ((array single-float) 'S) - ((array double-float) 'D) - ((array (complex single-float)) 'C) - ((array (complex double-float)) 'Z)))) - (%make-matrix rows cols data-type data))) - -(defun copy-matrix-storage (v) - "Copy a matrix storage vector V suitable for the DATA slot on a MATRIX structure." - (declare (type matrix-storage v)) - (copy-seq v)) - -(defun matrix-element-type (m) - "Return the element type of the matrix M." - (array-element-type (matrix-data m))) - -(defun pprint-matrix (stream matrix) - "Pretty-print a matrix MATRIX to the stream STREAM." - (flet ((print-real (x) - (format stream "~6,3f" x)) - (print-complex (z) - (format stream "~6,3f ~:[+~;-~]~6,3fj" - (realpart z) - (minusp (imagpart z)) - (abs (imagpart z))))) - (let* ((rows (matrix-rows matrix)) - (cols (matrix-cols matrix)) - (type (matrix-data-type matrix)) - (print-entry (ecase type - ((S D) #'print-real) - ((C Z) #'print-complex)))) - (pprint-logical-block (stream nil) - (print-unreadable-object (matrix stream :type t) - (format stream "[~A] ~Dx~D:" (symbol-name type) rows cols) - (dotimes (r rows) - (pprint-newline :mandatory stream) - (dotimes (c cols) - (funcall print-entry (ref matrix r c)) - (unless (= c (1- cols)) - (write-string " " stream))))))))) - -(set-pprint-dispatch 'matrix 'pprint-matrix) - -(defun make-complex-foreign-vector (entries) - "Makes a complex double FNV out ENTRIES, a list of complex numbers." - (let* ((len (length entries)) - (v (make-Z-storage len))) - (loop :for i :below len - :for e :in entries - :do (setf (aref v i) (coerce e '(complex double-float))) - :finally (return v)))) - -(defun vector-to-list (v) - "Make a list from a fnv." - (coerce v 'list)) - -(defun make-complex-vector (entries) - "Makes a complex column vector out of ENTRIES, a list of complex numbers." - (funcall #'make-complex-matrix (length entries) 1 entries)) - -(defun make-complex-matrix (m n entries) - "Makes an M-by-N matrix assuming ENTRIES is a list of complex numbers in column major order." - (check-type m matrix-dimension) - (check-type n matrix-dimension) - (let ((entries-size (length entries)) - (expected-size (* m n))) - (assert (= entries-size expected-size) - () - "Length of entries is ~D, is not ~D * ~D = ~D" - entries-size m n expected-size) - (make-matrix :rows m - :cols n - :data (funcall #'make-complex-foreign-vector entries)))) - -(defun make-zero-matrix (rows cols) - (make-matrix - :rows rows - :cols cols - :data (make-Z-storage (* rows cols)))) - -(defun square-matrix-p (m) - "Is M a square matrix?" - (= (matrix-rows m) - (matrix-cols m))) - -(defparameter *default-zero-comparison-epsilon* 1d-10 - "The default absolute radius about zero considered to still be zero.") - -(defun identityp (m) - "Is the matrix M an identity matrix up to the epsilon *DEFAULT-ZERO-COMPARISON-EPSILON*?" - ;; Sorry, this function is a little bit convoluted. - (flet ((within-epsilon (z) - (<= (abs z) *default-zero-comparison-epsilon*))) - (and (square-matrix-p m) - (progn - (map-indexes (matrix-rows m) - (matrix-cols m) - (lambda (r c) - (cond - ;; Diagonal is 1. - ((= r c) - (when (not (within-epsilon (- 1.0d0 (ref m r c)))) - (return-from identityp nil))) - ;; Everything else is 0. - ((not (within-epsilon (ref m r c))) - (return-from identityp nil))))) - t)))) - -(defun unitaryp (m) - "Is the matrix M is unitary up to the epsilon *DEFAULT-ZERO-COMPARISON-EPSILON*." - (if (numberp m) - (<= (abs (- 1.0 (abs m))) *default-zero-comparison-epsilon*) - (identityp (multiply-complex-matrices m (conjugate-transpose m))))) - -(defun map-indexes (rows cols f) - "Map the binary function F(i, j) across all matrix indexes 0 <= i < ROWS and 0 <= j < COLS. Return NIL." - (dotimes (c cols) - (dotimes (r rows) - (funcall f r c)))) - -(defun tabulate (rows cols f) - "Generate a matrix of ROWS x COLS by calling the function F(i, j) for 0 <= i < ROWS and 0 <= j < COLS." - (let ((m (make-zero-matrix rows cols))) - (map-indexes rows cols (lambda (r c) (setf (ref m r c) (funcall f r c)))) - m)) - -(defun lift-unary-function (function) - "Produces a unitary function that takes a matrix and returns a -matrix of the same dimension where each element of the output matrix -is the result of applying the unitary FUNCTION to the corresponding -element in the input matrix." - (check-type function function) - (lambda (matrix) - (check-type matrix matrix) - (tabulate (matrix-rows matrix) (matrix-cols matrix) - (lambda (i j) (funcall function (ref matrix i j)))))) - -(setf (symbol-function 'inc-matrix) (lift-unary-function #'1+)) -(setf (documentation #'inc-matrix 'function) - "Returns matrix with each element + 1.") -(setf (symbol-function 'dec-matrix) (lift-unary-function #'1-)) -(setf (documentation #'dec-matrix 'function) - "Returns matrix with each element - 1.") - -(defun lift-binary-function (function) - "Produces a binary function that takes a matrix and returns a -matrix of the same dimension where each element of the output matrix -is the result of applying the binary FUNCTION to the corresponding -elements in the input matrices." - (check-type function function) - (lambda (a b) - (check-type a matrix) - (check-type b matrix) - (tabulate (matrix-rows a) (matrix-cols a) - (lambda (i j) (funcall function - (ref a i j) - (ref b i j)))))) - -(defgeneric add-matrix-generic (a b)) -(defgeneric sub-matrix-generic (a b)) - -(let ((lifted-+ (lift-binary-function #'+)) - (lifted-- (lift-binary-function #'-))) - (defmethod add-matrix-generic ((a matrix) (b matrix)) - (funcall lifted-+ a b)) - - (defmethod add-matrix-generic ((a number) (b matrix)) - (let ((mat-a (tabulate (matrix-rows b) (matrix-cols b) - (lambda (i j) (declare (ignore i j)) a)))) - (add-matrix-generic mat-a b))) - - (defmethod add-matrix-generic ((b matrix) (a number)) - (add-matrix-generic a b)) - - (defmethod sub-matrix-generic ((a matrix) (b matrix)) - (funcall lifted-- a b)) - - (defmethod sub-matrix-generic ((a number) (b matrix)) - (let ((mat-a (tabulate (matrix-rows b) (matrix-cols b) - (lambda (i j) (declare (ignore i j)) a)))) - (sub-matrix-generic mat-a b))) - - (defmethod sub-matrix-generic ((b matrix) (a number)) - (let ((mat-a (tabulate (matrix-rows b) (matrix-cols b) - (lambda (i j) (declare (ignore i j)) a)))) - (sub-matrix-generic b mat-a)))) - -(defun add-matrix (matrix &rest more-matrices) - "Element-wise addition of input matrices." - (reduce #'add-matrix-generic more-matrices :initial-value matrix)) - -(defun sub-matrix (matrix &rest more-matrices) - "Element-wise subtraction of input matrices." - (reduce #'sub-matrix-generic more-matrices :initial-value matrix)) - -(defun make-identity-matrix (dimension) - "Make an identity matrix of dimension DIMENSION." - (tabulate dimension dimension (lambda (i j) (if (= i j) 1 0)))) - -(defun diag (m n entries) - "Creates a matrix with ENTRIES along the diagonal" - (let ((entries-size (length entries)) - (expected-size (min m n))) - (assert (= entries-size expected-size) () - "Min dimension is ~S but number of entries is ~S" expected-size entries-size) - (let ((mat (make-zero-matrix m n))) - (dotimes (i entries-size mat) - (setf (ref mat i i) (nth i entries)))))) - -(defun matrix-diagonal (m) - "Get the diagonal elements of the matrix M as a list." - (loop :for i :below (min (matrix-rows m) (matrix-cols m)) - :collect (ref m i i))) - -(declaim (inline column-major-index) - (ftype (function (matrix-dimension matrix-index matrix-index) matrix-index) - column-major-index)) -(defun column-major-index (rows i j) - "Give the linear index of a matrix element addressed by the I'th row and J'th column, for a matrix with ROWS rows." - (+ i (* j rows))) - -(defun ref (m i j) - "Accessor method for the element in the I-th row and J-th column of a matrix M, assuming zero indexing." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m)) - (data (matrix-data m))) - (check-type i matrix-index) - (check-type j matrix-index) - (assert (< -1 i rows) () "row index ~D is out of range" i) - (assert (< -1 j cols) () "col index ~D is out of range" j) - (aref data (column-major-index rows i j)))) - -(defun ptr-ref (m base i j) - "Accessor method for the pointer to the element in the I-th row and J-th column of a matrix M, assuming zero indexing." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m)) - (data (matrix-data m))) - (check-type i matrix-index) - (check-type j matrix-index) - (assert (< -1 i rows) () "row index ~D is out of range" i) - (assert (< -1 j cols) () "col index ~D is out of range" j) - (let ((idx (column-major-index rows i j))) - (etypecase (row-major-aref data 0) - (single-float (cffi:mem-aptr base :float idx)) - (double-float (cffi:mem-aptr base :double idx)) - ((complex single-float) (cffi:mem-aptr base :float (* 2 idx))) - ((complex double-float) (cffi:mem-aptr base :double (* 2 idx))))))) - -(defparameter *type-strictness* nil - "Be strict about types when setting values in a matrix.") - -(defun (setf ref) (new-value m i j) - "Set the value of M_IJ to NEW-VALUE." - (let* ((rows (matrix-rows m)) - (cols (matrix-cols m)) - (data (matrix-data m)) - (type (array-element-type data))) - (check-type i integer) - (check-type j integer) - (assert (< -1 i rows) () "row index ~S is out of range" i) - (assert (< -1 j cols) () "col index ~S is out of range" j) - (setf (aref data (column-major-index rows i j)) - (if *type-strictness* - new-value - (coerce new-value type))))) - -(defun multiply-complex-matrices (ma mb) - "Multiplies two complex marices MA and MB, returning MA*MB. If MA is M x KA and MB is KB x N, -it must be that KA = KB, and the resulting matrix is M x N." - (assert (subtypep (matrix-element-type ma) (matrix-element-type mb))) - (assert (subtypep (matrix-element-type mb) (matrix-element-type ma))) - (assert (subtypep (matrix-element-type ma) '(complex double-float))) - (let ((m (matrix-rows ma)) - (ka (matrix-cols ma)) - (kb (matrix-rows mb)) - (n (matrix-cols mb))) - (assert (= ka kb) () - "Matrix A has ~D columns while matrix B has ~D rows" ka kb) - (let ((a (matrix-data ma)) - (b (matrix-data mb)) - (transa "N") - (transb "N") - (alpha #C(1.0d0 0.0d0)) - (beta #C(0.0d0 0.0d0)) - (c (make-Z-storage (* m n)))) - (magicl.blas-cffi::%zgemm transa transb m n ka alpha a m b kb beta c m) - (make-matrix :rows m :cols n :data c)))) - -(defun conjugate-entrywise (m) - "Computes the conjugate of each entry of matrix M." - (check-type m matrix) - (let ((m* (make-zero-matrix (matrix-cols m) - (matrix-rows m)))) - (dotimes (i (matrix-rows m) m*) - (dotimes (j (matrix-cols m)) - (setf (ref m* i j) (conjugate (ref m i j))))))) - -(defun transpose (m) - "Computes the transpose of the matrix M." - (check-type m matrix) - (let ((mT (make-zero-matrix (matrix-cols m) - (matrix-rows m)))) - (dotimes (i (matrix-rows m) mT) - (dotimes (j (matrix-cols m)) - (setf (ref mT j i) (ref m i j)))))) - -(defun conjugate-transpose (m) - "Computes the conjugate transpose of the matrix M." - (check-type m matrix) - (let ((mT (make-zero-matrix (matrix-cols m) - (matrix-rows m)))) - (dotimes (i (matrix-rows m) mT) - (dotimes (j (matrix-cols m)) - (setf (ref mT j i) (conjugate (ref m i j))))))) - - -(defun scale (alpha x) - "Scale a complex double matrix X by a complex double ALPHA, i.e. return ALPHA*X." - (let* ((zx (copy-matrix-storage (matrix-data x))) - (za (coerce alpha '(complex double-float))) - (n (matrix-storage-size zx))) - (magicl.blas-cffi::%zscal n za zx 1) - (make-matrix :rows (matrix-rows x) :cols (matrix-cols x) :data zx))) - -(defun qr (m) - "Finds the QR factorization of the matrix M." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m))) - (multiple-value-bind (a tau) (lapack-unitary-triangular-decomposition m "QR") - (let* ((amat (make-matrix :rows rows :cols cols :data a)) - (r (get-square-triangular amat T cols)) - (q (unitary-triangular-helper-get-q amat tau "QR"))) - ;; change signs if diagonal elements of r are negative - (dotimes (j cols) - (let ((diag-elt (ref r j j))) - (assert (zerop (imagpart diag-elt)) - () "Diagonal element R_~D~D=~A is not real" j j diag-elt) - (setf diag-elt (realpart diag-elt)) - (when (minusp diag-elt) - (dotimes (i rows) - (when (<= j i (1- cols)) - (setf (ref r j i) (- (ref r j i)))) - (setf (ref q i j) (- (ref q i j))))))) - (values q r))))) - -(defun ql (m) - "Finds the QL factorization of the matrix M." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m))) - (multiple-value-bind (a tau) (lapack-unitary-triangular-decomposition m "QL") - (let* ((amat (make-matrix :rows rows :cols cols :data a)) - (l (get-square-triangular amat nil cols)) - (q (unitary-triangular-helper-get-q amat tau "QL"))) - ;; change signs if diagonal elements of L are negative - (dotimes (j cols) - (let ((diag-elt (ref l j j))) - (assert (zerop (imagpart diag-elt)) - () "Diagonal element L_~D~D=~A is not real" j j diag-elt) - (setf diag-elt (realpart diag-elt)) - (when (minusp diag-elt) - (dotimes (i rows) - (when (<= i j) - (setf (ref l j i) (- (ref l j i)))) - (setf (ref q i j) (- (ref q i j))))))) - (values q l))))) - -(defun rq (m) - "Finds the RQ factorization of the matrix M." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m))) - (multiple-value-bind (a tau) (lapack-unitary-triangular-decomposition m "RQ") - (let* ((amat (make-matrix :rows rows :cols cols :data a)) - (r (get-square-triangular amat T rows)) - (q (unitary-triangular-helper-get-q amat tau "RQ"))) - ;; change signs if diagonal elements of r are negative - (dotimes (i rows) - (let ((diag-elt (ref r i i))) - (assert (zerop (imagpart diag-elt)) - () "Diagonal element R_~D~D=~A is not real" i i diag-elt) - (setf diag-elt (realpart diag-elt)) - (when (minusp diag-elt) - (dotimes (j cols) - (when (<= j i) - (setf (ref r j i) (- (ref r j i)))) - (setf (ref q i j) (- (ref q i j))))))) - (values q r))))) - -(defun lq (m) - "Finds the LQ factorization of the matrix M." - (let ((rows (matrix-rows m)) - (cols (matrix-cols m))) - (multiple-value-bind (a tau) (lapack-unitary-triangular-decomposition m "LQ") - (let* ((amat (make-matrix :rows rows :cols cols :data a)) - (r (get-square-triangular amat nil rows)) - (q (unitary-triangular-helper-get-q amat tau "LQ"))) - ;; change signs if diagonal elements of l are negative - (dotimes (i rows) - (let ((diag-elt (ref r i i))) - (assert (zerop (imagpart diag-elt)) - () "Diagonal element R_~D~D=~A is not real" i i diag-elt) - (setf diag-elt (realpart diag-elt)) - (when (minusp diag-elt) - (dotimes (j cols) - (when (<= i j (1- rows)) - (setf (ref r j i) (- (ref r j i)))) - (setf (ref q i j) (- (ref q i j))))))) - (values q r))))) - -(defun lapack-unitary-triangular-decomposition (m option) - "Finds the QR/QL/RQ/LQ factorization of the matrix M to the intermediate representation, as given by the LAPACK ZGE(QR/QL/RQ/LQ)F subroutine, depending on the string option OPTION." - (let ((lapack-func (alexandria:eswitch (option :test #'string=) - ("QR" #'magicl.lapack-cffi::%zgeqrf) - ("QL" #'magicl.lapack-cffi::%zgeqlf) - ("RQ" #'magicl.lapack-cffi::%zgerqf) - ("LQ" #'magicl.lapack-cffi::%zgelqf))) - (rows (matrix-rows m)) - (cols (matrix-cols m)) - (a (copy-matrix-storage (matrix-data m))) - (lwork -1) - (info 0)) - (let ((lda rows) - (tau (make-Z-storage (min rows cols))) - (work (make-Z-storage (max 1 lwork)))) - ;; run it once as a workspace query - (funcall lapack-func rows cols a lda tau work lwork info) - (setf lwork (round (realpart (aref work 0)))) - (setf work (make-Z-storage (max 1 lwork))) - ;; run it again with optimal workspace size - (funcall lapack-func rows cols a lda tau work lwork info) - (values a tau)))) - -(defun get-square-triangular (a upper ord) - "Creates a square matrix from a triangular or trapezoidal portion of A. The square matrix is upper triangular and taken from the upper portion of A if and only if UPPER is T. The order of the square matrix is given by ORD." - (check-type upper boolean) - (check-type ord (integer 1 *) "a positive integer") - (let ((m (matrix-rows a)) - (n (matrix-cols a))) - (assert (<= ord (max m n)) () "ORD, given as ~D, is greater than the maximum dimension of A, ~D." ord (max m n)) - - (let ((tri (make-zero-matrix ord ord))) - (if (> m n) - (if upper - (loop for i from 0 to (1- ord) - do (loop for j from (max 0 (+ (- n ord) i)) to (1- n) - do (setf (ref tri i (+ j (- ord n))) (ref a i j)))) - (loop for i from (- m ord) to (1- m) - do (loop for j from 0 to (min (+ (- ord m) i) (1- n)) - do (setf (ref tri (- i (- m ord)) j) (ref a i j))))) - (if upper - (loop for j from (- n ord) to (1- n) - do (loop for i from 0 to (min (+ (- ord n) j) (1- m)) - do (setf (ref tri i (- j (- n ord))) (ref a i j)))) - (loop for j from 0 to (1- ord) - do (loop for i from (max 0 (+ (- m ord) j)) to (1- m) - do (setf (ref tri (+ i (- ord m)) j) (ref a i j)))))) - (values tri)))) - -(defun unitary-triangular-helper-get-q (amat tau option) - "Finds the unitary matrix Q from QR/QL/RQ/LQ factorization of the matrix M, given the reflectors and intermediate representation provided by the LAPACK ZGE(QR/QL/RQ/LQ)F subroutine. Whether the factorization is QR, QL, RQ or LQ is given by string option OPTION." - (let ((lapack-func (alexandria:eswitch (option :test #'string=) - ("QR" #'magicl.lapack-cffi::%zungqr) - ("QL" #'magicl.lapack-cffi::%zungql) - ("RQ" #'magicl.lapack-cffi::%zungrq) - ("LQ" #'magicl.lapack-cffi::%zunglq))) - (m (matrix-rows amat)) - (n (matrix-cols amat)) - (a (matrix-data amat)) - (k (matrix-storage-size tau)) - (lwork -1) - (info 0)) - (let ((lda m) - (work (make-Z-storage (max 1 lwork)))) - ;; run it once as a workspace query - (funcall lapack-func m n k a lda tau work lwork info) - (setf lwork (round (realpart (aref work 0)))) - (setf work (make-Z-storage (max 1 lwork))) - ;; run it again with optimal workspace size - (funcall lapack-func m n k a lda tau work lwork info) - (make-matrix :rows m :cols n :data a)))) - -(defun qr-helper-get-q (a tau n) - "Get the matrix Q as a product of reflectors, from results given by ZGEQRF." - (let ((m (/ (matrix-storage-size a) n)) - (k (matrix-storage-size tau)) - (lwork -1) - (info 0)) - (let ((lda m) - (work (make-Z-storage (max 1 lwork)))) - ;; run it once as a workspace query - (magicl.lapack-cffi::%zungqr m n k a lda tau work lwork info) - (setf lwork (truncate (realpart (aref work 0)))) - (setf work (make-Z-storage (max 1 lwork))) - ;; run it again with optimal workspace size - (magicl.lapack-cffi::%zungqr m n k a lda tau work lwork info) - (make-matrix :rows m :cols n :data a)))) - -(defun svd (m &key reduced) - "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt. If REDUCED is non-NIL, return the reduced SVD (where either U or V are just partial isometries and not necessarily unitary matrices)." - (let* ((jobu (if reduced "S" "A")) - (jobvt (if reduced "S" "A")) - (rows (matrix-rows m)) - (cols (matrix-cols m)) - (a (copy-matrix-storage (matrix-data m))) - (lwork -1) - (info 0) - (k (min rows cols)) - (u-cols (if reduced k rows)) - (vt-rows (if reduced k cols)) - (lda rows) - (s (make-D-storage k)) - (ldu rows) - (ldvt vt-rows) - (work1 (make-Z-storage (max 1 lwork))) - (work nil) - (rwork (make-D-storage (* 5 (min rows cols)))) - (u (make-Z-storage (* ldu u-cols))) - (vt (make-Z-storage (* ldvt cols)))) - ;; run it once as a workspace query - (magicl.lapack-cffi::%zgesvd jobu jobvt rows cols a lda s u ldu vt ldvt - work1 lwork rwork info) - (setf lwork (round (realpart (aref work1 0)))) - (setf work (make-Z-storage (max 1 lwork))) - ;; run it again with optimal workspace size - (magicl.lapack-cffi::%zgesvd jobu jobvt rows cols a lda s u ldu vt ldvt - work lwork rwork info) - (let ((smat (make-D-storage (* u-cols vt-rows)))) - (dotimes (i k) - (setf (aref smat (column-major-index u-cols i i)) - (aref s i))) - (values (make-matrix :rows rows :cols u-cols :data u) - (make-matrix :rows u-cols :cols vt-rows :data smat) - (make-matrix :rows vt-rows :cols cols :data vt))))) - -(defun slice (m rmin rmax cmin cmax) - "Get the subarray of M containing all elements M_IJ, where RMIN<=I 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)))) + (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)) + (case transa + (: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)) + (case transb + (:n :t) + (:t :n) + (:c (error "Specifying TRANSB to be :C is not allowed if B is ROW-MAJOR"))) + transb)) + (target (or target + (empty + (list m n) + :type ',type)))) + (,blas-function + (ecase ta + (:t "T") + (:c "C") + (:n "N")) + (ecase tb + (:t "T") + (:c "C") + (:n "N")) + m + n + k + alpha + (storage a) + (if (eql :n ta) m k) + (storage b) + (if (eql :n tb) k n) + beta + (storage target) + m) + target)))))) + +(defun generate-lapack-lu-for-type (class type lu-function) + (declare (ignore type)) + `(progn + (defmethod lu ((m ,class)) + (lapack-lu m)) + + (defmethod lapack-lu ((a ,class)) + (let* ((a-tensor (deep-copy-tensor a)) + (a (storage a-tensor)) + (m (nrows a-tensor)) + (n (ncols a-tensor)) + (lda m) + (ipiv-tensor (empty (list (max m n)) :type '(signed-byte 32))) + (ipiv (storage ipiv-tensor)) + (info 0)) + (when (eql :row-major (layout a-tensor)) (transpose! a-tensor)) + (,lu-function + m + n + a + lda + ipiv + info) + (values a-tensor ipiv-tensor))))) + +(defun generate-lapack-inv-for-type (class type lu-function inv-function) + `(progn + (defmethod inv ((m ,class)) + (lapack-inv m)) + + (defmethod lapack-inv ((a ,class)) + (let ((a-tensor (deep-copy-tensor a))) + (when (eql :row-major (layout a-tensor)) (transpose! a-tensor)) + (let* ((a (storage a-tensor)) + (m (nrows a-tensor)) + (n (ncols a-tensor)) + (lda m) + (ipiv-tensor (empty (list (max m n)) :type '(signed-byte 32))) + (ipiv (storage ipiv-tensor)) + (info 0)) + (,lu-function + m + n + a + lda + ipiv + info) + ;; TODO: This check is already performed by the LU + ;; function, however INFO is not being returned + ;; correctly. When the bindings are fixed this should + ;; just check if INFO is non-zero + (assert (cl:notany + (lambda (x) + (cl:= x 0)) + (diag a-tensor)) + () "The provided matrix is singular and cannot be inverted.") + (let* ((lwork -1) + (work1 (make-array (max 1 lwork) :element-type ',type)) + (work nil) + (info 0)) + ;; Perform work size query with work of length 1 + (,inv-function + n + a + m + ipiv + work1 + lwork + info) + (setf lwork (round (realpart (aref work1 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; Perform actual operation with correct work size + (,inv-function + n + a + m + ipiv + work + lwork + info)) + (values a-tensor)))))) + +(defun generate-lapack-svd-for-type (class type svd-function &optional real-type) + `(progn + (defmethod svd ((m ,class) &key reduced) + (lapack-svd m :reduced reduced)) + + (defmethod lapack-svd ((m ,class) &key reduced) + "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt. If REDUCED is non-NIL, return the reduced SVD (where either U or V are just partial isometries and not necessarily unitary matrices)." + (let* ((jobu (if reduced "S" "A")) + (jobvt (if reduced "S" "A")) + (rows (nrows m)) + (cols (ncols m)) + (a (alexandria:copy-array (storage (if (eql :row-major (layout m)) (transpose m) m)))) + (lwork -1) + (info 0) + (k (min rows cols)) + (u-cols (if reduced k rows)) + (vt-rows (if reduced k cols)) + (lda rows) + (s (make-array (min rows cols) :element-type ',(or real-type type))) + (ldu rows) + (ldvt vt-rows) + (work1 (make-array (max 1 lwork) :element-type ',type)) + (work nil) + ,@(when real-type + `((rwork (make-array (* 5 (min rows cols)) :element-type ',real-type))))) + (let ((u (make-array (* ldu rows) :element-type ',type)) + (vt (make-array (* ldvt cols) :element-type ',type))) + ;; run it once as a workspace query + (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt + work1 lwork ,@(when real-type `(rwork)) info) + (setf lwork (round (realpart (aref work1 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,svd-function jobu jobvt rows cols a lda s u ldu vt ldvt + work lwork ,@(when real-type `(rwork)) info) + (let ((smat (make-array (* u-cols vt-rows) :element-type ',(or real-type type)))) + (dotimes (i k) + (setf (aref smat (matrix-column-major-index i i u-cols vt-rows)) + (aref s i))) + (values (from-array u (list rows u-cols) :layout :column-major) + (from-array smat (list u-cols vt-rows) :layout :column-major) + (from-array vt (list vt-rows cols) :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) + `(progn + (defmethod eig ((m ,class)) + (lapack-eig m)) + + (defmethod lapack-eig ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m))) + (let ((rows (nrows m)) + (cols (ncols m)) + (a-tensor (deep-copy-tensor m))) + (when (eql :row-major (layout m)) (transpose! a-tensor)) + (let ((jobvl "N") + (jobvr "V") + (a (storage a-tensor)) + ,@(if real-type + `((w (make-array rows :element-type ',type))) + `((wr (make-array rows :element-type ',type)) + (wi (make-array rows :element-type ',type)))) + (vl (make-array rows :element-type ',type)) + (vr (make-array (* rows rows) :element-type ',type)) + (lwork -1) + (info 0) + ,@(when real-type + `((rwork (make-array (* 2 rows) :element-type ',real-type))))) + (let ((work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,eig-function jobvl jobvr rows a rows ,@(if real-type `(w) `(wr wi)) + vl 1 vr rows work lwork ,@(when real-type `(rwork)) info) + (setf lwork (truncate (realpart (row-major-aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; 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) :layout :column-major))))))))) + +(defun generate-lapack-hermitian-eig-for-type (class type eig-function real-type) + `(progn + (defmethod hermitian-eig ((m ,class)) + (lapack-hermitian-eig m)) + + (defmethod lapack-hermitian-eig ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m)) + (assertion (hermitian-matrix-p m))) + (let ((rows (nrows m)) + (a-tensor (deep-copy-tensor m))) + (when (eql :row-major (layout m)) (transpose! a-tensor)) + (let ((jobz "V") + (uplo "U") + (n rows) + (a (storage a-tensor)) + (lda rows) + (w (make-array rows :element-type ',real-type)) + (work (make-array 1 :element-type ',type)) + (lwork -1) + (rwork (make-array (- (* 3 rows) 2) :element-type ',real-type)) + (info 0)) + ;; run it once as a workspace query + (,eig-function jobz uplo n a lda w work lwork rwork info) + (setf lwork (truncate (realpart (row-major-aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,eig-function jobz uplo n a lda w work lwork rwork info) + (values (coerce w 'list) a-tensor))))))) + +;; TODO: implement row-major checks in these functions +(defun generate-lapack-ql-qr-rq-lq-for-type (class type + ql-function qr-function rq-function lq-function + ql-q-function qr-q-function rq-q-function lq-q-function) + `(progn + (defmethod qr ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m))) + (let ((rows (nrows m)) + (cols (ncols m))) + (multiple-value-bind (a tau) (lapack-qr m) + (let* ((r (upper-triangular a cols)) + (q (lapack-qr-q a tau))) + ;; change signs if diagonal elements of r are negative + (dotimes (j cols) + (let ((diag-elt (tref r j j))) + (assert (zerop (imagpart diag-elt)) + () "Diagonal element R_~D~D=~A is not real" j j diag-elt) + (setf diag-elt (realpart diag-elt)) + (when (minusp diag-elt) + (dotimes (i rows) + (when (<= j i (1- cols)) + (setf (tref r j i) (- (tref r j i)))) + (setf (tref q i j) (- (tref q i j))))))) + (values q r)))))) + + (defmethod ql ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m))) + (let ((rows (nrows m)) + (cols (ncols m))) + (multiple-value-bind (a tau) (lapack-ql m) + (let* ((l (lower-triangular a cols)) + (q (lapack-ql-q a tau))) + ;; change signs if diagonal elements of l are negative + (dotimes (j cols) + (let ((diag-elt (tref l j j))) + (assert (zerop (imagpart diag-elt)) + () "Diagonal element L_~D~D=~A is not real" j j diag-elt) + (setf diag-elt (realpart diag-elt)) + (when (minusp diag-elt) + (dotimes (i rows) + (when (<= i j) + (setf (tref l j i) (- (tref l j i)))) + (setf (tref q i j) (- (tref q i j))))))) + (values q l)))))) + + (defmethod rq ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m))) + (let ((rows (nrows m)) + (cols (ncols m))) + (multiple-value-bind (a tau) (lapack-rq m) + (let* ((r (upper-triangular a rows)) + (q (lapack-rq-q a tau))) + ;; change signs if diagonal elements of r are negative + (dotimes (i rows) + (let ((diag-elt (tref r i i))) + (assert (zerop (imagpart diag-elt)) + () "Diagonal element R_~D~D=~A is not real" i i diag-elt) + (setf diag-elt (realpart diag-elt)) + (when (minusp diag-elt) + (dotimes (j cols) + (when (<= j i) + (setf (tref r j i) (- (tref r j i)))) + (setf (tref q i j) (- (tref q i j))))))) + (values q r)))))) + + (defmethod lq ((m ,class)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p m))) + (let ((rows (nrows m)) + (cols (ncols m))) + (multiple-value-bind (a tau) (lapack-lq m) + (let* ((l (lower-triangular a rows)) + (q (lapack-lq-q a tau))) + ;; change signs if diagonal elements of l are negative + (dotimes (i rows) + (let ((diag-elt (tref l i i))) + (assert (zerop (imagpart diag-elt)) + () "Diagonal element L_~D~D=~A is not real" i i diag-elt) + (setf diag-elt (realpart diag-elt)) + (when (minusp diag-elt) + (dotimes (j cols) + (when (<= i j (1- rows)) + (setf (tref l j i) (- (tref l j i)))) + (setf (tref q i j) (- (tref q i j))))))) + (values q l)))))) + + (defmethod lapack-qr ((m ,class)) + (let* ((rows (nrows m)) + (cols (ncols m)) + (a-tensor (deep-copy-tensor m)) + (a (storage a-tensor)) + (lwork -1) + (info 0)) + (when (eql :row-major (layout m)) + (transpose! a-tensor)) + (let ((lda rows) + (tau (make-array (min rows cols) :element-type ',type)) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,qr-function rows cols a lda tau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,qr-function rows cols a lda tau work lwork info) + (values a-tensor + (from-array tau (list (min rows cols)) + :type ',type + :layout :column-major))))) + + (defmethod lapack-ql ((m ,class)) + (let* ((rows (nrows m)) + (cols (ncols m)) + (a-tensor (deep-copy-tensor m)) + (a (storage a-tensor)) + (lwork -1) + (info 0)) + (when (eql :row-major (layout m)) + (transpose! a-tensor)) + (let ((lda rows) + (tau (make-array (min rows cols) :element-type ',type)) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,ql-function rows cols a lda tau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,ql-function rows cols a lda tau work lwork info) + (values a-tensor + (from-array tau (list (min rows cols)) + :type ',type + :layout :column-major))))) + + (defmethod lapack-rq ((m ,class)) + (let* ((rows (nrows m)) + (cols (ncols m)) + (a-tensor (deep-copy-tensor m)) + (a (storage a-tensor)) + (lwork -1) + (info 0)) + (when (eql :row-major (layout m)) + (transpose! a-tensor)) + (let ((lda rows) + (tau (make-array (min rows cols) :element-type ',type)) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,rq-function rows cols a lda tau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,rq-function rows cols a lda tau work lwork info) + (values a-tensor + (from-array tau (list (min rows cols)) + :type ',type + :layout :column-major))))) + + (defmethod lapack-lq ((m ,class)) + (let* ((rows (nrows m)) + (cols (ncols m)) + (a-tensor (deep-copy-tensor m)) + (a (storage a-tensor)) + (lwork -1) + (info 0)) + (when (eql :row-major (layout m)) + (transpose! a-tensor)) + (let ((lda rows) + (tau (make-array (min rows cols) :element-type ',type)) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,lq-function rows cols a lda tau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,lq-function rows cols a lda tau work lwork info) + (values a-tensor + (from-array tau (list (min rows cols)) + :type ',type + :layout :column-major))))) + + (defmethod lapack-qr-q ((m ,class) tau) + (let ((m (nrows m)) + (n (ncols m)) + (a (storage m)) + (k (size tau)) + (atau (storage tau)) + (lwork -1) + (info 0)) + ;; n replaced with rank (k) to fulfil req (m >= n > 0) + (let ((lda m) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,qr-q-function m n k a lda atau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,qr-q-function m n k a lda atau work lwork info) + (from-array a (list m k) :layout :column-major)))) + + (defmethod lapack-ql-q ((m ,class) tau) + (let ((m (nrows m)) + (n (ncols m)) + (a (storage m)) + (k (size tau)) + (atau (storage tau)) + (lwork -1) + (info 0)) + (let ((lda m) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,ql-q-function m n k a lda atau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,ql-q-function m n k a lda atau work lwork info) + (from-array a (list m n) :layout :column-major)))) + + (defmethod lapack-rq-q ((m ,class) tau) + (let ((m (nrows m)) + (n (ncols m)) + (a (storage m)) + (k (size tau)) + (atau (storage tau)) + (lwork -1) + (info 0)) + (let ((lda m) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,rq-q-function m n k a lda atau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,rq-q-function m n k a lda atau work lwork info) + (from-array a (list m n) :layout :column-major)))) + + (defmethod lapack-lq-q ((m ,class) tau) + (let ((m (nrows m)) + (n (ncols m)) + (a (storage m)) + (k (size tau)) + (atau (storage tau)) + (lwork -1) + (info 0)) + (let ((lda m) + (work (make-array (max 1 lwork) :element-type ',type))) + ;; run it once as a workspace query + (,lq-q-function m n k a lda atau work lwork info) + (setf lwork (round (realpart (aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type ',type)) + ;; run it again with optimal workspace size + (,lq-q-function m n k a lda atau work lwork info) + (from-array a (list m n) :layout :column-major)))))) diff --git a/src/high-level/matrix.lisp b/src/high-level/matrix.lisp new file mode 100644 index 00000000..2511214a --- /dev/null +++ b/src/high-level/matrix.lisp @@ -0,0 +1,555 @@ +;;;; matrix.lisp +;;;; +;;;; Author: Cole Scott +;;;; Robert Smith + +(in-package #:magicl) + +(deftype matrix-storage (&optional type) + `(simple-array ,type (*))) + +(defstruct (matrix (:include abstract-tensor) + (:constructor nil) + (:copier nil)) + (nrows 0 :type alexandria:non-negative-fixnum) + (ncols 0 :type alexandria:non-negative-fixnum) + (size 0 :type alexandria:non-negative-fixnum :read-only t) + (layout :column-major :type (member :row-major :column-major))) + +(defmethod nrows ((m matrix)) + (matrix-nrows m)) + +(defmethod ncols ((m matrix)) + (matrix-ncols m)) + +(defmethod size ((m matrix)) + (matrix-size m)) + +(defmethod layout ((m matrix)) + (matrix-layout m)) + +;;; Specfic matrix classes +(defmacro defmatrix (name type tensor-class) + "Define a new matrix subclass with the specified NAME and element TYPE, +compatible with TENSOR-CLASS, as well as the abstract-tensor methods +required not specified by the generic MATRIX class (MAKE-TENSOR, +ELEMENT-TYPE, CAST, COPY-TENSOR, DEEP-COPY-TENSOR, TREF, SETF TREF)" + (let ((constructor-sym (intern (format nil "MAKE-~:@(~A~)" name))) + (copy-sym (intern (format nil "COPY-~:@(~A~)" name))) + (storage-sym (intern (format nil "~:@(~A~)-STORAGE" name)))) + `(progn + (defstruct (,name (:include matrix) + (:constructor ,constructor-sym + (nrows ncols size layout storage)) + (:copier ,copy-sym)) + (storage nil :type (matrix-storage ,type))) + #+sbcl (declaim (sb-ext:freeze-type ,name)) + + (defmethod storage ((m ,name)) + (,storage-sym m)) + + (defmethod element-type ((m ,name)) + (declare (ignore m)) + ',type) + + (defmethod make-storage ((class (eql ',name)) size initial-element) + (apply #'make-array + size + :element-type ',type + (if initial-element + (list :initial-element (coerce initial-element ',type)) + nil))) + + (defmethod make-tensor ((class (eql ',name)) shape &key initial-element layout storage) + (declare (type list shape) + (optimize (speed 3) (safety 0))) ;; This is probably not what you want... + (policy-cond:with-expectations (> speed safety) + ((type shape shape) + (assertion (cl:= 2 (length shape)))) + (let ((rows (first shape)) + (cols (second shape))) + (declare (type fixnum rows cols)) + (let ((size (the fixnum (* rows cols)))) + (funcall #',constructor-sym + rows + cols + size + (or layout :column-major) + (or + storage + (make-storage ',name size initial-element))))))) + + (defmethod cast ((tensor ,name) (class (eql ',name))) + (declare (ignore class)) + tensor) + (defmethod cast :before ((tensor ,tensor-class) (class (eql ',name))) + (declare (ignore class)) + (assert (cl:= 2 (order tensor)) + () + "Cannot change non-2 dimensional tensor to matrix.")) + (defmethod cast ((tensor ,tensor-class) (class (eql ',name))) + (declare (ignore class)) + (make-tensor ',name (shape tensor) + :storage (storage tensor) + :layout (layout tensor))) + (defmethod cast ((tensor ,name) (class (eql ',tensor-class))) + (declare (ignore class)) + (make-tensor ',tensor-class (shape tensor) + :storage (storage tensor) + :layout (layout tensor))) + + ;; TODO: This does not allow for args. Make this allow for args. + (defmethod copy-tensor ((m ,name) &rest args) + (declare (ignore args)) + (let ((new-m (,copy-sym m))) + (setf (,storage-sym new-m) + (make-array (matrix-size m) :element-type (element-type m))) + new-m)) + + (defmethod deep-copy-tensor ((m ,name) &rest args) + (declare (ignore args)) + (let ((new-m (,copy-sym m))) + (setf (,storage-sym new-m) + (copy-seq (,storage-sym m))) + new-m)) + + (defmethod tref ((matrix ,name) &rest pos) + (declare (dynamic-extent pos)) + (let ((numrows (matrix-nrows matrix)) + (numcols (matrix-ncols matrix))) + (declare (type fixnum numrows numcols)) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-matrix-index-p pos numrows numcols))) + (let ((row (first pos)) + (col (second pos))) + (declare (type fixnum row col)) + (let ((index (ecase (matrix-layout matrix) + (:row-major (+ col (the fixnum (* row numcols)))) + (:column-major (+ row (the fixnum (* col numrows))))))) + (declare (type alexandria:array-index index)) + (aref (,storage-sym matrix) index)))))) + + (defmethod (setf tref) (new-value (matrix ,name) &rest pos) + (declare (dynamic-extent pos)) + (let ((numrows (matrix-nrows matrix)) + (numcols (matrix-ncols matrix))) + (declare (type alexandria:non-negative-fixnum numrows numcols)) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-matrix-index-p pos numrows numcols))) + (let ((row (first pos)) + (col (second pos))) + (declare (type alexandria:non-negative-fixnum row col)) + (let ((index (ecase (matrix-layout matrix) + (:row-major (+ col (the fixnum (* row numcols)))) + (:column-major (+ row (the fixnum (* col numrows))))))) + (declare (type alexandria:array-index index)) + (setf (aref (,storage-sym matrix) index) + (coerce new-value ',type))))))) + + (defmethod into! ((function function) (matrix ,name)) + (let ((index 0) ;; TODO: make if less iffy + (storage (storage matrix))) + (if (eql :row-major (layout matrix)) + (loop :for j :below (nrows matrix) :do + (loop :for i :below (ncols matrix) :do + (setf (aref storage index) + (coerce (funcall function j i) ',type)) + (incf index))) + (loop :for j :below (ncols matrix) :do + (loop :for i :below (nrows matrix) :do + (setf (aref storage index) + (coerce (funcall function i j) ',type)) + (incf index)))) + matrix))))) + +(defun pprint-matrix (stream matrix &optional colon-p at-sign-p) + "Pretty-print a matrix MATRIX to the stream STREAM." + (declare (ignore colon-p) + (ignore at-sign-p)) + (flet ((print-real (x) + (format stream "~6,3f" x)) + (print-complex (z) + (format stream "~6,3f ~:[+~;-~]~6,3fj" + (realpart z) + (minusp (imagpart z)) + (abs (imagpart z)))) + (print-int (x) + (format stream "~3d" x))) + (let* ((rows (nrows matrix)) + (cols (ncols matrix)) + (type (element-type matrix)) + (print-entry + (cond + ((subtypep type 'complex) #'print-complex) + ((subtypep type 'integer) #'print-int) + (t #'print-real)))) + (pprint-logical-block (stream nil) + (print-unreadable-object (matrix stream :type t) + (format stream "(~Dx~D):" rows cols) + (dotimes (r rows) + (pprint-newline :mandatory stream) + (dotimes (c cols) + (funcall print-entry (tref matrix r c)) + (unless (cl:= c (1- cols)) + (write-string " " stream))))))))) + +(set-pprint-dispatch 'matrix 'pprint-matrix) + +(defun square-matrix-p (matrix) + "Whether the MATRIX is square" + (cl:= (nrows matrix) (ncols matrix))) + +(defun identity-matrix-p (matrix &optional (epsilon *double-comparison-threshold*)) + "Whether MATRIX is an idenity matrix" + (unless (square-matrix-p matrix) (return-from identity-matrix-p nil)) + (map-indexes (shape matrix) + (lambda (r c) + (unless (<= (abs (- (tref matrix r c) + (if (cl:= r c) + 1 0))) + epsilon) + (return-from identity-matrix-p nil)))) + t) + +(defun unitary-matrix-p (matrix &optional (epsilon *double-comparison-threshold*)) + "Whether MATRIX is a unitary matrix" + (identity-matrix-p (@ matrix (conjugate-transpose matrix)) epsilon)) + +(defun hermitian-matrix-p (matrix &optional (epsilon *double-comparison-threshold*)) + "Whether MATRIX is a unitary matrix" + (= matrix (conjugate-transpose matrix) epsilon)) + +(defmacro assert-square-matrix (&rest matrices) + `(progn + ,@(loop :for matrix in matrices + :collect `(assert (square-matrix-p ,matrix) + () + ,"The shape of ~a is ~a, which is not a square" ,(symbol-name matrix) (shape ,matrix))))) + +;;; Required abstract-tensor methods + +(defmethod order ((m matrix)) + (declare (ignore m)) + 2) + +(defmethod shape ((m matrix)) + (list (matrix-nrows m) (matrix-ncols m))) + +(defmethod (setf shape) (new-value (m matrix)) + (policy-cond:with-expectations (> speed safety) + ((type shape new-value) + (assertion (cl:= 2 (length new-value)))) + (setf (matrix-nrows m) (first new-value) + (matrix-ncols m) (second new-value)))) + +;; Specific constructors + +(defgeneric random-unitary (shape &key type) + (:documentation "Generate a uniformly random element of U(n).") + (:method (shape &key (type 'double-float)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-shape-p shape))) + (multiple-value-bind (q r) (qr (rand shape :type type :distribution #'alexandria:gaussian-random)) + (let ((d (diag r))) + (setf d (cl:map 'list (lambda (di) (/ di (sqrt (* di (conjugate di))))) d)) + (@ q (funcall #'from-diag d))))))) + +;; TODO: This should be generic to abstract-tensor +(defgeneric ptr-ref (m base i j) + (:documentation + "Accessor method for the pointer to the element in the I-th row and J-th column of a matrix M, assuming zero indexing.") + (:method ((m matrix) base i j) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-index-p (list i j) (shape m)))) + (let ((type (element-type m))) + ;; TODO: compensate for layout + (let ((idx (apply #'matrix-column-major-index i j (shape m)))) + (cond + ((subtypep type 'single-float) (cffi:mem-aptr base :float idx)) + ((subtypep type 'double-float) (cffi:mem-aptr base :double idx)) + ((subtypep type '(complex single-float)) (cffi:mem-aptr base :float (* 2 idx))) + ((subtypep type '(complex double-float)) (cffi:mem-aptr base :double (* 2 idx))) + (t (error "Incompatible element type ~a." type)))))))) + +(defgeneric row (matrix index) + (:documentation "Get row vector from a matrix") + (:method ((m matrix) index) + (check-type index alexandria:non-negative-fixnum) + (slice m + (list index 0) + (list (1+ index) (ncols m))))) + +(defgeneric column (matrix index) + (:documentation "Get column vector from a matrix") + (:method ((m matrix) index) + (slice m + (list 0 index) + (list (nrows m) (1+ index))))) + +(defgeneric mult (a b &key target alpha beta transa transb) + (:documentation "Multiply matrix a by matrix b, storing in target or creating a new matrix if target is not specified. +Target cannot be the same as a or b.")) + +(defgeneric @ (matrix &rest matrices) + (:documentation "Multiplication of matrices") + (:method (matrix &rest matrices) + (reduce #'mult matrices + :initial-value matrix))) + +;;; Generic matrix methods + +(defgeneric direct-sum (a b) + (:method ((a matrix) (b matrix)) + "Compute the direct sum of A and B" + (let* ((arows (nrows a)) + (acols (ncols a)) + (brows (nrows b)) + (bcols (ncols b)) + (rrows (+ arows brows)) + (rcols (+ acols bcols)) + (result (magicl:empty (list rrows rcols) + :type '(complex double-float)))) + (loop :for r :below arows :do + (loop :for c :below acols :do + (setf (tref result r c) (tref a r c)))) + (loop :for r :from arows :below rrows :do + (loop :for c :from acols :below rcols :do + (setf (tref result r c) (tref b (- r arows) (- c acols))))) + result))) + +(defgeneric kron (a b &rest rest) + (:documentation "Compute the Kronecker product of A and B") + (:method (a b &rest rest) + (let ((ma (nrows a)) + (mb (nrows b)) + (na (ncols a)) + (nb (ncols b))) + (flet ((calc-i-j (i j) (* (tref a (floor i mb) (floor j nb)) + (tref b (mod i mb) (mod j nb))))) + (reduce #'kron rest :initial-value (into! #'calc-i-j + (empty (list (* ma mb) (* na nb)) + :type '(complex double-float)))))))) + +(defgeneric transpose! (matrix &key fast) + (:documentation "Transpose MATRIX, replacing the elements of MATRIX, optionally performing a faster change of layout if FAST is specified") + (:method ((matrix matrix) &key fast) + "Transpose a matrix by copying values. +If FAST is t then just change layout. Fast can cause problems when you want to multiply specifying transpose." + (if fast + (progn (rotatef (matrix-ncols matrix) (matrix-nrows matrix)) + (setf (matrix-layout matrix) (ecase (matrix-layout matrix) + (:row-major :column-major) + (:column-major :row-major)))) + (let ((index-function + (ecase (matrix-layout matrix) + (:row-major #'matrix-row-major-index) + (:column-major #'matrix-column-major-index))) + (shape (shape matrix))) + (loop :for row :below (matrix-nrows matrix) + :do (loop :for col :from row :below (matrix-ncols matrix) + :do (rotatef + (aref (storage matrix) (apply index-function row col shape)) + (aref (storage matrix) (apply index-function col row shape))))) + (rotatef (matrix-ncols matrix) (matrix-nrows matrix)))) + matrix)) + +(defgeneric transpose (matrix) + (:documentation "Create a new matrix containing the transpose of MATRIX") + (:method ((matrix matrix)) + "Transpose a matrix by copying values. +If fast is t then just change layout. Fast can cause problems when you want to multiply specifying transpose." + (let ((new-matrix (copy-tensor matrix))) + (setf (matrix-ncols new-matrix) (matrix-nrows matrix) + (matrix-nrows new-matrix) (matrix-ncols matrix)) + (let ((index-function + (ecase (layout matrix) + (:row-major #'matrix-row-major-index) + (:column-major #'matrix-column-major-index))) + (shape (shape matrix))) + (loop :for row :below (matrix-nrows matrix) + :do (loop :for col :below (matrix-ncols matrix) + :do (let ((index1 (apply index-function row col shape)) + (index2 (apply index-function col row (shape new-matrix)))) + (setf (aref (storage new-matrix) index2) (aref (storage matrix) index1)))))) + new-matrix))) + +;; TODO: allow setf on matrix diag +(defgeneric diag (matrix) + (:documentation "Get a list of the diagonal elements of MATRIX") + (:method ((matrix matrix)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p matrix))) + (let ((rows (nrows matrix))) + (loop :for i :below rows + :collect (tref matrix i i)))))) + +(defgeneric trace (matrix) + (:documentation "Get the trace of MATRIX (sum of diagonals)") + (:method ((matrix matrix)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p matrix))) + (loop :for i :below (nrows matrix) + :sum (tref matrix i i))))) + +(defgeneric det (matrix) + (:documentation "Compute the determinant of a square matrix MATRIX") + (:method ((matrix matrix)) + (policy-cond:with-expectations (> speed safety) + ((assertion (square-matrix-p matrix))) + (let ((d 1)) + (multiple-value-bind (a ipiv) (lu matrix) + (dotimes (i (nrows matrix)) + (setq d (* d (tref a i i)))) + (dotimes (i (size ipiv)) + (unless (cl:= (1+ i) (tref ipiv i)) + (setq d (- d)))) + d))))) + +(defgeneric upper-triangular (matrix &optional order) + (:documentation "Get the upper triangular portion of the matrix") + (:method ((matrix matrix) &optional (order (ncols matrix))) + (assert-square-matrix matrix) + (let ((m (nrows matrix)) + (n (ncols matrix))) + (policy-cond:with-expectations (> speed safety) + ((assertion (<= order (max (nrows matrix) (ncols matrix))))) + (let ((target (empty (list order order) :layout (layout matrix) :type (element-type matrix)))) + (if (> m n) + (loop for i from 0 to (1- order) + do (loop for j from (max 0 (+ (- n order) i)) to (1- n) + do (setf (tref target i (+ j (- order n))) (tref matrix i j)))) + (loop for j from (- n order) to (1- n) + do (loop for i from 0 to (min (+ (- order n) j) (1- m)) + do (setf (tref target i (- j (- n order))) (tref matrix i j))))) + target))))) +;;; Synonym for upper-triangular +(setf (fdefinition 'triu) #'upper-triangular) + +(defgeneric lower-triangular (matrix &optional order) + (:documentation "Get the lower triangular portion of the matrix") + (:method ((matrix matrix) &optional (order (ncols matrix))) + (assert-square-matrix matrix) + (let ((m (nrows matrix)) + (n (ncols matrix))) + (policy-cond:with-expectations (> speed safety) + ((assertion (<= order (max (nrows matrix) (ncols matrix))))) + (let ((target (empty (list order order) :layout (layout matrix) :type (element-type matrix)))) + (if (> m n) + (loop for i from (- m order) to (1- m) + do (loop for j from 0 to (min (+ (- order m) i) (1- n)) + do (setf (tref target (- i (- m order)) j) (tref matrix i j)))) + (loop for j from 0 to (1- order) + do (loop for i from (max 0 (+ (- m order) j)) to (1- m) + do (setf (tref target (+ i (- order m)) j) (tref matrix i j))))) + target))))) +;;; Synonym for lower-triangular +(setf (fdefinition 'tril) #'lower-triangular) + +(defgeneric conjugate-transpose (matrix) + (:documentation "Compute the conjugate transpose of a matrix") + (:method ((matrix matrix)) + (map #'conjugate (transpose matrix)))) +(setf (fdefinition 'dagger) #'conjugate-transpose) + +(defgeneric conjugate-transpose! (matrix) + (:documentation "Compute the conjugate transpose of a matrix, replacing the elements") + (:method ((matrix matrix)) + (map! #'conjugate (transpose! matrix)))) +(setf (fdefinition 'dagger!) #'conjugate-transpose!) + +(defgeneric eig (matrix) + (:documentation "Find the (right) eigenvectors and corresponding eigenvalues of a square matrix M. Returns a list and a tensor (EIGENVALUES, EIGENVECTORS)") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "EIG is not defined for the generic matrix type."))) + +(defgeneric lu (matrix) + (:documentation "Get the LU decomposition of MATRIX. Returns two tensors (LU, IPIV)") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "LU is not defined for the generic matrix type."))) + +;; TODO: Make this one generic and move to lapack-macros +;; This is being blocked by the ZUNCSD shenanigans +(defgeneric csd (matrix p q) + (:documentation "Find the Cosine-Sine Decomposition of a matrix X given that it is to be partitioned with upper left block of dimension P-by-Q. Returns the CSD elements (VALUES U SIGMA VT) such that X=U*SIGMA*VT.") + (:method ((matrix matrix) 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))))) + +(defgeneric inv (matrix) + (:documentation "Get the inverse of the matrix") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "INV is not defined for the generic matrix type."))) + +(defgeneric svd (matrix &key reduced) + (:documentation "Find the SVD of a matrix M. Return (VALUES U SIGMA Vt) where M = U*SIGMA*Vt") + (:method ((matrix matrix) &key reduced) + (declare (ignore matrix reduced)) + (error "SVD is not defined for the generic matrix type."))) + +(defgeneric qr (matrix) + (:documentation "Finds the QL factorization of the matrix M. Returns two tensors (Q, R) +NOTE: Only square matrices supported") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "QR is not defined for the generic matrix type."))) + +(defgeneric ql (matrix) + (:documentation "Finds the QL factorization of the matrix M. Returns two tensors (Q, L) +NOTE: Only square matrices supported") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "QL is not defined for the generic matrix type."))) + +(defgeneric rq (matrix) + (:documentation "Finds the RQ factorization of the matrix M. Returns two tensors (R, Q) +NOTE: Only square matrices supported") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "RQ is not defined for the generic matrix type."))) + +(defgeneric lq (matrix) + (:documentation "Finds the LQ factorization of the matrix M. Returns two tensors (L, Q) +NOTE: Only square matrices supported") + (:method ((matrix matrix)) + (declare (ignore matrix)) + (error "LQ is not defined for the generic matrix type."))) diff --git a/src/high-level/polynomial-solver.lisp b/src/high-level/polynomial-solver.lisp index 675cc60b..e0576588 100644 --- a/src/high-level/polynomial-solver.lisp +++ b/src/high-level/polynomial-solver.lisp @@ -49,11 +49,11 @@ It is assumed that POLYNOMIAL is well-formed in the sense that its leading coeff (loop :with polynomial := (polynomial-normalize polynomial) :with coefficients := (polynomial-coefficients polynomial) :with n := (1- (length coefficients)) - :with matrix := (make-zero-matrix n n) + :with matrix := (zeros (list n n) :type '(complex double-float)) :for i :below n :for coefficient :of-type (complex double-float) := (aref coefficients i) - :when (plusp i) :do (setf (ref matrix i (1- i)) 1.0d0) - :do (setf (ref matrix i (1- n)) (- coefficient)) + :when (plusp i) :do (setf (tref matrix i (1- i)) 1.0d0) + :do (setf (tref matrix i (1- n)) (- coefficient)) :finally (return matrix)))) (nth-value 0 (eig (make-companion-matrix polynomial))))) @@ -68,7 +68,7 @@ It is assumed that POLYNOMIAL is well-formed in the sense that its leading coeff (c (aref coefficients 0) (aref coefficients i)) (x (complex 1.0d0) (* x value)) (y (* x c) (+ y (* x c)))) - ((= i (1- (length coefficients))) y) + ((cl:= i (1- (length coefficients))) y) (declare (type fixnum i) (type (complex double-float) x y))))) diff --git a/src/high-level/random.lisp b/src/high-level/random.lisp deleted file mode 100644 index 5780c7ba..00000000 --- a/src/high-level/random.lisp +++ /dev/null @@ -1,30 +0,0 @@ -;;;; random.lisp -;;;; -;;;; Author: Robert Smith - -(in-package #:magicl) - -;;; Some routines to create random matrices. - -(defun random-matrix (rows cols) - "Create a Z matrix of size ROWS x COLS with uniformly random complex entries in [0,1) + [0,1)i." - (tabulate rows cols - (lambda (i j) - (declare (ignore i j)) - (complex (random 1.0d0) (random 1.0d0))))) - -(defun random-gaussian-matrix (rows cols) - "Create a Z matrix of size ROWS x COLS with normally distributed random complex entries." - (tabulate rows cols - (lambda (i j) - (declare (ignore i j)) - (complex (alexandria:gaussian-random) - (alexandria:gaussian-random))))) - -(defun random-unitary (n) - "Generate a uniformly random element of U(n)." - (declare (optimize (speed 0) (debug 3))) - (multiple-value-bind (q r) (qr (random-gaussian-matrix n n)) - (let ((d (matrix-diagonal r))) - (map-into d (lambda (di) (/ di (sqrt (* di (conjugate di))))) d) - (multiply-complex-matrices q (funcall #'diag n n d))))) diff --git a/src/high-level/shape.lisp b/src/high-level/shape.lisp new file mode 100644 index 00000000..309b7485 --- /dev/null +++ b/src/high-level/shape.lisp @@ -0,0 +1,61 @@ +;;;; shapes.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +;;; Shapes + +;; Predicates + +(declaim (inline valid-shape-p)) +(defun valid-shape-p (shape) + (and (typep shape 'list) + (plusp (length shape)) + (cl:every (lambda (x) (typep x 'alexandria:positive-fixnum)) shape))) + +(declaim (inline square-shape-p)) +(defun square-shape-p (shape) + (and (valid-shape-p shape) + (apply #'cl:= shape))) + +(declaim (inline valid-index-p)) +(defun valid-index-p (index &optional shape) + (declare (notinline valid-index-p)) + (if (null shape) + (and (typep index 'list) + (plusp (length index)) + (cl:every (lambda (x) (typep x 'alexandria:non-negative-fixnum)) index)) + (and (valid-index-p index) + (cl:= (length index) (length shape)) + (cl:every #'< index shape)))) + +(declaim (inline valid-matrix-index-p)) +(defun valid-matrix-index-p (index &optional nrows ncols) + (if (or (null nrows) (null ncols)) + (and (typep index 'list) + (cl:= 2 (length index)) + (cl:every (lambda (x) (typep x 'alexandria:non-negative-fixnum)) index)) + (and (typep index 'list) + (cl:= 2 (length index)) + (cl:every (lambda (x) (typep x 'alexandria:non-negative-fixnum)) index) + (< (first index) nrows) + (< (second index) ncols)))) + +;; Types +(deftype shape (&optional order) + `(satisfies valid-shape-p)) + +(deftype index () + `(satisfies valid-index-p)) + +;; Assertions +(defmacro assert-square-shape (&rest shapes) + `(progn + ,@(loop :for shape in shapes + :collect `(assert (square-shape-p ,shape) + () + "The value of ~a is ~a, which is not a square SHAPE" ,(symbol-name shape) ,shape)))) + +(defun fixnum-to-shape (num &optional (order 2)) + (make-list order :initial-element num)) diff --git a/src/high-level/specialize-constructor.lisp b/src/high-level/specialize-constructor.lisp new file mode 100644 index 00000000..446be179 --- /dev/null +++ b/src/high-level/specialize-constructor.lisp @@ -0,0 +1,58 @@ +;;;; specialize-constructor.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +;; NOTE: I hate this +(defun infer-tensor-type (type shape val) + (declare (type list shape) + (optimize (speed 3) (safety 0))) + (if type + (cond + ((cl:= 1 (length shape)) + (cond + ((subtypep type 'single-float) 'vector/single-float) + ((subtypep type 'double-float) 'vector/double-float) + ((subtypep type '(complex single-float)) 'vector/complex-single-float) + ((subtypep type '(complex double-float)) 'vector/complex-double-float) + ((subtypep type '(signed-byte 32)) 'vector/int32) + (t (error "no compatible tensor constructor for type ~a" type)))) + ((cl:= 2 (length shape)) + (cond + ((subtypep type 'single-float) 'matrix/single-float) + ((subtypep type 'double-float) 'matrix/double-float) + ((subtypep type '(complex single-float)) 'matrix/complex-single-float) + ((subtypep type '(complex double-float)) 'matrix/complex-double-float) + ((subtypep type '(signed-byte 32)) 'matrix/int32) + (t (error "no compatible tensor constructor for type ~a" type)))) + (t + (cond + ((subtypep type 'single-float) 'tensor/single-float) + ((subtypep type 'double-float) 'tensor/double-float) + ((subtypep type '(complex single-float)) 'tensor/complex-single-float) + ((subtypep type '(complex double-float)) 'tensor/complex-double-float) + ((subtypep type '(signed-byte 32)) 'tensor/int32) + (t (error "no compatible tensor constructor for type ~a" type))))) + (cond + ((cl:= 1 (length shape)) + (etypecase val + (single-float 'vector/single-float) + (double-float 'vector/double-float) + ((complex single-float) 'vector/complex-single-float) + ((complex double-float) 'vector/complex-double-float) + ((signed-byte 32) 'vector/int32))) + ((cl:= 2 (length shape)) + (etypecase val + (single-float 'matrix/single-float) + (double-float 'matrix/double-float) + ((complex single-float) 'matrix/complex-single-float) + ((complex double-float) 'matrix/complex-double-float) + ((signed-byte 32) 'matrix/int32))) + (t (etypecase val + (single-float 'tensor/single-float) + (double-float 'tensor/double-float) + ((complex single-float) 'tensor/complex-single-float) + ((complex double-float) 'tensor/complex-double-float) + ((signed-byte 32) 'tensor/int32)))))) + diff --git a/src/high-level/specialize-tensor.lisp b/src/high-level/specialize-tensor.lisp new file mode 100644 index 00000000..a785bb17 --- /dev/null +++ b/src/high-level/specialize-tensor.lisp @@ -0,0 +1,31 @@ +;;;; specialize-tensor.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +;; XXX: This is a little overkill for the small amount of compatible +;; tensors we are defining. Possible to change this to just +;; defining SPECIALIZE-TENSOR and GENERALIZE-TENSOR for a class +(defmacro defcompatible (function &rest tensors) + "Generate COMPATIBLE-TENSORS methods specifying compatible abstract-tensor classes for a given tensor. + +FUNCTION takes the tensor as an argument and returns a list of +compatible classes. +TENSORS specifies all the classes this function will be defined for" + `(progn + ,@(loop :for tensor :in tensors + :collect `(defmethod compatible-tensors ((tensor ,tensor)) (funcall ,function tensor))))) + +(defgeneric compatible-tensors (tensor) + (:documentation "Get a list of all compatible classes for a tensor with most-specific first")) + +(defgeneric cast (tensor class) + (:documentation "Cast a tensor to CLASS")) + +(defmethod specialize-tensor ((tensor abstract-tensor)) + (cast tensor (first (compatible-tensors tensor)))) + +(defmethod generalize-tensor ((tensor abstract-tensor)) + (cast tensor (car (last (compatible-tensors tensor))))) + diff --git a/src/high-level/tensor.lisp b/src/high-level/tensor.lisp new file mode 100644 index 00000000..9eea9ab6 --- /dev/null +++ b/src/high-level/tensor.lisp @@ -0,0 +1,132 @@ +;;;; tensor.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftype tensor-storage (&optional type) + `(simple-array ,type (*))) + +(defstruct (tensor (:include abstract-tensor) + (:constructor nil) + (:copier nil)) + (order 0 :type alexandria:non-negative-fixnum) + (shape '(0) :type list) + (size 0 :type alexandria:positive-fixnum :read-only t) + (layout :column-major :type (member :row-major :column-major))) + +(defmethod size ((a tensor)) + (tensor-size a)) + +(defmethod layout ((a tensor)) + (tensor-layout a)) + +(defmethod order ((a tensor)) + (tensor-order a)) + +(defmethod shape ((a tensor)) + (tensor-shape a)) + +(defmethod (setf shape) (new-value (a tensor)) + (reshape a new-value)) + +;;; Specfic tensor classes +(defmacro deftensor (name type) + "Define a new tensor subclass with the specified NAME and element +TYPE as well as the abstract-tensor methods required not specified by +the generic TENSOR class (MAKE-TENSOR, ELEMENT-TYPE, CAST, +COPY-TENSOR, DEEP-COPY-TENSOR, TREF, SETF TREF)" + (let ((constructor-sym (intern (format nil "MAKE-~:@(~A~)-STRUCT" name))) + (copy-sym (intern (format nil "COPY-~:@(~A~)" name))) + (storage-sym (intern (format nil "~:@(~A~)-STORAGE" name)))) + `(progn + (defstruct (,name (:include tensor) + (:constructor ,constructor-sym + (order shape size layout storage)) + (:copier ,copy-sym)) + (storage nil :type (tensor-storage ,type))) + #+sbcl (declaim (sb-ext:freeze-type ,name)) + + (defmethod storage ((m ,name)) + (,storage-sym m)) + + (defmethod element-type ((m ,name)) + (declare (ignore m)) + ',type) + + (defmethod make-tensor ((class (eql ',name)) shape &key initial-element layout storage) + (policy-cond:with-expectations (> speed safety) + ((type shape shape)) + (let ((size (reduce #'* shape))) + (funcall #',constructor-sym + (length shape) + shape + size + (or layout :column-major) + (or + storage + (apply #'make-array + size + :element-type ',type + (if initial-element + (list :initial-element (coerce initial-element ',type)) + nil))))))) + (defmethod cast ((tensor ,name) (class (eql ',name))) + (declare (ignore class)) + tensor) + + ;; TODO: This does not allow for args. Make this allow for args. + (defmethod copy-tensor ((m ,name) &rest args) + (declare (ignore args)) + (let ((new-m (,copy-sym m))) + (setf (,storage-sym new-m) + (make-array (tensor-size m) :element-type (element-type m))) + new-m)) + + (defmethod deep-copy-tensor ((m ,name) &rest args) + (declare (ignore args)) + (let ((new-m (,copy-sym m))) + (setf (,storage-sym new-m) + (copy-seq (,storage-sym m))) + new-m)) + + (defmethod tref ((tensor ,name) &rest pos) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-index-p pos (tensor-shape tensor)))) + (let ((index (case (tensor-layout tensor) + (:row-major (row-major-index pos (tensor-shape tensor))) + (:column-major (column-major-index pos (tensor-shape tensor)))))) + (aref (,storage-sym tensor) index)))) + + (defmethod (setf tref) (new-value (tensor ,name) &rest pos) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-index-p pos (tensor-shape tensor)))) + (let ((index (case (tensor-layout tensor) + (:row-major (row-major-index pos (tensor-shape tensor))) + (:column-major (column-major-index pos (tensor-shape tensor)))))) + (setf (aref (,storage-sym tensor) index) + new-value))))))) + +(defun pprint-tensor (stream tensor &optional colon-p at-sign-p) + "Pretty-print a matrix MATRIX to the stream STREAM." + (declare (ignore colon-p) + (ignore at-sign-p)) + (pprint-logical-block (stream nil) + (print-unreadable-object (tensor stream :type t) + (format stream "(~{~D~^x~})" (shape tensor))))) + +(set-pprint-dispatch 'tensor 'pprint-tensor) + +;;; Generic tensor methods + +(defgeneric reshape (tensor shape) + (:documentation "Change the shape of the tensor. +WARNING: This method acts differently depending on the layout of the tensor. Do not expect row-major to act the same as column-major.") + (:method ((tensor tensor) shape) + (policy-cond:with-expectations (> speed safety) + ((assertion (cl:= (tensor-size tensor) (reduce #'* shape)))) + (setf (tensor-shape tensor) shape) + (setf (tensor-order tensor) (length shape)) + (specialize-tensor tensor))) + (:method ((tensor abstract-tensor) shape) + (reshape (generalize-tensor tensor) shape))) diff --git a/src/high-level/types/complex-double-float.lisp b/src/high-level/types/complex-double-float.lisp new file mode 100644 index 00000000..f42b6a32 --- /dev/null +++ b/src/high-level/types/complex-double-float.lisp @@ -0,0 +1,243 @@ +;;;; complex-double-float.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftensor tensor/complex-double-float (complex double-float)) + +(defmatrix matrix/complex-double-float (complex double-float) tensor/complex-double-float) + +(defvector vector/complex-double-float (complex double-float) tensor/complex-double-float) + +(defcompatible + (lambda (tensor) + (case (order tensor) + (1 '(vector/complex-double-float + tensor/complex-double-float)) + (2 '(matrix/complex-double-float + tensor/complex-double-float)) + (t '(tensor/complex-double-float)))) + tensor/complex-double-float + matrix/complex-double-float + vector/complex-double-float) + +(defmethod dot ((vector1 vector/complex-double-float) (vector2 vector/complex-double-float)) + (assert (cl:= (size vector1) (size vector2)) + () "Vectors must have the same size. The first vector is size ~a and the second vector is size ~a." + (size vector1) (size vector2)) + (loop :for i :below (size vector1) + :sum (* (tref vector1 i) (conjugate (tref vector2 i))))) + +(defmethod = ((tensor1 tensor/complex-double-float) (tensor2 tensor/complex-double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 matrix/complex-double-float) (tensor2 matrix/complex-double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 vector/complex-double-float) (tensor2 vector/complex-double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) + +;; ZUNCSD is broken in magicl lapack bindings (Issue #72) +(COMMON-LISP:DEFUN %ZUNCSD-XPOINTERS + (JOBU1 JOBU2 JOBV1T JOBV2T TRANS SIGNS M P Q X11 LDX11 X12 + LDX12 X21 LDX21 X22 LDX22 THETA U1 LDU1 U2 LDU2 V1T LDV1T + V2T LDV2T WORK LWORK RWORK LRWORK IWORK INFO) + (CFFI:WITH-FOREIGN-OBJECTS ((M-REF71665 ':INT32) (P-REF71666 ':INT32) + (Q-REF71667 ':INT32) (LDX11-REF71669 ':INT32) + (LDX12-REF71671 ':INT32) (LDX21-REF71673 ':INT32) + (LDX22-REF71675 ':INT32) (LDU1-REF71678 ':INT32) + (LDU2-REF71680 ':INT32) (LDV1T-REF71682 ':INT32) + (LDV2T-REF71684 ':INT32) (LWORK-REF71686 ':INT32) + (LRWORK-REF71688 ':INT32) (INFO-REF71690 ':INT32)) + (COMMON-LISP:SETF (CFFI:MEM-REF M-REF71665 :INT32) M) + (COMMON-LISP:SETF (CFFI:MEM-REF P-REF71666 :INT32) P) + (COMMON-LISP:SETF (CFFI:MEM-REF Q-REF71667 :INT32) Q) + (COMMON-LISP:SETF (CFFI:MEM-REF LDX11-REF71669 :INT32) LDX11) + (COMMON-LISP:SETF (CFFI:MEM-REF LDX12-REF71671 :INT32) LDX12) + (COMMON-LISP:SETF (CFFI:MEM-REF LDX21-REF71673 :INT32) LDX21) + (COMMON-LISP:SETF (CFFI:MEM-REF LDX22-REF71675 :INT32) LDX22) + (COMMON-LISP:SETF (CFFI:MEM-REF LDU1-REF71678 :INT32) LDU1) + (COMMON-LISP:SETF (CFFI:MEM-REF LDU2-REF71680 :INT32) LDU2) + (COMMON-LISP:SETF (CFFI:MEM-REF LDV1T-REF71682 :INT32) LDV1T) + (COMMON-LISP:SETF (CFFI:MEM-REF LDV2T-REF71684 :INT32) LDV2T) + (COMMON-LISP:SETF (CFFI:MEM-REF LWORK-REF71686 :INT32) LWORK) + (COMMON-LISP:SETF (CFFI:MEM-REF LRWORK-REF71688 :INT32) LRWORK) + (COMMON-LISP:SETF (CFFI:MEM-REF INFO-REF71690 :INT32) INFO) + (MAGICL.CFFI-TYPES:WITH-ARRAY-POINTERS ((THETA-REF71676 THETA) + (U1-REF71677 U1) (U2-REF71679 U2) + (V1T-REF71681 V1T) + (V2T-REF71683 V2T) + (WORK-REF71685 WORK) + (RWORK-REF71687 RWORK) + (IWORK-REF71689 IWORK)) + (MAGICL.LAPACK-CFFI::%%ZUNCSD + JOBU1 JOBU2 JOBV1T JOBV2T TRANS SIGNS M-REF71665 P-REF71666 + Q-REF71667 X11 LDX11-REF71669 X12 LDX12-REF71671 + X21 LDX21-REF71673 X22 LDX22-REF71675 THETA-REF71676 + U1-REF71677 LDU1-REF71678 U2-REF71679 LDU2-REF71680 V1T-REF71681 + LDV1T-REF71682 V2T-REF71683 LDV2T-REF71684 WORK-REF71685 LWORK-REF71686 + RWORK-REF71687 LRWORK-REF71688 IWORK-REF71689 INFO-REF71690)))) + +(defmethod lapack-csd ((x matrix/complex-double-float) p q) + (let* ((m (nrows x)) + (layout (layout x)) + (xcopy (deep-copy-tensor x))) + (assert-square-matrix x) + (check-type p integer) + (check-type q integer) + (assert (<= 1 p (1- m)) () "P = ~D is out of range" p) + (assert (<= 1 q (1- m)) () "Q = ~D is out of range" q) + (let ((jobu1 "Y") + (jobu2 "Y") + (jobv1t "Y") + (jobv2t "Y") + (trans + (if (eql :row-major layout) + "T" + "F")) + (signs "D") + ;; leading dimension is M because full MATRIX array will be used + (ldx11 m) + (ldx12 m) + (ldx21 m) + (ldx22 m) + (r (min p (- m p) q (- m q))) + (ldu1 p) + (ldu2 (- m p)) + (ldv1t q) + (ldv2t (- m q)) + (lwork -1) + (work (make-array 1 :element-type '(complex double-float))) + (lrwork -1) + (rwork (make-array 1 :element-type 'double-float)) + (info 0)) + ;; rather than slice up matrix, use full array with pointers to head of blocks + ;; + ;; WARNING: THIS ABYSS IS WHERE DRAGONS LIVE. HERE THE GARBAGE + ;; COLLECTOR MUST BE TAMED SO WE DON'T SCREW UP THE POINTERS + ;; INTO THE LISP HEAP. + ;; + ;; HOURS WASTED HERE: 10 + (magicl.cffi-types:with-array-pointers ((xcopy-ptr (storage xcopy))) + (let ((x11 xcopy-ptr) + (x12 (ptr-ref xcopy xcopy-ptr 0 q)) + (x21 (ptr-ref xcopy xcopy-ptr p 0)) + (x22 (ptr-ref xcopy xcopy-ptr p q)) + (theta (make-array r + :element-type 'double-float)) + (u1 (make-array (* ldu1 p) + :element-type '(complex double-float))) + (u2 (make-array (* ldu2 (- m p)) + :element-type '(complex double-float))) + (v1t (make-array (* ldv1t q) + :element-type '(complex double-float))) + (v2t (make-array (* ldv2t (- m q)) + :element-type '(complex double-float))) + (iwork (make-array (- m r) + :element-type '(signed-byte 32)))) + ;; run it once as a workspace query + (%ZUNCSD-XPOINTERS jobu1 jobu2 jobv1t jobv2t + trans signs m p q + x11 ldx11 x12 ldx12 x21 ldx21 x22 ldx22 + theta u1 ldu1 u2 ldu2 v1t ldv1t v2t ldv2t + work lwork rwork lrwork iwork info) + (setf lwork (truncate (realpart (row-major-aref work 0)))) + (setf work (make-array (max 1 lwork) :element-type '(complex double-float))) + (setf lrwork (truncate (row-major-aref rwork 0))) + (setf rwork (make-array (max 1 lrwork) :element-type 'double-float)) + ;; run it again with optimal workspace size + (%ZUNCSD-XPOINTERS jobu1 jobu2 jobv1t jobv2t + trans signs m p q + x11 ldx11 x12 ldx12 x21 ldx21 x22 ldx22 + theta u1 ldu1 u2 ldu2 v1t ldv1t v2t ldv2t + work lwork rwork lrwork iwork info) + (values (from-array u1 (list p p) :layout :column-major) + (from-array u2 (list (- m p) (- m p)) :layout :column-major) + (from-array v1t (list q q) :layout :column-major) + (from-array v2t (list (- m q) (- m q)) :layout :column-major) + (coerce theta 'list))))))) + +(defmethod csd-2x2-basic ((unitary-matrix-2x2 matrix/complex-double-float) p q) + "Returns the Cosine-Sine decomposition of an equipartitioned UNITARY-MATRIX-2x2. The values of P and Q are assumed to be equal to one and ignored. See the documentation of LAPACK-CSD for details about the returned values." + ;; This function is meant to be used within LAPACK-CSD in the base case + ;; where the unitary matrix is 2x2 (i.e., it is equivalent to a ZYZ + ;; decomposition). It computes the CS decomposition in Lisp faster (more + ;; than three times) and with less memory overhead (conses less than half + ;; as much) than the corresponding LAPACK wrapper. + (declare (values matrix/complex-double-float + matrix/complex-double-float + matrix/complex-double-float + matrix/complex-double-float + list) + (ignorable p q) + (optimize (speed 3) (safety 0) (debug 0) (space 0))) + + (let* ((data (storage unitary-matrix-2x2)) + (a1 (aref data 0)) + (a2 (aref data 1)) + (a3 (aref data 2)) + (a4 (aref data 3)) + + (c (abs a1)) + (u1 (cis (phase a1))) + (s (abs a2)) + (u2 (cis (phase a2)))) + + (declare (type (simple-array (complex double-float) (4)) data) + (type (double-float -1.0d0 1.0d0) c s) + (type (complex double-float) u1 u2) + (dynamic-extent c s u1 u2)) + + (let ((v2h (conjugate (/ 1.0d0 (- (* c (conjugate u2) a4) + (* s (conjugate u1) a3))))) + (mu1 (empty '(1 1))) + (mu2 (empty '(1 1))) + (mv1h (empty '(1 1))) + (mv2h (empty '(1 1)))) + + (macrolet ((matrix-1x1-data (matrix) + `(the (simple-array (complex double-float) (1)) (storage ,matrix)))) + + (setf (aref (matrix-1x1-data mu1) 0) u1 + (aref (matrix-1x1-data mu2) 0) u2 + (aref (matrix-1x1-data mv1h) 0) #C(1.0d0 0.0d0) + (aref (matrix-1x1-data mv2h) 0) v2h)) + + (values mu1 mu2 mv1h mv2h (list (atan s c)))))) diff --git a/src/high-level/types/complex-single-float.lisp b/src/high-level/types/complex-single-float.lisp new file mode 100644 index 00000000..10daa44b --- /dev/null +++ b/src/high-level/types/complex-single-float.lisp @@ -0,0 +1,75 @@ +;;;; complex-single-float.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftensor tensor/complex-single-float (complex single-float)) + +(defmatrix matrix/complex-single-float (complex single-float) tensor/complex-single-float) + +(defvector vector/complex-single-float (complex single-float) tensor/complex-single-float) + +(defcompatible + (lambda (tensor) + (case (order tensor) + (1 '(vector/complex-single-float + tensor/complex-single-float)) + (2 '(matrix/complex-single-float + tensor/complex-single-float)) + (t '(tensor/complex-single-float)))) + tensor/complex-single-float + matrix/complex-single-float + vector/complex-single-float) + +(defmethod dot ((vector1 vector/complex-single-float) (vector2 vector/complex-single-float)) + (assert (cl:= (size vector1) (size vector2)) + () "Vectors must have the same size. The first vector is size ~a and the second vector is size ~a." + (size vector1) (size vector2)) + (loop :for i :below (size vector1) + :sum (* (tref vector1 i) (conjugate (tref vector2 i))))) + +(defmethod = ((tensor1 tensor/complex-single-float) (tensor2 tensor/complex-single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 matrix/complex-single-float) (tensor2 matrix/complex-single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 vector/complex-single-float) (tensor2 vector/complex-single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (and (<= (abs (- (realpart (apply #'tref tensor1 pos)) + (realpart (apply #'tref tensor2 pos)))) + epsilon) + (<= (abs (- (imagpart (apply #'tref tensor1 pos)) + (imagpart (apply #'tref tensor2 pos)))) + epsilon)) + (return-from = nil)))) + t) diff --git a/src/high-level/types/double-float.lisp b/src/high-level/types/double-float.lisp new file mode 100644 index 00000000..37bec3a1 --- /dev/null +++ b/src/high-level/types/double-float.lisp @@ -0,0 +1,59 @@ +;;;; double-float.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftensor tensor/double-float double-float) + +(defmatrix matrix/double-float double-float tensor/double-float) + +(defvector vector/double-float double-float tensor/double-float) + +(defcompatible + (lambda (tensor) + (case (order tensor) + (1 '(vector/double-float + tensor/double-float)) + (2 '(matrix/double-float + tensor/double-float)) + (t '(tensor/double-float)))) + tensor/double-float + matrix/double-float + vector/double-float) + +(defmethod = ((tensor1 tensor/double-float) (tensor2 tensor/double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 matrix/double-float) (tensor2 matrix/double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 vector/double-float) (tensor2 vector/double-float) &optional (epsilon *double-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) diff --git a/src/high-level/types/int32.lisp b/src/high-level/types/int32.lisp new file mode 100644 index 00000000..54675f0e --- /dev/null +++ b/src/high-level/types/int32.lisp @@ -0,0 +1,23 @@ +;;;; double-float.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftensor tensor/int32 (signed-byte 32)) + +(defmatrix matrix/int32 (signed-byte 32) tensor/int32) + +(defvector vector/int32 (signed-byte 32) tensor/int32) + +(defcompatible + (lambda (tensor) + (case (order tensor) + (1 '(vector/int32 + tensor/int32)) + (2 '(matrix/int32 + tensor/int32)) + (t '(tensor/int32)))) + tensor/int32 + matrix/int32 + vector/int32) diff --git a/src/high-level/types/single-float.lisp b/src/high-level/types/single-float.lisp new file mode 100644 index 00000000..f98698bc --- /dev/null +++ b/src/high-level/types/single-float.lisp @@ -0,0 +1,59 @@ +;;;; single-float.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftensor tensor/single-float single-float) + +(defmatrix matrix/single-float single-float tensor/single-float) + +(defvector vector/single-float single-float tensor/single-float) + +(defcompatible + (lambda (tensor) + (case (order tensor) + (1 '(vector/single-float + tensor/single-float)) + (2 '(matrix/single-float + tensor/single-float)) + (t '(tensor/single-float)))) + tensor/single-float + matrix/single-float + vector/single-float) + +(defmethod = ((tensor1 tensor/single-float) (tensor2 tensor/single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 matrix/single-float) (tensor2 matrix/single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) + +(defmethod = ((tensor1 vector/single-float) (tensor2 vector/single-float) &optional (epsilon *float-comparison-threshold*)) + (unless (equal (shape tensor1) (shape tensor2)) + (return-from = nil)) + (map-indexes + (shape tensor1) + (lambda (&rest pos) + (unless (<= (abs (- (apply #'tref tensor1 pos) + (apply #'tref tensor2 pos))) + epsilon) + (return-from = nil)))) + t) diff --git a/src/high-level/util.lisp b/src/high-level/util.lisp new file mode 100644 index 00000000..69db5eef --- /dev/null +++ b/src/high-level/util.lisp @@ -0,0 +1,97 @@ +;;;; util.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(defvar *float-comparison-threshold* single-float-epsilon) +(defvar *double-comparison-threshold* double-float-epsilon) + +(defun row-major-index (pos dims) + (declare (type index pos) + (type shape dims) + (optimize (speed 3) (safety 0))) + (if (and (cdr dims) (not (cddr dims))) ;; Quick test for length 2 + (matrix-row-major-index (first pos) (second pos) (first dims) (second dims)) + (loop :for i :of-type fixnum in pos + :for d :of-type fixnum in dims + :for acc :of-type fixnum := i + :then (+ i (the fixnum (* d acc))) + :finally (return acc)))) + +(defun column-major-index (pos dims) + (declare (type index pos) + (type shape dims)) + (if (and (cdr dims) (not (cddr dims))) ;; Quick test for length 2 + (matrix-column-major-index (first pos) (second pos) (first dims) (second dims)) + (loop :for i :of-type fixnum in (reverse pos) + :for d :of-type fixnum in (reverse dims) + :for acc :of-type fixnum := i + :then (+ i (the fixnum (* d acc))) + :finally (return acc)))) + +(declaim (inline matrix-row-major-index)) +(defun matrix-row-major-index (row col numrows numcols) + (declare (optimize (speed 3) (safety 0)) + (ignore numrows) + (type fixnum row col numcols) + (values fixnum)) + (+ col (the fixnum (* row numcols)))) + +(declaim (inline matrix-column-major-index)) +(defun matrix-column-major-index (row col numrows numcols) + (declare (optimize (speed 3) (safety 0)) + (ignore numcols) + (type fixnum row col numrows) + (values fixnum)) + (+ row (the fixnum (* col numrows)))) + +(defun from-row-major-index (index dims) + (check-type index fixnum) + (check-type dims shape) + (reverse + (loop :for d :in (reverse dims) + :with acc := index + :collect + (multiple-value-bind (a r) + (floor acc d) + (setf acc a) + r)))) + +(defun from-column-major-index (index dims) + (check-type index fixnum) + (check-type dims shape) + (loop :for d :in dims + :with acc := index + :collect + (multiple-value-bind (a r) + (floor acc d) + (setf acc a) + r))) + +(defun map-indexes (dims f) + "Call a function F for all indexes in shape DIMS, going in row-major order" + (declare (type function f) + (type shape dims)) + (let ((ihead (make-list (list-length dims)))) + (declare (dynamic-extent ihead)) + (labels ((rec (dims itail) + (declare (type list dims itail)) + (cond + ((endp dims) (apply f ihead)) + (t (dotimes (d (the alexandria:non-negative-fixnum (car dims))) + (rplaca itail d) + (rec (cdr dims) (cdr itail))))))) + (rec dims ihead)))) + +(defun map-column-indexes (dims f) + "Call a function F for all indexes in shape DIMS, going in column-major order" + (check-type f function) + (check-type dims shape) + (map-indexes + dims + (lambda (index) + (funcall f (reverse index))))) + +(defun normalize-type (type) + (upgraded-array-element-type type)) diff --git a/src/high-level/vector.lisp b/src/high-level/vector.lisp new file mode 100644 index 00000000..5f58ddef --- /dev/null +++ b/src/high-level/vector.lisp @@ -0,0 +1,168 @@ +;;;; vector.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl) + +(deftype vector-storage (&optional type) + `(simple-array ,type (*))) + +(defstruct (vector (:include abstract-tensor) + (:constructor nil) + (:copier nil)) + (size 0 :type alexandria:positive-fixnum)) + +(defmacro defvector (name type tensor-class) + "Define a new vector subclass with the specified NAME and element TYPE, +compatible with TENSOR-CLASS, as well as the abstract-tensor methods +required not specified by the generic VECTOR class (MAKE-TENSOR, +ELEMENT-TYPE, CAST, COPY-TENSOR, DEEP-COPY-TENSOR, TREF, SETF TREF)" + (let ((constructor-sym (intern (format nil "MAKE-~:@(~A~)" name))) + (copy-sym (intern (format nil "COPY-~:@(~A~)" name))) + (storage-sym (intern (format nil "~:@(~A~)-STORAGE" name)))) + `(progn + (defstruct (,name (:include vector) + (:constructor ,constructor-sym + (size storage)) + (:copier ,copy-sym)) + (storage nil :type (vector-storage ,type))) + #+sbcl (declaim (sb-ext:freeze-type ,name)) + + (defmethod storage ((v ,name)) + (,storage-sym v)) + + (defmethod element-type ((v ,name)) + (declare (ignore v)) + ',type) + + (defmethod make-tensor ((class (eql ',name)) shape &key initial-element layout storage) + (declare (ignore layout)) + (policy-cond:with-expectations (> speed safety) + ((type shape shape) + (assertion (cl:= 1 (length shape)))) + (let ((size (reduce #'* shape))) + (funcall #',constructor-sym + size + (or + storage + (apply #'make-array + size + :element-type ',type + (if initial-element + (list :initial-element (coerce initial-element ',type)) + nil))))))) + + (defmethod cast ((tensor ,name) (class (eql ',name))) + (declare (ignore class)) + tensor) + (defmethod cast :before ((tensor ,tensor-class) (class (eql ',name))) + (declare (ignore class)) + (assert (cl:= 1 (order tensor)) + () + "Cannot change non-1 dimensional tensor to vector.")) + (defmethod cast ((tensor ,tensor-class) (class (eql ',name))) + (declare (ignore class)) + (make-tensor ',name (shape tensor) + :storage (storage tensor))) + (defmethod cast ((tensor ,name) (class (eql ',tensor-class))) + (declare (ignore class)) + (make-tensor ',tensor-class (shape tensor) + :storage (storage tensor))) + + ;; TODO: This does not allow for args. Make this allow for args. + (defmethod copy-tensor ((v ,name) &rest args) + (declare (ignore args)) + (let ((new-v (,copy-sym v))) + (setf (,storage-sym new-v) + (make-array (vector-size v) :element-type (element-type v))) + new-v)) + + (defmethod deep-copy-tensor ((v ,name) &rest args) + (declare (ignore args)) + (let ((new-v (,copy-sym v))) + (setf (,storage-sym new-v) + (copy-seq (,storage-sym v))) + new-v)) + + (defmethod tref ((vector ,name) &rest pos) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-index-p pos (list (vector-size vector))))) + (aref (,storage-sym vector) (first pos)))) + + (defmethod (setf tref) (new-value (vector ,name) &rest pos) + (policy-cond:with-expectations (> speed safety) + ((assertion (valid-index-p pos (list (vector-size vector))))) + (setf (aref (,storage-sym vector) (first pos)) new-value)))))) + + +(defun pprint-vector (stream vector) + "Pretty-print a vector VECTOR to the stream STREAM." + (flet ((print-real (x) + (format stream "~6,3f" x)) + (print-complex (z) + (format stream "~6,3f ~:[+~;-~]~6,3fj" + (realpart z) + (minusp (imagpart z)) + (abs (imagpart z)))) + (print-int (x) + (format stream "~3d" x))) + (let* ((size (size vector)) + (type (element-type vector)) + (print-entry + (cond + ((subtypep type 'complex) #'print-complex) + ((subtypep type 'integer) #'print-int) + (t #'print-real)))) + (pprint-logical-block (stream nil) + (print-unreadable-object (vector stream :type t) + (format stream "(~D):" size) + (dotimes (e size) + (pprint-newline :mandatory stream) + (funcall print-entry (tref vector e)))))))) + +(set-pprint-dispatch 'vector 'pprint-vector) + +;;; Required abstract-tensor methods + +(defmethod order ((vector vector)) + (declare (ignore vector)) + 1) + +(defmethod shape ((vector vector)) + (list (vector-size vector))) + +(defmethod (setf shape) (new-value (vector vector)) + (policy-cond:with-expectations (> speed safety) + ((type shape new-value) + (assertion (cl:= 1 (length new-value)))) + (setf (vector-size vector) (first new-value)))) + +(defgeneric dot (vector1 vector2) + (:documentation "Compute the dot product of two vectors") + (:method ((vector1 vector) (vector2 vector)) + (policy-cond:with-expectations (> speed safety) + ((assertion (cl:= (size vector1) (size vector2)))) + (loop :for i :below (size vector1) + :sum (* (tref vector1 i) (tref vector2 i)))))) + +(deftype p-norm-type () + `(or (member :inf :infinity :positive-infinity) + (integer 1))) + +(defun norm (vector &optional (p 2)) + "Compute the p-norm of a vector. + +If P is not specified then the Euclidean norm is computed. +Special values of P are: :INFINITY :INF :POSITIVE-INFINITY" + (declare (type vector vector) + (type p-norm-type p)) + (case p + (1 + (reduce #'+ (storage vector) :key #'abs)) + ((:infinity :inf :positive-infinity) + (abs (reduce (lambda (x y) (max (abs x) (abs y))) (storage vector)))) + (t + (expt (reduce (lambda (x y) (+ x (abs (expt y p)))) + (storage vector) + :initial-value 0) + (/ p))))) diff --git a/src/packages.lisp b/src/packages.lisp index 6849e6e4..f7a2227d 100644 --- a/src/packages.lisp +++ b/src/packages.lisp @@ -38,62 +38,130 @@ (defpackage #:magicl (:use #:common-lisp - #:cffi) + #:cffi + #:abstract-classes) #+package-local-nicknames (:local-nicknames (#:blas #:magicl.blas-cffi) (#:lapack #:magicl.lapack-cffi) (#:expokit #:magicl.expokit-cffi)) (:import-from #:magicl.foreign-libraries #:print-availability-report) - (:export #:*type-strictness* ; VARIABLE - #:S #:D #:C #:Z ; SYMBOLS - #:matrix ; TYPE, FUNCTION - #:make-matrix ; FUNCTION - #:matrix-rows ; READER - #:matrix-cols ; READER - #:matrix-element-type ; FUNCTION - #:make-zero-matrix ; FUNCTION - #:square-matrix-p ; FUNCTION - #:identityp ; FUNCTION - #:unitaryp ; FUNCTION - #:map-indexes ; FUNCTION - #:tabulate ; FUNCTION - #:make-identity-matrix ; FUNCTION - #:print-availability-report - #:with-blapack - #:make-complex-matrix - #:make-complex-vector - #:conjugate-entrywise - #:transpose - #:conjugate-transpose + (:shadow #:vector + #:= + #:map + #:trace + #:every + #:some + #:notevery + #:notany) + (:export #:with-blapack + + ;; abstract-tensor protocol + #:specialize-tensor + #:generalize-tensor + #:shape + #:tref + #:order + #:size + #:element-type + #:lisp-array + + #:every + #:some + #:notevery + #:notany + + #:map + #:map! + + #:reshape #:slice - #:matrix-column - #:matrix-row - #:qr - #:ql - #:rq - #:lq - #:svd - #:multiply-complex-matrices + + ;; Classes + #:tensor + #:matrix + + ;; Accessors + #:nrows + #:ncols + + ;; Subtypes + #:tensor/single-float + #:tensor/double-float + #:tensor/complex-single-float + #:tensor/complex-double-float + #:matrix/single-float + #:matrix/double-float + #:matrix/complex-single-float + #:matrix/complex-double-float + #:vector/single-float + #:vector/double-float + #:vector/complex-single-float + #:vector/complex-double-float + + ;; Constructors + #:make-tensor + #:empty + #:const + #:rand + #:eye + #:arange + #:from-array + #:from-list + #:from-diag + #:zeros + #:ones + + #:random-unitary + + ;; Operators + #:binary-operator + #:.+ + #:.- + #:.* + #:./ + #:.^ + #:= + #:map + + ;; Matrix operators + #:square-matrix-p + #:identity-matrix-p + #:unitary-matrix-p + #:row + #:column + #:@ + #:mult + #:kron #:scale - #:ref - #:csd - #:lapack-csd + #:scale! + #:diag #:det - #:inv - #:dagger + #:upper-triangular + #:lower-triangular + #:triu + #:tril + #:transpose + #:transpose! + #:orthonormalize + #:orthonormalize! + #:trace #:direct-sum - #:diag - #:matrix-diagonal + #:conjugate-transpose + #:conjugate-transpose! + #:dagger + #:dagger! #:eig #:hermitian-eig - #:logm - #:kron - #:exptm - #:inc-matrix - #:dec-matrix - #:add-matrix - #:sub-matrix + #:inv + #:lu + #:csd + #:svd + #:ql + #:qr + #:rq + #:lq + #:polynomial #:make-polynomial #:polynomial-coefficients @@ -101,10 +169,22 @@ #:polynomial-eval #:polynomial-diff #:polynomial-newton-iteration - ) - ;; random.lisp - (:export #:random-matrix ; FUNCTION - #:random-gaussian-matrix ; FUNCTION - #:random-unitary ; FUNCTION + ;; Vector operators + #:dot + #:norm + + ;; LAPACK stuff + #:lapack-eig + #:lapack-lu + #:lapack-csd + #:lapack-svd + #:lapack-ql + #:lapack-qr + #:lapack-rq + #:lapack-lq + #:lapack-ql-q + #:lapack-qr-q + #:lapack-rq-q + #:lapack-lq-q )) diff --git a/src/with-array-pointers.lisp b/src/with-array-pointers.lisp index 385eff89..204e8e57 100644 --- a/src/with-array-pointers.lisp +++ b/src/with-array-pointers.lisp @@ -25,7 +25,7 @@ "Like LET, but binds with dynamic extent the pointers to the data of numeric simple arrays being bound. WARNING: Do not close over these pointers or otherwise store them outside of the extent of this macro." - (assert (every (lambda (binding) + (assert (cl:every (lambda (binding) (and (listp binding) (= 2 (length binding)) (symbolp (first binding)))) diff --git a/tests/abstract-tensor-tests.lisp b/tests/abstract-tensor-tests.lisp new file mode 100644 index 00000000..9c881db6 --- /dev/null +++ b/tests/abstract-tensor-tests.lisp @@ -0,0 +1,34 @@ +;;;; abstract-tensor-tests.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl-tests) + +(deftest test-tensor-equality () + "Test that tensor equality is sane for tensors of dimension 1 to 8" + (loop :for dimensions :on '(8 7 6 5 4 3 2 1) :do + (loop :for type :in +magicl-types+ :do + (is (magicl:= (magicl:const 1 dimensions :type type) + (magicl:const 1 dimensions :type type))) + (is (not (magicl:= (magicl:const 1 dimensions :type type) + (magicl:const 2 dimensions :type type))))))) + +(deftest test-tensor-shape () + "Test that the shape is returned correctly for tensors of dimension 1 to 8" + (loop :for dimensions :on '(8 7 6 5 4 3 2 1) :do + (is (equal dimensions (magicl:shape (magicl:empty dimensions)))))) + +(deftest test-tensor-order () + "Test that the order is returned correctly for tensors of dimension 1 to 8" + (loop :for dimensions :on '(8 7 6 5 4 3 2 1) + :for i :from 8 :downto 1 :do + (is (equal i (magicl:order (magicl:empty dimensions)))))) + +(deftest test-tensor-tref () + (let ((tensor (magicl:from-list '(1 2 3 4 5 + 6 7 8 9 10 + 11 12 13 14 15) + '(3 5) + :type 'double-float))) + (loop :for i :below 15 + :do (= i (apply #'magicl:tref tensor (magicl::from-row-major-index i '(3 5))))))) diff --git a/tests/constants.lisp b/tests/constants.lisp new file mode 100644 index 00000000..d7ac3041 --- /dev/null +++ b/tests/constants.lisp @@ -0,0 +1,30 @@ +;;;; constants.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl-tests) + +(alexandria:define-constant +magicl-types+ + '(single-float + double-float + (complex single-float) + (complex double-float) + (signed-byte 32)) + :test #'equal) + +(alexandria:define-constant +magicl-float-types+ + '(single-float + double-float + (complex single-float) + (complex double-float)) + :test #'equal) + +(alexandria:define-constant +magicl-matrix-classes+ + '(magicl::matrix/single-float + magicl::matrix/double-float + magicl::matrix/complex-single-float + magicl::matrix/complex-double-float + magicl::matrix/int32) + :test #'equal) + +(defconstant +double-comparison-threshold-loose+ (* 256 double-float-epsilon)) diff --git a/tests/high-level-tests.lisp b/tests/high-level-tests.lisp index c269bc59..6679813e 100644 --- a/tests/high-level-tests.lisp +++ b/tests/high-level-tests.lisp @@ -2,86 +2,109 @@ ;;;; ;;;; Author: Joseph Lin ;;;; Robert Smith +;;;; Cole Scott (in-package #:magicl-tests) (deftest test-determinant () "Test that DET works." - (let* ((x (magicl::make-complex-matrix 3 3 (list 6 4 2 1 -2 8 1 5 7))) - (d (magicl::det x))) - (is (= d -306.0d0)))) + (let* ((x (magicl:from-list '(6 4 2 1 -2 8 1 5 7) '(3 3) + :type 'double-float)) + (d (magicl:det x))) + (is (= d -306d0)))) + +(deftest test-p-norm () + "Test that the p-norm of vectors returns sane values." + ;; Basic 3-4-5 + (is (= 5 (magicl:norm (magicl:from-list '(3 4) '(2))))) + ;; One element vector should return element for all + (loop :for val :in '(-3 0 10) :do + (let ((x (magicl:from-list (list val) '(1)))) + (is (= (abs val) (magicl:norm x 1))) + (is (= (abs val) (magicl:norm x 2))) + (is (= (abs val) (magicl:norm x 3))) + (is (= (abs val) (magicl:norm x :infinity))))) + ;; Test known values + (let ((x (magicl:from-list '(1 -2 3 4 5 -6) '(6)))) + (is (= 6 (magicl:norm x :infinity))) + (is (= 21 (magicl:norm x 1))) + (is (= 9.539392 (magicl:norm x 2))))) (deftest test-examples () "Run all of the examples. Does not check for their correctness." (is (magicl-examples:all-examples))) -(deftest test-identity-checking () - "Check that we can make and identify identities." - (loop :for i :from 1 :to 128 :do - (is (magicl:identityp (magicl:make-identity-matrix i))))) - (deftest test-random-unitary () "Check that we can make and identify unitaries." (loop :for i :from 1 :to 128 :do - (is (magicl:unitaryp (magicl:random-unitary i))))) + (is (magicl:unitary-matrix-p (magicl:random-unitary (list i i) :type '(complex double-float)) +double-comparison-threshold-loose+)))) (deftest test-logm () "Check that the matrix logarithm is the inverse of the matrix exponential." - (let ((x (magicl:random-unitary 4))) - (loop :for i :from 0 :to (1- (magicl:matrix-cols x)) - :for j :from 0 :to (1- (magicl:matrix-rows x)) - :do (let ((diff (- (magicl:ref x i j) (magicl:ref (magicl-transcendental:expm - (magicl-transcendental:logm x)) - i j))) - (eps .00001f0)) + (let ((x (magicl:random-unitary '(4 4) :type '(complex double-float)))) + (loop :for i :from 0 :to (1- (magicl:ncols x)) + :for j :from 0 :to (1- (magicl:nrows x)) + :do (let ((diff (- (magicl:tref x i j) + (magicl:tref (magicl-transcendental:expm + (magicl-transcendental:logm x)) + i j))) + (eps .00001f0)) (is (< (abs diff) eps)))))) (deftest test-kron () "Test a few properties of the kronecker product." (let* ((matrix-dim 2) - (eye2 (make-identity-matrix matrix-dim)) + (eye2 (magicl:eye '(2 2))) (eye4 (magicl:kron eye2 eye2)) - (x (magicl:make-complex-matrix matrix-dim matrix-dim '(0 1 1 0))) + (x (magicl:from-list '(0.0 1.0 1.0 0.0) '(2 2))) (xxx (magicl:kron x x x)) - (xxx-cols (matrix-cols xxx))) + (xxx-cols (magicl:ncols xxx))) ;; Check that the kronecker product of two small identities is a larger identity. - (is (magicl:identityp eye4)) + (is (magicl:identity-matrix-p eye4)) + ;; Check that XXX is evaluated correctly. - (is (loop :for i :below (matrix-rows xxx) :always - (loop :for j :below xxx-cols - :always (if (= (- xxx-cols 1) (+ i j)) - (= 1 (ref xxx i j)) - (= 0 (ref xxx i j)))))) + (is (loop :for i :below (magicl:nrows xxx) + :always + (loop :for j :below xxx-cols + :always (if (cl:= (cl:- xxx-cols 1) (cl:+ i j)) + (cl:= 1 (magicl:tref xxx i j)) + (cl:= 0 (magicl:tref xxx i j)))))) ;; Check that IX /= XI. - (is (not (equalp (magicl:kron eye2 x) (kron x eye2)))) + (is (not (magicl:= (magicl:kron eye2 x) (magicl:kron x eye2)))) ;; Check that it is correctly working on more than two inputs. XXX /= XXI - (is (not (equalp (magicl:kron x x x ) (kron x x eye2)))) + (is (not (magicl:= (magicl:kron x x x) (magicl:kron x x eye2)))) ;; Check that one of the two gives the correct analytic result. - (is (equalp (kron x eye2) (make-complex-matrix 4 4 '(0 0 1 0 - 0 0 0 1 - 1 0 0 0 - 0 1 0 0)))) - ;; Check that it yields teh right dimensions. - (is (= (matrix-rows xxx) (matrix-cols xxx) (expt matrix-dim 3))) + (is (magicl:= (magicl:kron x eye2) + (magicl:from-list '(0 0 1 0 + 0 0 0 1 + 1 0 0 0 + 0 1 0 0) + '(4 4) + :type '(complex double-float)))) + ;; Check that it yields the right dimensions. + (is (cl:= (magicl:nrows xxx) (magicl:ncols xxx) (expt matrix-dim 3))) )) + (deftest test-svd () "Test the full and reduced SVDs." (labels ((mul-diag-times-gen (diag matrix) "Returns a newly allocated matrix resulting from the product of DIAG (a diagonal real matrix) with MATRIX (a complex matrix)." + #+ignore (declare (type matrix diag matrix) (values matrix)) - (let* ((m (matrix-rows diag)) - (k (matrix-cols matrix)) - (result (make-zero-matrix m k))) - (dotimes (i (min m (matrix-cols diag)) result) - (let ((dii (ref diag i i))) + (let* ((m (magicl:nrows diag)) + (k (magicl:ncols matrix)) + (result (magicl:empty (list m k)))) + (dotimes (i (min m (magicl:ncols diag)) result) + (let ((dii (magicl:tref diag i i))) (dotimes (j k) - (setf (ref result i j) (* dii (ref matrix i j)))))))) + (setf (magicl:tref result i j) + (* dii (magicl:tref matrix i j)))))))) (norm-inf (matrix) "Return the infinity norm of vec(MATRIX)." - (let ((data (magicl::matrix-data matrix))) + (let ((data (magicl::storage matrix))) (reduce #'max data :key #'abs))) (zero-p (matrix &optional (tolerance (* 1.0d2 double-float-epsilon))) @@ -90,57 +113,60 @@ (check-full-svd (matrix) "Validate full SVD of MATRIX." - (let ((m (matrix-rows matrix)) - (n (matrix-cols matrix))) + (let ((m (magicl:nrows matrix)) + (n (magicl:ncols matrix))) (multiple-value-bind (u sigma vh) - (svd matrix) - (is (= (matrix-rows u) (matrix-cols u) m)) - (is (and (= (matrix-rows sigma) m) (= (matrix-cols sigma) n))) - (is (= (matrix-rows vh) (matrix-cols vh) n)) - (is (zero-p (sub-matrix matrix (multiply-complex-matrices u (mul-diag-times-gen sigma vh)))))))) + (magicl:svd matrix) + (is (= (magicl:nrows u) (magicl:ncols u) m)) + (is (and (= (magicl:nrows sigma) m) (= (magicl:ncols sigma) n))) + (is (= (magicl:nrows vh) (magicl:ncols vh) n)) + (is (zero-p (magicl:.- matrix (magicl:@ u (mul-diag-times-gen sigma vh)))))))) (check-reduced-svd (matrix) "Validate reduced SVD of MATRIX." - (let* ((m (matrix-rows matrix)) - (n (matrix-cols matrix)) + (let* ((m (magicl:nrows matrix)) + (n (magicl:ncols matrix)) (k (min m n))) (multiple-value-bind (u sigma vh) - (svd matrix :reduced t) - (is (and (= (matrix-rows u) m) - (= (matrix-cols u) k))) - (is (= (matrix-rows sigma) (matrix-cols sigma) k)) - (is (and (= (matrix-rows vh) k) - (= (matrix-cols vh) n))) - (is (zero-p (sub-matrix matrix (multiply-complex-matrices u (mul-diag-times-gen sigma vh))))))))) - - (let ((tall-thin-matrix (random-matrix 8 2))) + (magicl:svd matrix :reduced t) + (is (and (= (magicl:nrows u) m) + (= (magicl:ncols u) k))) + (is (= (magicl:nrows sigma) (magicl:ncols sigma) k)) + (is (and (= (magicl:nrows vh) k) + (= (magicl:ncols vh) n))) + (is (zero-p (magicl:.- matrix (magicl:@ u (mul-diag-times-gen sigma vh))))))))) + + (let ((tall-thin-matrix (magicl:rand '(8 2)))) (check-full-svd tall-thin-matrix) (check-reduced-svd tall-thin-matrix)) - (let ((short-fat-matrix (random-matrix 2 8))) + (let ((short-fat-matrix (magicl:rand '(2 8)))) (check-full-svd short-fat-matrix) (check-reduced-svd short-fat-matrix)))) + (deftest test-csd-2x2-basic () "Test CS decomposition of an equipartitioned 2x2 unitary matrix." - (let ((x (magicl:random-unitary 2)) + (let ((x (magicl:random-unitary '(2 2) + :type '(complex double-float))) (tol (* 1.0d2 double-float-epsilon))) (multiple-value-bind (u1 u2 v1h v2h theta) (magicl::csd-2x2-basic x 1 1) (multiple-value-bind (u1* u2* v1h* v2h* theta*) - (lapack-csd x 1 1) - (is (< (abs (- (ref u1 0 0) (ref u1* 0 0))) tol)) - (is (< (abs (- (ref u2 0 0) (ref u2* 0 0))) tol)) - (is (< (abs (- (ref v1h 0 0) (ref v1h* 0 0))) tol)) - (is (< (abs (- (ref v2h 0 0) (ref v2h* 0 0))) tol)) + (magicl:lapack-csd x 1 1) + (is (< (abs (- (magicl:tref u1 0 0) (magicl:tref u1* 0 0))) tol)) + (is (< (abs (- (magicl:tref u2 0 0) (magicl:tref u2* 0 0))) tol)) + (is (< (abs (- (magicl:tref v1h 0 0) (magicl:tref v1h* 0 0))) tol)) + (is (< (abs (- (magicl:tref v2h 0 0) (magicl:tref v2h* 0 0))) tol)) (is (< (abs (- (first theta) (first theta*))) tol)))))) (deftest test-polynomial-solver () "Test univariate polynomial solver." ;; Test random polynomials with favorable coefficients. (flet ((make-random-polynomial (degree) - (let ((c (magicl::matrix-data (magicl:random-matrix 1 degree)))) + (let ((c (magicl::storage (magicl:rand (list 1 degree) + :type '(complex double-float))))) (setf (aref c (1- degree)) (complex 1.0d0)) (magicl::%make-polynomial :coefficients c)))) (dotimes (i 10) diff --git a/tests/matrix-tests.lisp b/tests/matrix-tests.lisp new file mode 100644 index 00000000..19f49f07 --- /dev/null +++ b/tests/matrix-tests.lisp @@ -0,0 +1,72 @@ +;;;; tests/matrix-tests.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl-tests) + +(defmacro is-matrix (&rest rest) + `(progn ,@(loop :for m :in rest + :collect `(is (subtypep (type-of ,m) 'matrix))))) + +(defmacro is-not-matrix (&rest rest) + `(progn ,@(loop :for m :in rest + :collect `(is (not (subtypep (type-of ,m) 'matrix)))))) + +(deftest test-identity-matrix-p () + "Test that identity matrices can be identified by IDENTITY-MATRIX-P for all types of matrixes from 1x1 to 64x64" + (loop :for i :from 1 :to 64 + :do (loop :for type :in +magicl-types+ + :do (is (magicl:identity-matrix-p (magicl:eye (list i i) :type type))) + (is (not (magicl:identity-matrix-p (magicl:eye (list i i) :value 2 :type type)))) + (is (not (magicl:identity-matrix-p (magicl:const 0 (list i i) :type type))))))) + +(deftest test-square-matrix-p () + "Test that square matrices can be identified by IDENTITY-MATRIX-P for all types of matrixes from 1x1 to 64x64" + (loop :for i :from 1 :to 64 + :do (loop :for type :in +magicl-types+ + :do (is (magicl:square-matrix-p (magicl:empty (list i i) :type type))) + (is (not (magicl:square-matrix-p (magicl:empty (list i (* 2 i)) :type type))))))) + +;; Multiplication + +(deftest test-matrix-multiplication-errors () + (signals simple-error (magicl:@ + (magicl:empty '(3 3)) + (magicl:empty '(1 1)))) + (signals simple-error (magicl:@ + (magicl:empty '(1 2)) + (magicl:empty '(1 2)))) + (signals simple-error (magicl:@ + (magicl:empty '(5 2)) + (magicl:empty '(2 3)) + (magicl:empty '(2 3))))) + +(deftest test-complex-matrix-multiplication-results () + "Test a few basic complex matrix multiplications" + (let* ((m-old (magicl:from-list '(#C(1d0 2d0) #C(3d0 4d0) #C(5d0 6d0) #C(7d0 8d0)) '(2 2) :layout :row-major)) + (m (magicl:from-list '(#C(1d0 2d0) #C(3d0 4d0) #C(5d0 6d0) #C(7d0 8d0)) '(2 2) :layout :row-major)) + (x-old (magicl:from-list '(#C(1d0 2d0) #C(3d0 4d0)) '(2 1))) + (x (magicl:from-list '(#C(1d0 2d0) #C(3d0 4d0)) '(2 1))) + (expected (magicl:from-list '(#C(-10d0 28d0) #C(-18d0 68d0)) '(2 1)))) + ;; Check that the multiplication is correct and does not throw any errors + (is (magicl:= expected (magicl:@ m x))) + + ;; Check that the multiplication did not modify the inputs + (is (magicl:= m-old m)) + (is (magicl:= x-old x)) + + ;; Check that doing 2x1 @ 2x2 errors + (signals error (magicl:@ x m)))) + +(deftest test-random-unitary-properties () + "Test calls to RANDOM-UNITARY for all float types and sizes 1x1 to 64x64 to check properties" + (loop :for type :in +magicl-float-types+ + :do (loop :for i :from 1 :to 64 + :do (let ((m (magicl:random-unitary (list i i) :type type))) + (is (> 5e-5 (abs (cl:- + (abs (magicl:det m)) + 1)))) + (is (magicl:= + (magicl:eye (list i i) :type type) + (magicl:@ m (magicl:transpose m)) + 5e-5)))))) diff --git a/tests/packages.lisp b/tests/packages.lisp index e51cc1aa..c85341ec 100644 --- a/tests/packages.lisp +++ b/tests/packages.lisp @@ -3,8 +3,6 @@ ;;;; Author: Joseph Lin (fiasco:define-test-package #:magicl-tests - (:use #:magicl) - ;; suite.lisp (:export #:run-magicl-tests)) diff --git a/tests/specialization-tests.lisp b/tests/specialization-tests.lisp new file mode 100644 index 00000000..9117c93f --- /dev/null +++ b/tests/specialization-tests.lisp @@ -0,0 +1,31 @@ +;;;; specialization-tests.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl-tests) + +(defmacro is-subtype (type &rest rest) + `(progn ,@(loop :for m :in rest + :collect `(is (subtypep (type-of ,m) ',type))))) + +(defmacro is-not-subtype (type &rest rest) + `(progn ,@(loop :for m :in rest + :collect `(is (not (subtypep (type-of ,m) ',type)))))) + +(deftest test-tensor-shape-specialization () + (let ((vector (magicl:empty '(1))) + (matrix (magicl:empty '(1 2))) + (tensor (magicl:empty '(1 2 3)))) + (is-subtype magicl::vector vector) + (is-subtype magicl::matrix matrix) + (is-subtype magicl::tensor tensor) + (is-not-subtype magicl::vector matrix tensor) + (is-not-subtype magicl::matrix vector tensor) + (is-not-subtype magicl::tensor vector matrix))) + +(deftest test-tensor-type-specialization () + (loop :for type :in +magicl-types+ + :for class :in +magicl-matrix-classes+ + :do (is (equalp + class + (class-name (class-of (magicl:empty '(2 2) :type type))))))) diff --git a/tests/util-tests.lisp b/tests/util-tests.lisp new file mode 100644 index 00000000..812fb6cc --- /dev/null +++ b/tests/util-tests.lisp @@ -0,0 +1,60 @@ +;;;; util-tests.lisp +;;;; +;;;; Author: Cole Scott + +(in-package #:magicl-tests) + +(deftest test-row-major-index () + (let ((rows 3) + (cols 5)) + (loop :for i :below (* rows cols) + :do (is (cl:= i (magicl::row-major-index + (magicl::from-row-major-index i (list rows cols)) + (list rows cols))))) + (let ((examples '((0 (0 0)) + (1 (0 1)) + (4 (0 4)) + (5 (1 0)) + (7 (1 2))))) + (loop :for ex :in examples + :do (is (cl:= (first ex) + (magicl::row-major-index + (second ex) + (list rows cols)))) + (is (equal (second ex) + (magicl::from-row-major-index + (first ex) + (list rows cols)))))))) + +(deftest test-column-major-index () + (let ((rows 3) + (cols 5)) + (loop :for i :below (* rows cols) + :do (is (cl:= i (magicl::column-major-index + (magicl::from-column-major-index i (list rows cols)) + (list rows cols))))) + (let ((examples '((0 (0 0)) + (1 (1 0)) + (4 (1 1)) + (5 (2 1)) + (7 (1 2))))) + (loop :for ex :in examples + :do (is (cl:= (first ex) + (magicl::column-major-index + (second ex) + (list rows cols)))) + (is (equal (second ex) + (magicl::from-column-major-index + (first ex) + (list rows cols)))))))) + +(deftest test-map-indexes () + (let* ((rows 3) + (cols 5) + (arr (make-array (list rows cols) :element-type 'double-float :initial-element 0d0))) + (magicl::map-indexes (list rows cols) + (lambda (&rest pos) + (incf (apply #'aref arr pos)))) + (loop :for r :below rows + :do (loop :for c :below cols + :do (is (cl:= 1 (aref arr r c))))))) diff --git a/transcendental/package.lisp b/transcendental/package.lisp index 740377cb..89052da4 100644 --- a/transcendental/package.lisp +++ b/transcendental/package.lisp @@ -11,7 +11,6 @@ (defpackage #:magicl-transcendental (:use #:common-lisp - #:cffi - #:magicl) + #:cffi) (:export #:expm #:logm)) diff --git a/transcendental/transcendental.lisp b/transcendental/transcendental.lisp index 284ec4d7..ea47ab7d 100644 --- a/transcendental/transcendental.lisp +++ b/transcendental/transcendental.lisp @@ -7,15 +7,15 @@ (defun expm (m) "Finds the exponential of a square matrix M." (let ((ideg 6) - (rows (matrix-rows m)) + (rows (magicl:nrows m)) (tcoef (coerce 1.0 'double-float)) - (h (magicl::copy-matrix-storage (magicl::matrix-data m))) + (h (copy-seq (magicl::storage m))) (iexph 0) (ns 0) (iflag 0)) (let ((lwsp (+ (* 4 rows rows) ideg 1)) - (ipiv (magicl::make-int32-storage rows))) - (let ((wsp (magicl::make-Z-storage lwsp))) + (ipiv (magicl::make-array (list rows) :element-type '(signed-byte 32)))) + (let ((wsp (magicl::make-array (list lwsp) :element-type '(complex double-float)))) ;; Requires direct foreign function call due to need to access a pointer ;; to an integer (IEXPH). (CFFI:WITH-FOREIGN-OBJECTS ((IDEG-REF103 ':INT32) (M-REF104 ':INT32) @@ -39,18 +39,18 @@ IPIV-ptr IEXPH-REF111 NS-REF112 IFLAG-REF113)) (setf iexph (CFFI:MEM-REF IEXPH-REF111 :INT32))) - (let ((exph (magicl::make-Z-storage (* rows rows)))) + (let ((exph (magicl::make-array (list (* rows rows)) :element-type '(complex double-float)))) (dotimes (i (* rows rows)) (setf (row-major-aref exph i) (row-major-aref wsp (+ i (1- iexph))))) - (values (make-matrix :rows rows :cols rows :data exph))))))) + (values (magicl:from-array exph (list rows rows) :layout :column-major))))))) (defun logm (m) "Finds the matrix logarithm of a given square matrix M assumed to be diagonalizable, with nonzero eigenvalues" - (multiple-value-bind (vals vects) (eig m) + (multiple-value-bind (vals vects) (magicl:eig m) (let ((new-log-diag (let ((log-vals (mapcar #'log vals))) - (diag (matrix-cols m) (matrix-rows m) log-vals)))) - (multiply-complex-matrices - vects - (multiply-complex-matrices new-log-diag (inv vects)))))) + (magicl:from-diag log-vals)))) + (magicl:@ vects + new-log-diag + (magicl:inv vects)))))