Is this page helpful?

Code example (expectation value)#

The following code example illustrates how to define a tensor network state, a tensor network operator, and then compute the expectation value of the tensor network operator over the tensor network state. The full code can be found in the NVIDIA/cuQuantum repository (here).

Headers and error handling#

 7#include <cstdlib>
 8#include <cstdio>
 9#include <cassert>
10#include <complex>
11#include <vector>
12#include <bitset>
13#include <iostream>
14
15#include <cuda_runtime.h>
16#include <cutensornet.h>
17
18
19#define HANDLE_CUDA_ERROR(x) \
20{ const auto err = x; \
21  if( err != cudaSuccess ) \
22  { printf("CUDA error %s in line %d\n", cudaGetErrorString(err), __LINE__); fflush(stdout); std::abort(); } \
23};
24
25#define HANDLE_CUTN_ERROR(x) \
26{ const auto err = x; \
27  if( err != CUTENSORNET_STATUS_SUCCESS ) \
28  { printf("cuTensorNet error %s in line %d\n", cutensornetGetErrorString(err), __LINE__); fflush(stdout); std::abort(); } \
29};
30
31
32int main()
33{
34  static_assert(sizeof(size_t) == sizeof(int64_t), "Please build this sample on a 64-bit architecture!");
35
36  constexpr std::size_t fp64size = sizeof(double);

Define the tensor network state#

Let’s define a tensor network state corresponding to a 16-qubit quantum circuit.

40  // Quantum state configuration
41  constexpr int32_t numQubits = 16; // number of qubits
42  const std::vector<int64_t> qubitDims(numQubits,2); // qubit dimensions
43  std::cout << "Quantum circuit: " << numQubits << " qubits\n";

Initialize the cuTensorNet library handle#

47  // Initialize the cuTensorNet library
48  HANDLE_CUDA_ERROR(cudaSetDevice(0));
49  cutensornetHandle_t cutnHandle;
50  HANDLE_CUTN_ERROR(cutensornetCreate(&cutnHandle));
51  std::cout << "Initialized cuTensorNet library on GPU 0\n";

Define quantum gates on GPU#

55  // Define necessary quantum gate tensors in Host memory
56  const double invsq2 = 1.0 / std::sqrt(2.0);
57  //  Hadamard gate
58  const std::vector<std::complex<double>> h_gateH {{invsq2, 0.0},  {invsq2, 0.0},
59                                                   {invsq2, 0.0}, {-invsq2, 0.0}};
60  //  Pauli X gate
61  const std::vector<std::complex<double>> h_gateX {{0.0, 0.0}, {1.0, 0.0},
62                                                   {1.0, 0.0}, {0.0, 0.0}};
63  //  Pauli Y gate
64  const std::vector<std::complex<double>> h_gateY {{0.0, 0.0}, {0.0, -1.0},
65                                                   {0.0, 1.0}, {0.0, 0.0}};
66  //  Pauli Z gate
67  const std::vector<std::complex<double>> h_gateZ {{1.0, 0.0}, {0.0, 0.0},
68                                                   {0.0, 0.0}, {-1.0, 0.0}};
69  //  CX gate
70  const std::vector<std::complex<double>> h_gateCX {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
71                                                    {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
72                                                    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0},
73                                                    {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}};
74
75  // Copy quantum gates to Device memory
76  void *d_gateH{nullptr}, *d_gateX{nullptr}, *d_gateY{nullptr}, *d_gateZ{nullptr}, *d_gateCX{nullptr};
77  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateH, 4 * (2 * fp64size)));
78  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateX, 4 * (2 * fp64size)));
79  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateY, 4 * (2 * fp64size)));
80  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateZ, 4 * (2 * fp64size)));
81  HANDLE_CUDA_ERROR(cudaMalloc(&d_gateCX, 16 * (2 * fp64size)));
82  std::cout << "Allocated quantum gate memory on GPU\n";
83  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateH, h_gateH.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
84  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateX, h_gateX.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
85  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateY, h_gateY.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
86  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateZ, h_gateZ.data(), 4 * (2 * fp64size), cudaMemcpyHostToDevice));
87  HANDLE_CUDA_ERROR(cudaMemcpy(d_gateCX, h_gateCX.data(), 16 * (2 * fp64size), cudaMemcpyHostToDevice));
88  std::cout << "Copied quantum gates to GPU memory\n";

