tebdol

Simulation of ultracold atoms in optical lattices
git clone https://miroslavurbanek.com/git/tebdol.git
Log | Files | Refs | README

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))))