Parallel Computing in QuSpin¶
This example shows how to speed up QuSpin code via multi-threading by using
OpenMP,
Intel’s MKL library for NumPy/SciPy (which is Anaconda’s default NumPy version, starting from Anaconda 2.5 onwards).
To install quspin with OpenMP support using anaconda (see also Parallel Computing Support), run
$ conda install -c weinbe58 omp quspin
The example below demonstrates how to use the OpenMP version of quspin for parallel computing. It is set up in such a way that the number of OpenMP and MKL threads is controlled from the command line [cf. code lines 8,9]. To run the script, run
$ python example12.py ${OMP_NUM_THREADS} ${MKL_NUM_THREADS}
You can directly compare the speed for different values of the number of threads [make sure your machine’s processor has more than one core]
$ python example12.py 1 1 # single-threaded computation
$ python example12.py 4 1 # multi-threaded OpenMP computation, speedup for basis functions, evolution, and matrix-vector multiplication
$ python example12.py 1 4 # multi-threaded MKL computation, speedup for diagonalization-like routines
$ python example12.py 4 4 # simulaneous OpenMP and MKL multi-threading speedup
Notice how, as explained in Parallel Computing Support, OMP_NUM_THREADS improves the speed of the basis computation (code line 43), and the time evolution (code line 65), while MKL_NUM_THREADS improves the speed of the exact diagonalization step (code line 60).
Note: there is a common problem with using OpenMP on OSX with anaconda packages for Python 3, which may induce an error unrelated to QuSpin:
$ OMP: Error #15: Initializing libiomp5.dylib, but found libomp.dylib already initialized.
However, this error can be disabled [at one’s own risk!] until it is officially fixed, see code line 7 below.
Script¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 | from __future__ import print_function, division
import sys,os
# line 4 and line 5 below are for development purposes and can be removed
qspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,qspin_path)
#
os.environ['KMP_DUPLICATE_LIB_OK']='True' # uncomment this line if omp error occurs on OSX for python 3
os.environ['OMP_NUM_THREADS']=str(int(sys.argv[1])) # set number of OpenMP threads to run in parallel
os.environ['MKL_NUM_THREADS']=str(int(sys.argv[2])) # set number of MKL threads to run in parallel
#
###########################################################################
# example 12 #
# In this script we show how to use QuSpin's OpenMP and MKL capabilities. #
###########################################################################
from quspin.operators import hamiltonian
from quspin.basis import spin_basis_general
from quspin.operators._make_hamiltonian import _consolidate_static
import numpy as np
from scipy.special import comb
np.random.seed(1) #fixes seed of rng
from time import time # timing package
def run_computation():
#
###### define model parameters ######
J1=1.0 # spin=spin interaction
J2=0.5 # magnetic field strength
Omega=8.0 # drive frequency
Lx, Ly = 4, 4 # linear dimension of spin-1/2 2d lattice
N_2d = Lx*Ly # number of sites for spin-1/2
#
###### setting up user-defined symmetry transformations for 2d lattice ######
sites = np.arange(N_2d) # sites [0,1,2,....]
x = sites%Lx # x positions for sites
y = sites//Lx # y positions for sites
#
T_x = (x+1)%Lx + Lx*y # translation along x-direction
T_y = x +Lx*((y+1)%Ly) # translation along y-direction
#
T_a = (x+1)%Lx + Lx*((y+1)%Ly) # translation along anti-diagonal
T_d = (x-1)%Lx + Lx*((y+1)%Ly) # translation along diagonal
#
###### setting up bases ######
basis_2d = spin_basis_general(N_2d,pauli=False) # making the basis: sped up by OpenMP if symmetries are on
print('finished computing basis')
#
###### setting up hamiltonian ######
# set up time-dependence
def drive(t,Omega):
return np.cos(Omega*t)
drive_args=[Omega,]
# setting up site-coupling lists
J1_list=[[J1,i,T_x[i]] for i in range(N_2d)] + [[J1,i,T_y[i]] for i in range(N_2d)]
J2_list=[[J2,i,T_d[i]] for i in range(N_2d)] + [[J2,i,T_a[i]] for i in range(N_2d)]
#
static =[ ["xx",J1_list],["yy",J1_list],["zz",J1_list] ]
dynamic=[ ["xx",J2_list,drive,drive_args],["yy",J2_list,drive,drive_args],["zz",J2_list,drive,drive_args] ]
# build hamiltonian
H=hamiltonian(static,dynamic,basis=basis_2d,dtype=np.float64,check_symm=False,check_herm=False)
# diagonalise H
E,V=H.eigsh(time=0.0,k=50,which='LA') # H.eigsh sped up by MKL
print('finished computing energies')
psi_0=V[:,0]
# evolve state
t=np.linspace(0.0,20*2*np.pi/Omega,21)
psi_t=H.evolve(psi_0,t[0],t,iterate=True) # H.evolve sped up by OpenMP
for j,psi in enumerate(psi_t):
E_t = H.expt_value(psi,time=t[j])
print("finished evolving up to time step {:d}".format(j) )
# time computation
ti = time() # start timer
run_computation()
print("simulation took {0:.4f} sec".format(time()-ti))
|