Allocate the scratch buffer on GPU#

 92  // Query the free memory on Device
 93  std::size_t freeSize{0}, totalSize{0};
 94  HANDLE_CUDA_ERROR(cudaMemGetInfo(&freeSize, &totalSize));
 95  const std::size_t scratchSize = (freeSize - (freeSize % 4096)) / 2; // use half of available memory with alignment
 96  void *d_scratch{nullptr};
 97  void *h_scratch{nullptr};  // host scratch workspace buffer (attached if the prepare call reports a need)
 98  int64_t hostWorksize{0};
 99  HANDLE_CUDA_ERROR(cudaMalloc(&d_scratch, scratchSize));
100  std::cout << "Allocated " << scratchSize << " bytes of scratch memory on GPU\n";

Create a pure tensor network state#

Now let’s create a pure tensor network state for a 16-qubit quantum circuit.

104  // Create the initial quantum state
105  cutensornetState_t quantumState;
106  HANDLE_CUTN_ERROR(cutensornetCreateState(cutnHandle, CUTENSORNET_STATE_PURITY_PURE, numQubits, qubitDims.data(),
107                    CUDA_C_64F, &quantumState));
108  std::cout << "Created the initial quantum state\n";

Apply quantum gates#

Let’s construct the GHZ quantum circuit by applying the corresponding quantum gates.

112  // Construct the final quantum circuit state (apply quantum gates) for the GHZ circuit
113  int64_t id;
114  HANDLE_CUTN_ERROR(cutensornetStateApplyTensorOperator(cutnHandle, quantumState, 1, std::vector<int32_t>{{0}}.data(),
115                    d_gateH, nullptr, 1, 0, 1, &id));
116  for(int32_t i = 1; i < numQubits; ++i) {
117    HANDLE_CUTN_ERROR(cutensornetStateApplyTensorOperator(cutnHandle, quantumState, 2, std::vector<int32_t>{{i-1,i}}.data(),
118                      d_gateCX, nullptr, 1, 0, 1, &id));
119  }
120  std::cout << "Applied quantum gates\n";

Construct a tensor network operator#

Now let’s create an empty tensor network operator for 16-qubit states and then append three components to it, where each component is a direct product of Pauli matrices scaled by some complex coefficient (like in the Jordan-Wigner representation).

124  // Create an empty tensor network operator
125  cutensornetNetworkOperator_t hamiltonian;
126  HANDLE_CUTN_ERROR(cutensornetCreateNetworkOperator(cutnHandle, numQubits, qubitDims.data(), CUDA_C_64F, &hamiltonian));
127  // Append component (0.5 * Z1 * Z2) to the tensor network operator
128  {
129    const int32_t numModes[] = {1, 1}; // Z1 acts on 1 mode, Z2 acts on 1 mode
130    const int32_t modesZ1[] = {1}; // state modes Z1 acts on
131    const int32_t modesZ2[] = {2}; // state modes Z2 acts on
132    const int32_t * stateModes[] = {modesZ1, modesZ2}; // state modes (Z1 * Z2) acts on
133    const void * gateData[] = {d_gateZ, d_gateZ}; // GPU pointers to gate data
134    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.5,0.0},
135                      2, numModes, stateModes, NULL, gateData, &id));
136  }
137  // Append component (0.25 * Y3) to the tensor network operator
138  {
139    const int32_t numModes[] = {1}; // Y3 acts on 1 mode
140    const int32_t modesY3[] = {3}; // state modes Y3 acts on
141    const int32_t * stateModes[] = {modesY3}; // state modes (Y3) acts on
142    const void * gateData[] = {d_gateY}; // GPU pointers to gate data
143    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.25,0.0},
144                      1, numModes, stateModes, NULL, gateData, &id));
145  }
146  // Append component (0.13 * Y0 X2 Z3) to the tensor network operator
147  {
148    const int32_t numModes[] = {1, 1, 1}; // Y0 acts on 1 mode, X2 acts on 1 mode, Z3 acts on 1 mode
149    const int32_t modesY0[] = {0}; // state modes Y0 acts on
150    const int32_t modesX2[] = {2}; // state modes X2 acts on
151    const int32_t modesZ3[] = {3}; // state modes Z3 acts on
152    const int32_t * stateModes[] = {modesY0, modesX2, modesZ3}; // state modes (Y0 * X2 * Z3) acts on
153    const void * gateData[] = {d_gateY, d_gateX, d_gateZ}; // GPU pointers to gate data
154    HANDLE_CUTN_ERROR(cutensornetNetworkOperatorAppendProduct(cutnHandle, hamiltonian, cuDoubleComplex{0.13,0.0},
155                      3, numModes, stateModes, NULL, gateData, &id));
156  }
157  std::cout << "Constructed a tensor network operator: (0.5 * Z1 * Z2) + (0.25 * Y3) + (0.13 * Y0 * X2 * Z3)" << std::endl;

