mps.lisp (6166B)
1 (defpackage :mps 2 (:use :common-lisp :blas :mpi :tensor) 3 (:export :make-submps 4 :submps-start 5 :submps-size 6 :submps-tensors 7 :submps-singvals 8 :left-open-boundary-tensor 9 :right-open-boundary-tensor 10 :mps-prepare 11 :mps-overlap 12 :mps-expectation-value 13 :save-metadata 14 :save-submps 15 :save-mps 16 :state-parameter 17 :state-tensor 18 :load-mps)) 19 20 (in-package :mps) 21 22 (defstruct submps 23 start 24 size 25 tensors 26 singvals) 27 28 (defun left-open-boundary-tensor (operator) 29 (let* ((i (first (tensor-indices operator))) 30 (d (segment-dimension (find-if #'list-zerop i :key #'segment-numbers)))) 31 (functional-tensor 32 (list (list (make-segment :numbers '(0) :dimension 1)) 33 (conjugate-index i)) 34 #'(lambda (s) 35 (let ((u (subscript-subscript (second s)))) 36 (if (= u (1- d)) 37 1 38 0)))))) 39 40 (defun right-open-boundary-tensor (operator) 41 (let* ((i (fourth (tensor-indices operator)))) 42 (functional-tensor 43 (list (conjugate-index i) 44 (list (make-segment :numbers '(0) :dimension 1))) 45 #'(lambda (s) 46 (if (zerop (subscript-subscript (first s))) 47 1 48 0))))) 49 50 (defun right-singvals-tensor (index) 51 (functional-tensor 52 (list 53 (mapcar #'(lambda (s) 54 (make-segment :numbers (list-- (segment-numbers s)) 55 :dimension (segment-dimension s))) 56 index) 57 index) 58 #'(lambda (s) (declare (ignore s)) 1))) 59 60 (defun mps-prepare (rank length submps maxdim) 61 (let* ((s (submps-start submps)) 62 (d (submps-size submps)) 63 (u (submps-tensors submps)) 64 (v (make-array (length u))) 65 l w r) 66 (when (< s length) 67 (setf r (if (zerop s) 68 (let ((m (list (make-segment :numbers '(0) :dimension 1)))) 69 (functional-tensor 70 (list m m) 71 #'(lambda (s) (declare (ignore s)) 1))) 72 (mpi-receive-object (1- rank)))) 73 (loop for i from 0 below d do 74 (setf (svref u i) (tensor-contract r 1 (svref u i) 0)) 75 (multiple-value-setq (l w r) (tensor-decompose (svref u i) '(0 1) maxdim :left)) 76 (setf (svref u i) l) 77 (setf (svref v i) w)) 78 (if (= (+ s d) length) 79 (setf (svref v (1- d)) (right-singvals-tensor (nth 2 (tensor-indices (svref u (1- d)))))) 80 (mpi-send-object r (1+ rank)))) 81 (setf (submps-singvals submps) v))) 82 83 (defun mps-overlap (size rank length submps1 submps2) 84 (let ((s (submps-start submps1)) 85 (d (submps-size submps1)) 86 (u1 (submps-tensors submps1)) 87 (u2 (submps-tensors submps2)) 88 l r) 89 (setf l (if (zerop s) 90 (let ((m (list (make-segment :numbers '(0) :dimension 1)))) 91 (functional-tensor (list m m) 92 #'(lambda (s) (declare (ignore s)) 1))) 93 (mpi-receive-object (1- rank)))) 94 (loop for i from 0 below d do 95 (setf l (tensor-contract 96 (tensor-conjugate (svref u2 i)) 97 '(0 1) 98 (tensor-contract l 1 (svref u1 i) 0) 99 '(0 1)))) 100 (if (and (not (zerop d)) 101 (= (+ s d) length)) 102 (progn 103 (setf r (functional-tensor 104 (list (nth 2 (tensor-indices (svref u2 (1- d)))) 105 (conjugate-index (nth 2 (tensor-indices (svref u1 (1- d)))))) 106 #'(lambda (s) (declare (ignore s)) 1))) 107 (setf l (tensor-contract l '(0 1) r '(0 1))))) 108 (if (= rank (1- size)) 109 (tensor-scalar l) 110 (mpi-send-object l (1+ rank))))) 111 112 (defun mps-expectation-value (size rank length mpo submps) 113 (let ((s (submps-start submps)) 114 (d (submps-size submps)) 115 (u (submps-tensors submps)) 116 l r) 117 (setf l (if (zerop s) 118 (let ((m (list (make-segment :numbers '(0) :dimension 1)))) 119 (functional-tensor (list m m m) 120 #'(lambda (s) (declare (ignore s)) 1))) 121 (mpi-receive-object (1- rank)))) 122 (loop 123 for i from 0 below d 124 for j from s 125 do 126 (setf l (tensor-contract 127 (tensor-conjugate (svref u i)) 128 '(0 1) 129 (tensor-contract (svref mpo j) '(0 2) 130 (tensor-contract l 2 (svref u i) 0) '(1 2)) 131 '(2 0)))) 132 (if (and (not (zerop d)) 133 (= (+ s d) length)) 134 (progn 135 (setf r (functional-tensor 136 (list (nth 2 (tensor-indices (svref u (1- d)))) 137 (list (make-segment :numbers '(0) :dimension 1)) 138 (conjugate-index (nth 2 (tensor-indices (svref u (1- d)))))) 139 #'(lambda (s) (declare (ignore s)) 1))) 140 (setf l (tensor-contract l '(0 1 2) r '(0 1 2))))) 141 (if (= rank (1- size)) 142 (tensor-scalar l) 143 (mpi-send-object l (1+ rank))))) 144 145 (defun save-metadata (dimension length particles pathname) 146 (with-open-file (s (ensure-directories-exist (merge-pathnames #p"metadata" pathname)) 147 :direction :output 148 :if-exists :supersede) 149 (prin1 (list 150 (cons :site-dimension dimension) 151 (cons :model-length length) 152 (cons :particle-number particles)) 153 s) 154 (terpri s))) 155 156 (defun save-tensor (tensor pathname index) 157 (with-open-file (s (merge-pathnames (pathname (format nil "tensor.~A" index)) pathname) 158 :direction :output 159 :if-exists :supersede) 160 (prin1 tensor s) 161 (terpri s))) 162 163 (defun save-submps (submps pathname) 164 (with-standard-io-syntax 165 (loop 166 repeat (submps-size submps) 167 for i from (submps-start submps) 168 for a across (submps-tensors submps) 169 do (save-tensor a pathname i)))) 170 171 (defun save-mps (mps dimension particles pathname) 172 (with-standard-io-syntax 173 (let ((*print-pretty* t)) 174 (when (or (not (probe-file pathname)) (y-or-n-p "Overwrite state ~S?" pathname)) 175 (save-metadata dimension (length mps) particles pathname) 176 (loop 177 for i below (length mps) 178 do (save-tensor (svref mps i) pathname i)))))) 179 180 (defun state-parameter (pathname parameter) 181 (cdr (assoc parameter 182 (with-open-file (s (merge-pathnames #p"metadata" pathname)) 183 (read s))))) 184 185 (defun state-tensor (pathname index) 186 (with-standard-io-syntax 187 (with-open-file (s (merge-pathnames (pathname (format nil "tensor.~A" index)) pathname)) 188 (read s)))) 189 190 (defun load-submps (pathname start end) 191 (loop 192 with d = (- end start) 193 with a = (make-array (1+ d)) 194 for i from 0 195 for j from start below end 196 do (setf (svref a i) (state-tensor pathname j)) 197 finally (return (make-submps :start start :size d :tensors a)))) 198 199 (defun load-mps (rank partition pathname) 200 (let ((s (if (zerop rank) 0 (nth (1- rank) partition))) 201 (e (nth rank partition))) 202 (load-submps pathname s e)))