GPU accelerated Trotter-Suzuki solver for quantum spin dynamics

This module performs the phase corrections depending on the Z-projection (up or down) and position of each spin in each state of the Ising basis. In this stage there are no cross terms between different basis states (all operations are diagonal), and therefore the parallelization is trivial. However, due to the two-spin terms, the time of evaluation of each phase correction increases quadratically with the number of spins.

Since the Hamiltonian is time-independent, the phase corrections can be computed in advance (i.e. before performing the actual sequence of partial evolutions) and stored in memory for the following evolutions, amortizing its high cost. In fact, Eq. 9 ensures that the phases only depend on the time step and the Hamiltonian parameters. The disadvantage of such pre-computed phase corrections lies on the amount of memory required to store three $2^N$-element vectors, increasing memory usage by 150% and limiting simulations to 27 spins on 6GB GPUs instead of 28. This technique, however, reduces phase correction time for large systems by 97%.

Additionally, since every phase correction is followed by a backward rotation that also operates on the whole state, the phase correction kernel is merged with the backward rotation kernel. Therefore, the phase correction is performed when the rotation kernel reads elements from memory and just before applying the rotation. This saves two global memory operations on each element and a kernel call.

Combined, these improvements reduce the overall simulation time by up to 52%.

The rotation module acts on pairs of basis states by means of the operators $X$ and $Y$ defined in Eq. 5. If a particular basis state has the $j$-th spin down, then it is paired to the basis state that has the $j$-th spin up and the same configuration for the rest of spins. Let $n_{j,\downarrow }$ be the binary representation of a particular basis state in which the $j$-th spin¹¹1Consistently with the binary representation, spins are enumerated from right to left. is down, i.e. its $j$-th bit is $0$: $n_{j,\downarrow }=\dots {\downarrow }_j\dots$. Then, $n_{j,\downarrow }$ is paired to $n_{j,\uparrow }=\dots {\uparrow }_j\dots =n_{j,\downarrow }+2^j$.

This pairing procedure must range the whole Hilbert space. Thus, for every value of $j$ the pairing must be repeated over the $2^{N-1}$ states that have the sought configuration in the $j$-th spin.

In order to show how the parallelization can be performed at this stage, let us exemplify one case of the pairing procedure. If we address the rotation of the second spin in a four-spin system, then the pairs of states under consideration are the following:

Notice that the number of pairings is $8=2^{N-1}$. From this example it is clear that there are no repetitions in any of the basis states when executing one specific rotation (i.e. a particular $j$). Thus, each rotation operation within the same $j$-spin results independent and can be parallelized.

As stated above, the rotation of the $j$-th spin involves $2^{N-1}$ independent operations. Accordingly, a GPU kernel that rotates the $j$-th spin on a single pair is launched on $2^{N-1}$ threads, one per pair. Since the rotation has an extremely low arithmetic intensity of 0.125 operations per byte transferred (i.e. 8 operations per pair of states read and written back to memory) compared to our GPU’s 27.46 double precision FLOPS/bandwidth ratio, the kernel is severely limited by memory bandwidth to feed each thread with data. To avoid this bottleneck, the optimized implementation reads each of these states from global memory and performs as many rotations as possible before writing them back.

Let us consider a subset of $M$ specific spins ($M\le N$) denoted by $F$ and fix its state configuration, i.e. with a particular choice of $\uparrow$ or $\downarrow$ for each of the $M$ spins. Then, the set of all states that satisfy the specified configuration for the spins in $F$ is closed under the pairing operation that flips any spin which is not in $F$. For example, in the case $N=4$ one can consider the $F$ set composed by spins 1 and 3 (i.e. $M=2$). Now, if we fix the configuration $\cdot {\uparrow }_3\cdot {\downarrow }_1$, the set of states that satisfy such configuration is $\{{0100}_b;{0110}_b;{1100}_b;{1110}_b\}$, and this set is closed under the pairing operation for spins 2 and 4 (those which are not in $F$).

The previous observation allows us to launch $2^M$ thread blocks on the GPU, and define their specific spin configurations (the state of the $M$ spins in $F$) using the binary encoding of each thread block’s unique index. This leaves each thread block with its own set of $2^{N-M}$ states to rotate and, by the observation above, every state’s peers for spins $\left\{\mathrm{1..}N\right\}-F$ are also within the set. Each thread block can then copy its set of states to the fast shared memory available in each multiprocessor within a GPU, and perform the rotations for all the spins in $\left\{\mathrm{1..}N\right\}-F$ without any further global memory access. Once these rotations are done, the results are copied back to global memory and the kernel finishes.