Create the expectation value object#

Once the quantum circuit and the tensor network operator have been constructed, let’s create the expectation value object that will compute the expectation value of the specified tensor network operator over the final state of the specified quantum circuit.

161  // Specify the quantum circuit expectation value
162  cutensornetStateExpectation_t expectation;
163  HANDLE_CUTN_ERROR(cutensornetCreateExpectation(cutnHandle, quantumState, hamiltonian, &expectation));
164  std::cout << "Created the specified quantum circuit expectation value\n";

Configure the expectation value calculation#

Now we can configure the expectation value object by setting the number of hyper-samples to be used by the tensor network contraction path finder.

168  // Configure the computation of the specified quantum circuit expectation value
169  const int32_t numHyperSamples = 8; // desired number of hyper samples used in the tensor network contraction path finder
170  HANDLE_CUTN_ERROR(cutensornetExpectationConfigure(cutnHandle, expectation,
171                    CUTENSORNET_EXPECTATION_CONFIG_NUM_HYPER_SAMPLES, &numHyperSamples, sizeof(numHyperSamples)));

Prepare the expectation value calculation#

Let’s create a workspace descriptor and prepare the computation of the desired expectation value.

175  // Prepare the specified quantum circuit expectation value for computation
176  cutensornetWorkspaceDescriptor_t workDesc;
177  HANDLE_CUTN_ERROR(cutensornetCreateWorkspaceDescriptor(cutnHandle, &workDesc));
178  std::cout << "Created the workspace descriptor\n";
179  HANDLE_CUTN_ERROR(cutensornetExpectationPrepare(cutnHandle, expectation, scratchSize, workDesc, 0x0));
180  std::cout << "Prepared the specified quantum circuit expectation value\n";
181  double flops {0.0};
182  HANDLE_CUTN_ERROR(cutensornetExpectationGetInfo(cutnHandle, expectation,
183                    CUTENSORNET_EXPECTATION_INFO_FLOPS, &flops, sizeof(flops)));
184  std::cout << "Total flop count = " << (flops/1e9) << " GFlop\n";

Set up the workspace#

Now we can set up the required workspace buffer.

