Up to 128 qubits: Matlab C++ tutorial
If you are simulating dilute systems with few excitations, simulating states of 128 spins or qubits might be possible. Do not expect to do it with, say, bare Matlab. Instead, you will need to code something in a lower level language that allows you to have greater control over memory that is used. In working with this particular problem I found it advantageous to code it in C++ and embed it in Matlab using the MEX interface.
To illustrate it, below I insert a program I use to simulate a particular problem with spins. It is a subroutine that creates a cell array with M sparse operators. Element (i,j) is a sparse matrix representing a σ- operator on the j-th site of a system with M spins, acting between the spaces with "N" and "N-1" excitations. Using this table or similar ideas you can easily code an XY model acting a space with a fixed number of excitations.
The key, however, is to be able to store such a large number of configurations. For this I use a trick from GCC which is to use 128-bit integers. That's implemented efficiently by the compiler without resorting to special libraries and at most has a 25% performance hit with respect to ordinary integers (64-bit).
/* This is a large integer, large enough to index into a vector * containing all possible states of a spin chain with a fixed * number of excitation. */ typedef mwIndex conf_index; /* This is an even larger integer, capable of holding the configuration * of up to 128 qubits as bits. We use GNU's extension for working * with 128-bit integers. */ typedef unsigned int uint128_t __attribute__((mode(TI))); typedef uint128_t configuration;
The actual program is linked here. You have to compile it and use it more or less as shown below. Note that the same code will teach you how to create sparse matrices, cell arrays and a lot of other stuff.
< M A T L A B (R) >
Copyright 1984-2013 The MathWorks, Inc.
R2013a (8.1.0.604) 64-bit (glnxa64)
February 15, 2013
To get started, type one of these: helpwin, helpdesk, or demo.
For product information, visit www.mathworks.com.
>> mex -largeArrayDims sm_operators_128.cc
Warning: You are using gcc version "4.9.2-6)". The version
currently supported with MEX is "4.4.x".
For a list of currently supported compilers see:
http://www.mathworks.com/support/compilers/current_release/
>> a = sm_operators_128(4,30);
>> a
a =
Columns 1 through 3
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 4 through 6
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 7 through 9
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 10 through 12
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 13 through 15
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 16 through 18
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 19 through 21
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 22 through 24
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 25 through 27
[4060x27405 double] [4060x27405 double] [4060x27405 double]
Columns 28 through 30
[4060x27405 double] [4060x27405 double] [4060x27405 double]
>> size(a)
ans =
1 30
>> a = sm_operators_128(4,70);
>> size(a)
ans =
1 70
Several warnings:
- While useful and extremely fast for building sparse operators, if you want to reach very large problems, you should write a function that applies the operator, not a function that constructs it. More on an upcoming post.
- Beware that this routine does not check memory requirements. Have a look at the size of the Hilbert space, which is roughly the binomial of the # of modes vs # of excitations. If that times 16 bytes looks large in comparison with your memory, try a different alternative.