While the previous procedure improves rotation time, the approach has a problem: when rotating “higher” spins, the lower bits of the state handled by a thread block are defined by its thread block index, which results in adjacent states being interleaved across different thread blocks. For example, when rotating spins 1 and 2 in a single kernel call on a 4-spin problem, the set of states rotated by each thread block is:

However, when rotating spins 3 and 4, the set of states rotated by each thread block is:

It is easy to notice that consecutive states (e.g. ${0000}_b$ and ${0001}_b$) are assigned to different thread blocks. Since GPU memory controllers fetch long bursts (currently 128 bytes long) of contiguous data to L1 cache, which is local to each multiprocessor (which execute few thread blocks at a time), this interleaving wastes memory bandwidth significantly.

In order to mitigate the mentioned problem, we never include the $s$ lowest spins $\left[\mathrm{1..}s\right]$ in $F$, so that each thread block always operates on runs of at least $2^s$ consecutive elements in memory. Only the first rotation kernel that is called in a full evolution rotates spins $\left[\mathrm{1..}s\right]$, while the rest of the kernels perform $s$ fewer rotations than in the previous approach because these lower spins still consume shared memory. The improved bandwidth efficiency allowed by coalesced memory accesses results in a $300\%$ speedup over the previous kernel on high-spin rotations for $s=3$, which matches the memory controller’s 128-byte line length. Figure 1 shows how a 6-spin rotation is performed using one kernel call with $s=0$ and two kernel calls with $s=3$.

We also tried storing these $2^s$ elements in registers instead of shared memory, to preserve the maximum amount of rotations per kernel call. This generated significant register pressure due to the size of complex double precision numbers, producing low multiprocessor occupancy, and resulted in the kernel performing worse than the kernel that uses shared memory to store every element.

We measured the rotation times for different values of $s$ and different thread block sizes, and found out that there are no optimal values that work for every rotation and every system size. Thus, we wrote a small tool that does a brute force search for the fastest series of rotations that perform a full evolution for every system size that fits in memory. These optimal configurations are stored in a file and are used by our rotation code in following executions.

We replicated the GPU scheme on a CPU using OpenMP intrinsics with few changes to the algorithm: we still precompute phase corrections, divide the system into blocks and perform as many rotations as possible, but each thread loops over whole blocks at a time. The rotation code is instantiated for many different block sizes and values of $s$ using C++ template metaprogramming, providing the compiler with as much compile-time data as possible to perform automatic vectorization. Instead of using shared memory to store intermediate values, since the blocking strategy has strong temporal and spatial locality we assume that the caches on the CPU will hold the values without hitting the memory bus. As in the GPU code, optimal rotation sizes and their tunable parameters are precomputed using a separate brute force search tool.

As a result of this implementation we plot in Fig. 2 the execution times for the original Fortran CPU code without optimizations, the backported CPU code running on one and six threads, and the three different GPU versions. The tests were run on an Intel Core i7-980 six-core processor with triple channel DDR3-1066 memory and a Tesla C2075 GPU. The CPU code was compiled with GCC 4.7.2, using -Ofast -march=native -mtune=native -fno-exceptions -fno-rtti -flto -fopenmp compiler flags, while the GPU code was compiled with CUDA 5.0 using -O -use_fast_math -gencode arch=compute_20,code=sm_20 compiler flags.

In Fig. 2, it can also be observed that for small systems, the original CPU code is faster than parallel implementations because computing the quadratic phase correction is faster than performing a table lookup for small values of $N$, and it does not suffer from parallelization overhead like its CPU siblings. Meanwhile, the GPU versions do not have enough work to feed all the GPU’s execution resources.

On the other hand, when $N$ is considerably large, the execution time of the GPU implementation increases exponentially in $N$. The inflection point around $N\simeq 14$ indicates that the GPU reaches its full capacity, and afterwards each GPU version scales as ${\alpha }^N$, with a factor $\alpha$ higher than $2$. In fact, this parameter can be obtained by a linear fitting of each curve in Fig. 2 for $N>14$. The original CPU implementation yields $\alpha =2.26$, while the CPU Backported (6 thread) yields $\alpha =2.13$. The naive GPU-Basic and the GPU-Optimized result $\alpha =2.13$ and $\alpha =2.05$ respectively. These small differences between the implementations are reflected in the behavior of the relative speedup, which is shown in Fig. 3-a and b for selected implementations.