188  // Attach the workspace buffer
189  int64_t worksize {0};
190  HANDLE_CUTN_ERROR(cutensornetWorkspaceGetMemorySize(cutnHandle,
191                                                      workDesc,
192                                                      CUTENSORNET_WORKSIZE_PREF_RECOMMENDED,
193                                                      CUTENSORNET_MEMSPACE_DEVICE,
194                                                      CUTENSORNET_WORKSPACE_SCRATCH,
195                                                      &worksize));
196  std::cout << "Required scratch GPU workspace size (bytes) = " << worksize << std::endl;
197  if(worksize <= scratchSize) {
198    HANDLE_CUTN_ERROR(cutensornetWorkspaceSetMemory(cutnHandle, workDesc, CUTENSORNET_MEMSPACE_DEVICE,
199                      CUTENSORNET_WORKSPACE_SCRATCH, d_scratch, worksize));
200  }else{
201    std::cout << "ERROR: Insufficient workspace size on Device!\n";
202    std::abort();
203  }
204
205  // Query the host scratch workspace size and attach a host scratch buffer if one is required
206  HANDLE_CUTN_ERROR(cutensornetWorkspaceGetMemorySize(cutnHandle, workDesc, CUTENSORNET_WORKSIZE_PREF_RECOMMENDED,
207                                                      CUTENSORNET_MEMSPACE_HOST, CUTENSORNET_WORKSPACE_SCRATCH, &hostWorksize));
208  std::cout << "Scratch Host workspace size (bytes) = " << hostWorksize << std::endl;
209  if(hostWorksize > 0) {
210    h_scratch = std::realloc(h_scratch, hostWorksize);
211    if(h_scratch == nullptr) { std::cout << "ERROR: Host scratch allocation failed!\n"; std::abort(); }
212    HANDLE_CUTN_ERROR(cutensornetWorkspaceSetMemory(cutnHandle, workDesc, CUTENSORNET_MEMSPACE_HOST,
213                      CUTENSORNET_WORKSPACE_SCRATCH, h_scratch, hostWorksize));
214  }
215  std::cout << "Set the workspace buffer\n";

Compute the requested expectation value#

Once everything has been set up, we compute the requested expectation value and print it. Note that the returned expectation value is not normalized. The 2-norm of the tensor network state is returned as a separate argument.

219  // Compute the specified quantum circuit expectation value
220  std::complex<double> expectVal{0.0,0.0}, stateNorm2{0.0,0.0};
221  HANDLE_CUTN_ERROR(cutensornetExpectationCompute(cutnHandle, expectation, workDesc,
222                    static_cast<void*>(&expectVal), static_cast<void*>(&stateNorm2), 0x0));
223  std::cout << "Computed the specified quantum circuit expectation value\n";
224  expectVal /= stateNorm2;
225  std::cout << "Expectation value = (" << expectVal.real() << ", " << expectVal.imag() << ")\n";
226  std::cout << "Squared 2-norm of the state = (" << stateNorm2.real() << ", " << stateNorm2.imag() << ")\n";

Free resources#

230  // Destroy the workspace descriptor
231  HANDLE_CUTN_ERROR(cutensornetDestroyWorkspaceDescriptor(workDesc));
232  std::cout << "Destroyed the workspace descriptor\n";
233
234  // Destroy the quantum circuit expectation value
235  HANDLE_CUTN_ERROR(cutensornetDestroyExpectation(expectation));
236  std::cout << "Destroyed the quantum circuit state expectation value\n";
237
238  // Destroy the tensor network operator
239  HANDLE_CUTN_ERROR(cutensornetDestroyNetworkOperator(hamiltonian));
240  std::cout << "Destroyed the tensor network operator\n";
241
242  // Destroy the quantum circuit state
243  HANDLE_CUTN_ERROR(cutensornetDestroyState(quantumState));
244  std::cout << "Destroyed the quantum circuit state\n";
245
246  std::free(h_scratch);
247  HANDLE_CUDA_ERROR(cudaFree(d_scratch));
248  HANDLE_CUDA_ERROR(cudaFree(d_gateCX));
249  HANDLE_CUDA_ERROR(cudaFree(d_gateZ));
250  HANDLE_CUDA_ERROR(cudaFree(d_gateY));
251  HANDLE_CUDA_ERROR(cudaFree(d_gateX));
252  HANDLE_CUDA_ERROR(cudaFree(d_gateH));
253  std::cout << "Freed memory on GPU\n";
254
255  // Finalize the cuTensorNet library
256  HANDLE_CUTN_ERROR(cutensornetDestroy(cutnHandle));
257  std::cout << "Finalized the cuTensorNet library\n";
258
259  return 0;
260}