exp.lisp (2531B)
1 (defpackage :exp 2 (:use :common-lisp :sb-ext :blas :array) 3 (:export :hermitian-matrix-exponential)) 4 5 (in-package :exp) 6 7 (defun hermitian-matrix-eigenvalues (matrix) 8 ;; do not destroy the original matrix 9 (setf matrix (array-copy matrix)) 10 (unless (and (= (array-rank matrix) 2) 11 (= (array-dimension matrix 0) (array-dimension matrix 1))) 12 (error "Matrix has incorrect dimensions.")) 13 (let* ((n (array-dimension matrix 0)) 14 (abstol (dlamch "s")) 15 (w (make-double-array n)) 16 (z (make-blas-array (list n n))) 17 (isuppz (make-integer-array (* 2 n))) 18 (work (make-blas-array 1)) 19 (lwork -1) 20 (rwork (make-double-array 1)) 21 (lrwork -1) 22 (iwork (make-integer-array 1)) 23 (liwork -1)) 24 (sb-sys:with-pinned-objects 25 ((array-storage-vector matrix) 26 (array-storage-vector w) 27 (array-storage-vector z) 28 (array-storage-vector isuppz)) 29 (let ((ba (double-array-alien matrix)) 30 (bw (double-array-alien w)) 31 (bz (double-array-alien z)) 32 (bisuppz (integer-array-alien isuppz))) 33 (macrolet ((m () 34 ;; todo: zheevr returns info parameter, check it 35 `(sb-sys:with-pinned-objects 36 ((array-storage-vector work) 37 (array-storage-vector rwork) 38 (array-storage-vector iwork)) 39 (let ((bwork (double-array-alien work)) 40 (brwork (double-array-alien rwork)) 41 (biwork (integer-array-alien iwork))) 42 (zheevr "v" "a" "u" n ba n 0d0 0d0 0 0 abstol 0 bw bz n 43 bisuppz bwork lwork brwork lrwork biwork liwork))))) 44 (m) 45 (setf lwork (floor (realpart (aref work 0)))) 46 (setf work (make-blas-array lwork)) 47 (setf lrwork (floor (aref rwork 0))) 48 (setf rwork (make-double-array lrwork)) 49 (setf liwork (aref iwork 0)) 50 (setf iwork (make-integer-array liwork)) 51 (m)))) 52 ;; m = v*e*v'; return (e v') 53 (values w z))) 54 55 ;; diagonal-matrix-* is a destructive operation 56 57 (defun diagonal-matrix-* (diagonal matrix) 58 (let ((m (array-dimension matrix 0)) 59 (n (array-dimension matrix 1))) 60 (dotimes (i m matrix) 61 (let ((x (aref diagonal i))) 62 (dotimes (j n) 63 (setf (aref matrix i j) (* (aref matrix i j) x))))))) 64 65 (defun hermitian-matrix-exponential (matrix &optional (parameter 1)) 66 (multiple-value-bind (e w) (hermitian-matrix-eigenvalues matrix) 67 (let* ((n (array-dimension e 0)) 68 (f (make-blas-array n)) 69 (v (array-conjugate w))) 70 (dotimes (i n) 71 (setf (aref f i) 72 (coerce (exp (* parameter (aref e i))) 73 '(complex double-float)))) 74 (diagonal-matrix-* f w) 75 (array-contract v 0 w 0))))