tebd.lisp (9175B)
1 (defpackage :tebd 2 (:use :common-lisp :blas :tensor :mpi :mps) 3 (:export :make-single-site-propagator 4 :make-double-site-propagator 5 :make-mph-propagator 6 :make-uniform-propagator 7 :tebd-evolve)) 8 9 (in-package :tebd) 10 11 (defun left-tensor (operator) 12 (let* ((i (first (tensor-indices operator))) 13 (d (segment-dimension (find-if #'list-zerop i :key #'segment-numbers)))) 14 (functional-tensor 15 (list i) 16 #'(lambda (s) 17 (let ((u (first s))) 18 (if (and (list-zerop (subscript-numbers u)) 19 (= (subscript-subscript u) (1- d))) 20 1 21 0)))))) 22 23 (defun right-tensor (operator) 24 (let* ((i (fourth (tensor-indices operator)))) 25 (functional-tensor 26 (list i) 27 #'(lambda (s) 28 (let ((u (first s))) 29 (if (and (list-zerop (subscript-numbers u)) 30 (= (subscript-subscript u)) 0) 31 1 32 0)))))) 33 34 (defun central-tensor (operator) 35 (let* ((i (first (tensor-indices operator))) 36 (d (segment-dimension (find-if #'list-zerop i :key #'segment-numbers)))) 37 (functional-tensor 38 (list i 39 (reverse 40 (mapcar #'(lambda (s) 41 (make-segment :numbers (list-- (segment-numbers s)) 42 :dimension (segment-dimension s))) 43 i))) 44 #'(lambda (s) 45 (let ((u (subscript-subscript (first s))) 46 (v (subscript-subscript (second s)))) 47 (if (= u v) 48 (if (= u (1- d)) 49 0 50 1) 51 0)))))) 52 53 (defun make-single-site-propagator (operator step) 54 (let* ((l (left-tensor operator)) 55 (r (right-tensor operator)) 56 (h (tensor-contract l 0 (tensor-contract operator 3 r 0) 0))) 57 (tensor-exponential h '((0) (1)) (* -1 +blas-i+ step)))) 58 59 (defun make-double-site-propagator (operator1 operator2 step) 60 (let* ((l (left-tensor operator1)) 61 (r (right-tensor operator2)) 62 (c (central-tensor operator2)) 63 (h (tensor-contract 64 l 65 0 66 (tensor-contract 67 operator1 68 3 69 (tensor-contract c 1 (tensor-contract operator2 3 r 0) 0) 70 0) 71 0))) 72 (tensor-exponential h '((0 2) (1 3)) (* -1 +blas-i+ step)))) 73 74 (defun make-single-propagators (mph step) 75 (let ((f (svref mph 0)) 76 (g (svref mph (1- (length mph))))) 77 (list (if f (make-single-site-propagator f step)) 78 (if g (make-single-site-propagator g step))))) 79 80 (defun make-double-propagators (mph step) 81 (loop 82 with l = (1- (length mph)) 83 with a = (make-array l) 84 for i below l 85 for f = (svref mph i) 86 for g = (svref mph (1+ i)) 87 when (and f g) 88 do (setf (svref a i) (make-double-site-propagator f g step)) 89 finally (return a))) 90 91 (defstruct propagator 92 single1 93 double1 94 single2 95 double2) 96 97 (defun make-mph-propagator (mph step) 98 (make-propagator 99 :single1 (make-single-propagators mph step) 100 :double1 (make-double-propagators mph step) 101 :single2 (make-single-propagators mph (/ step 2)) 102 :double2 (make-double-propagators mph (/ step 2)))) 103 104 (defun make-uniform-propagator (length hamiltonian step) 105 (let ((s1 (make-single-site-propagator hamiltonian step)) 106 (d1 (make-double-site-propagator hamiltonian hamiltonian step)) 107 (s2 (make-single-site-propagator hamiltonian (/ step 2))) 108 (d2 (make-double-site-propagator hamiltonian hamiltonian (/ step 2)))) 109 (make-propagator 110 :single1 (list s1 s1) 111 :double1 (make-array (1- length) :initial-element d1) 112 :single2 (list s2 s2) 113 :double2 (make-array (1- length) :initial-element d2)))) 114 115 (defun apply-single-site-operator (operator tensor) 116 (tensor-permute (tensor-contract operator 1 tensor 1) '(1 0 2))) 117 118 (defun apply-double-site-operator (operator tensor1 tensor2 singvals maxdim) 119 (let ((u (tensor-contract 120 operator 121 '(1 3) 122 (tensor-contract tensor1 2 tensor2 0) 123 '(1 2)))) 124 (multiple-value-bind (l s r) 125 (tensor-decompose (tensor-contract u 3 singvals 0) '(2 0) maxdim :left) 126 (declare (ignore r)) 127 (values l s (tensor-contract (tensor-conjugate l) '(0 1) u '(2 0)))))) 128 129 (defun apply-even-tebd-operators (length single double submps maxdim &key (subwt nil) (gc nil)) 130 (let* ((s (submps-start submps)) 131 (d (submps-size submps)) 132 (e (+ s d)) 133 (u (submps-tensors submps)) 134 (v (submps-singvals submps))) 135 (if (and (= e length) (oddp e) (not (zerop d))) 136 (let ((w (mpi-wtime))) 137 (setf (svref u (1- d)) (apply-single-site-operator (cadr single) (svref u (1- d)))) 138 (when subwt 139 (incf (svref subwt (1- d)) (- (mpi-wtime) w))))) 140 (loop 141 for i from (if (evenp s) 0 1) below (if (= e length) (1- d) d) by 2 142 for w = (mpi-wtime) 143 do 144 (when gc 145 (sb-ext:gc :full t)) 146 (multiple-value-bind (l s r) 147 (apply-double-site-operator 148 (svref double (+ s i)) 149 (svref u i) 150 (svref u (1+ i)) 151 (svref v (1+ i)) 152 maxdim) 153 (setf (svref u i) l) 154 (setf (svref u (1+ i)) r) 155 (setf (svref v i) s)) 156 (when subwt 157 (incf (svref subwt i) (- (mpi-wtime) w)))))) 158 159 (defun apply-odd-tebd-operators (length single double submps maxdim &key (subwt nil) (gc nil)) 160 (let* ((s (submps-start submps)) 161 (d (submps-size submps)) 162 (e (+ s d)) 163 (u (submps-tensors submps)) 164 (v (submps-singvals submps))) 165 (if (zerop s) 166 (let ((w (mpi-wtime))) 167 (setf (svref u 0) (apply-single-site-operator (car single) (svref u 0))) 168 (when subwt 169 (incf (svref subwt 0) (- (mpi-wtime) w))))) 170 (if (and (= e length) (evenp e) (not (zerop d))) 171 (let ((w (mpi-wtime))) 172 (setf (svref u (1- d)) (apply-single-site-operator (cadr single) (svref u (1- d)))) 173 (when subwt 174 (incf (svref subwt (1- d)) (- (mpi-wtime) w))))) 175 (loop 176 for i from (if (evenp s) 1 0) below (if (= e length) (1- d) d) by 2 177 for w = (mpi-wtime) 178 do 179 (when gc 180 (sb-ext:gc :full t)) 181 (multiple-value-bind (l s r) 182 (apply-double-site-operator 183 (svref double (+ s i)) 184 (svref u i) (svref u (1+ i)) 185 (svref v (1+ i)) 186 maxdim) 187 (setf (svref u i) l) 188 (setf (svref u (1+ i)) r) 189 (setf (svref v i) s)) 190 (when subwt 191 (incf (svref subwt i) (- (mpi-wtime) w)))))) 192 193 (defun tebd-exchange-phase-0 (rank length submps) 194 (let* ((s (submps-start submps)) 195 (d (submps-size submps)) 196 (e (+ s d)) 197 (u (submps-tensors submps)) 198 (v (submps-singvals submps))) 199 (when (not (zerop d)) 200 (let (l) 201 (sb-sys:without-gcing 202 (when (and (oddp s) (/= s 0)) 203 (push (mpi-issend-object (svref u 0) (1- rank)) l) 204 (push (mpi-issend-object (svref v 0) (1- rank)) l)) 205 (when (and (oddp e) (/= e length)) 206 (setf (svref u d) (mpi-receive-object (1+ rank))) 207 (setf (svref v d) (mpi-receive-object (1+ rank)))) 208 (mapc #'mpi-wait l)))))) 209 210 (defun tebd-exchange-phase-1 (rank length submps) 211 (let* ((s (submps-start submps)) 212 (d (submps-size submps)) 213 (e (+ s d)) 214 (u (submps-tensors submps)) 215 (v (submps-singvals submps))) 216 (when (not (zerop d)) 217 (let (l) 218 (sb-sys:without-gcing 219 (when (and (oddp e) (/= e length)) 220 (push (mpi-issend-object (svref u d) (1+ rank)) l)) 221 (when (and (evenp s) (/= s 0)) 222 (push (mpi-issend-object (svref u 0) (1- rank)) l) 223 (push (mpi-issend-object (svref v 0) (1- rank)) l)) 224 (when (and (oddp s) (/= s 0)) 225 (setf (svref u 0) (mpi-receive-object (1- rank)))) 226 (when (and (evenp e) (/= e length)) 227 (setf (svref u d) (mpi-receive-object (1+ rank))) 228 (setf (svref v d) (mpi-receive-object (1+ rank)))) 229 (mapc #'mpi-wait l)))))) 230 231 (defun tebd-exchange-phase-2 (rank length submps) 232 (let* ((s (submps-start submps)) 233 (d (submps-size submps)) 234 (e (+ s d)) 235 (u (submps-tensors submps)) 236 (v (submps-singvals submps))) 237 (when (not (zerop d)) 238 (let (l) 239 (sb-sys:without-gcing 240 (when (and (evenp e) (/= e length)) 241 (push (mpi-issend-object (svref u d) (1+ rank)) l)) 242 (when (and (oddp s) (/= s 0)) 243 (push (mpi-issend-object (svref u 0) (1- rank)) l) 244 (push (mpi-issend-object (svref v 0) (1- rank)) l)) 245 (when (and (evenp s) (/= s 0)) 246 (setf (svref u 0) (mpi-receive-object (1- rank)))) 247 (when (and (oddp e) (/= e length)) 248 (setf (svref u d) (mpi-receive-object (1+ rank))) 249 (setf (svref v d) (mpi-receive-object (1+ rank)))) 250 (mapc #'mpi-wait l)))))) 251 252 (defun tebd-exchange-phase-3 (rank length submps) 253 (let* ((s (submps-start submps)) 254 (d (submps-size submps)) 255 (e (+ s d)) 256 (u (submps-tensors submps))) 257 (when (not (zerop d)) 258 (let (l) 259 (sb-sys:without-gcing 260 (when (and (oddp e) (/= e length)) 261 (push (mpi-issend-object (svref u d) (1+ rank)) l)) 262 (when (and (oddp s) (/= s 0)) 263 (setf (svref u 0) (mpi-receive-object (1- rank)))) 264 (mapc #'mpi-wait l)))))) 265 266 (defun tebd-evolve (rank length propagator submps maxdim &key (subwt nil) (gc nil)) 267 (tebd-exchange-phase-0 rank length submps) 268 (apply-even-tebd-operators 269 length 270 (propagator-single2 propagator) 271 (propagator-double2 propagator) 272 submps 273 maxdim 274 :subwt subwt 275 :gc gc) 276 (tebd-exchange-phase-1 rank length submps) 277 (apply-odd-tebd-operators 278 length 279 (propagator-single1 propagator) 280 (propagator-double1 propagator) 281 submps 282 maxdim 283 :subwt subwt 284 :gc gc) 285 (tebd-exchange-phase-2 rank length submps) 286 (apply-even-tebd-operators 287 length 288 (propagator-single2 propagator) 289 (propagator-double2 propagator) 290 submps 291 maxdim 292 :subwt subwt 293 :gc gc) 294 (tebd-exchange-phase-3 rank length submps))