Since the memory bandwidth is the main limiting factor of the algorithm’s performance, we measured the bandwidth achieved by the multiple rotation kernel to verify the quality of our implementation. Our fastest CUDA kernel implementation uses 80% of the 120GBps reported by NVIDIA’s bandwidth test tool included in the CUDA SDK, while the backported CPU implementation uses between 50% and 70% of the 10 GBps reported by the STREAM benchmark McCalpin2007 () on our system. These results reflect our decision to focus our optimizations on GPU code and ensure that our CPU implementation avoids hitting the memory bus during a series of rotations. We estimate that better tuned CPU implementations can extract a similar percentage of the platform’s memory bandwidth. However, with current high-end dual-socket quad-channel DDR3 platforms reaching 100 GBps XeonBW () and current high-end GPUs doubling this quantity TeslaBW () we expect GPUs to maintain their dominance. The high bandwidth utilization on the GPU also implies that further improvements will require a different algorithm to perform rotations.

The last issue to be addressed concerns the accuracy of our TS implementation. Several comparisons were performed between the present method and exact diagonalization schemes Lea2 (); Lea3 (); Lea4 (), i.e. the solution of the Schrödinger equation as an eigenproblem. The evolution of local and non-local physical magnitudes was compared for one-dimensional spin systems described by Hamiltonians in the form of Eq. 1. For $N=16$ ²²2Standard exact diagonalization is limited by the size of systems that can be handled., a sequence of $5\times {10}^6$ TS steps yields a relative difference bounded by ${10}^{-8}$. A second accuracy test was provided by the spurious total magnetization observed in polarization-conserving Hamiltonians. In fact, this a strictly physical condition: if the initial state has zero total magnetization, then the evolved one should remain so. In general, we observe a linear increase of the total magnetization of the spin system, which is consistent with the general expectancies of a linear increase of the TS error as function of the number of steps. In particular, for $N=22$ and $5\times {10}^6$ TS steps, we observe a relative deviation less than ${10}^{-6}$, a precision which is good enough for practical purposes.

An alternative strategy for error quantification is provided by the Loschmidt echo (LE) scholarpedia (). It evaluates the reversibility of a system’s dynamics in the presence of uncontrolled degrees of freedom Zangara-PRA-2012 (). For a particular initial state $|{\Psi }_0⟩$, the standard LE definition jalpa () is:

where $H_0$ is a reversed Hamiltonian, and $\Sigma$ encloses uncontrolled non-reversed degrees of freedom. In the present case, we set $\Sigma \equiv 0$, which for an ideal perfectly accurate computation would yield $M_{LE}\equiv 1$, $\forall t$. In principle, this statement is completely irrespective of the TS approximation and its intrinsic accuracy (i.e. the validity of Eq. 4). As a matter of fact, exactly the same sequences of TS partial evolutions are applied in the forward and backwards dynamics, except for the change in the sign of $t$. Therefore, the deviations of $M_{LE}$ away from the unity could be intrinsically originated in the execution of floating point operations by the specific Hardware.

We consider initial states given by random superpositions of Ising states (e.g. Eq. 13), for a $N=16$ spin set. Two cases of different dynamical complexity are evaluated. The first case is a ring configuration described by a nearest neighbors Heisenberg Hamiltonian. As above, this interaction conserves total magnetization, and then the dynamics does not explore the whole Hilbert space. For this case, $5\times {10}^5$ TS steps yield $|1-M_{LE}|=5.206647\times {10}^{-8}$, while $5\times {10}^6$ TS steps yield $|1-M_{LE}|=5.2066504\times {10}^{-7}$. Again, the error appears to increase linearly with the number of TS steps. However, we stress here that these deviations cannot be associated to the TS decomposition. The second case is built from the first, adding double quantum terms $S_j^x{S}_k^x-S_j^y{S}_k^y$ up to third next nearest neighbors. In this situation, total magnetization is not conserved, and therefore dynamics effectively mixes all spin projection subspaces exploring the whole Hilbert space. It turns out that $5\times {10}^6$ TS steps yield $|1-M_{LE}|=5.2094531\times {10}^{-7}$, which is slightly larger than the first case. This may indicate that errors in computational operations can depend on the complexity of the many-body dynamics.

Once again, these examples evidence a satisfactory accuracy of our TS-GPU implementation. A systematic study of the LE as an error quantifier is well beyond the scope of this article. This would require ranging over $N$ and the number of TS steps, different Hamiltonian complexities and initial states, among many other factors. However, the LE turns to be a promising witness to address computational errors in GPU and CPU implementations.