Parallel computing support

1.2.1 Parallel support in the operator module: hamiltonian, quantum_operator and quantum_LinearOperator

One of the most ubiquitous operations in exact diagonalization codes is the matrix-vector product: the matrix represents a quantum operator and the vector – the quantum state being acted on by the operator. This is used pretty much everywhere except for full diagonalization of the matrix: from evolution to Lanczos methods.

In QuSpin, operators are represented by large sparse matrices and the quantum states are typically represented by dense vectors. In the computer science/mathematics literature, it is well known that this matrix-vector product is one of the most important operations done during a computation, so there has been a lot of work on trying to efficiently implement this operation in parallel. Most of the literature discusses only sparse-matrix – vector product as opposed to sparse-matrix – dense matrix products, and this is reflected in QuSpin’s implementation. Currently QuSpin supports multi-threading only when the multiplication is on a vector (even though multiplication by two-dimensional arrays is allowed as well, but the code switches to a single-threaded version).

We have specifically designed QuSpin to work very efficiently with the structure of the hamiltonian and quantum_operator classes. This low level code replaces the use of SciPy’s default functionality (which adds unnecessary overhead and slows down the code when pushing to large system sizes). This required to limit the number of supported matrix formats used by QuSpin’s operator classes. Currently, we support: csr, csc, dia and also dense matrices when constructing a hamiltonian or quantum_operator objects to allow for a broad range of applicability. For example, one can get a performance boost when transposing your hamiltonian or quantum_operator as csr <–> csc and dia <–> dia without having to copy any data. The dense matrices we fall back on NumPy’s library to do the calculation as it is specifically optimized for the kinds of calculations we need in QuSpin.

For the supported sparse-matrix formats csr, csc, and dia, we have implemented multi-threaded matrix-vector products (see tools.misc.matvec() function), all of which show very nearly linear scaling with increasing the number of cores on modern processors. Even though the performance gains are more modest on older CPU architectures, they can still be useful when simulating large system sizes as one typically needs to allocate a lot of memory space when submitting a job (which usually just means requesting more cores).

To sum up, whenever one can prefer matrix-vector products in the code, using QuSpin’s interface this will lead to the automatic use of multi-threading, when the OpenMP version is used. For instance, one commonly used function, which automatically benefits from multi-threading via the parallel matrix-vector product, is hamiltonian.evolve().

At the same time, in some places automatic multithreading is not so obvious: for instance if one is trying to find the ground state of a particular hamiltonian object H one might do the following:

E,V = H.eigsh(time=t0,k=1,which="SA")

The code just above will actually not use any multi-threading: this is because this code is actually equivilent to doing:

E,V = eigsh(H.tocsr(time=t0),k=1,which="SA")

However, one can still beneft from the multi-threaded matrix-vector product by using the H.aslinearoperator(time=t0) method:

E,V = eigsh(H.aslinearoperator(time=t0),k=1,which="SA")

Casting H as a LinearOperator object enables the use of the methods H.dot() and H.transpose().conj().dot(). These methods will be used to do the eigenvalue calculation, which will then benefit from multi-threading (note that one cannot use LinearOperator by default when calling H.eigsh() since it limits the code functionality).

1.2.2 Parallel support in the general basis classes *_basis_general

Starting from QuSpin 0.3.1, we have efficient implementation of parallel support for the methods in the *_basis_general classes. Additionally, we have also added an implementation of inplace_Op which is used to do ‘on the fly’ calculation of an operator acting on a state using multi-threading OpenMP speed-up (which can be accessed simply by using any general basis in the quantum_LinearOperator class).

Note that the *_basis_1d classes do not support OpenMP.

1.2.3 Parallel support in tools module

The function tools.misc.matvec() wraps an efficient version of various matrix-vector products based on a scheme which provides equal work load to all the threads, regardless of the sparsity structure of the matrix (see this paper for more details). This speedup will be automatically inherited by the function tools.evolution.expm_multiply_parallel(), which creates a more efficient multi-threaded version of SciPy’s SciPy.sparse.linalg.expm_multiply function, and also by the dot attribute of the classes in the Operators module (and hence, for instance, also in the evolve functions).

Notice that the Tools module functions would have to be explicitly used by the user in order for a calculation to gain speedup via OMP.