Parallel Computing in QuSpin

download script

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

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