CUDA-QFT/qft_host.cpp at main · zigpot/CUDA-QFT · GitHub

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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <algorithm>
using std::max;
using std::min;

#include <cuComplex.h>
#include "qft_host.h"


// basic host implementation


// process bits i and i^(1ul << target)
static void phase_shift_host(int control, int target, unsigned long i, cuDoubleComplex *v){
	float phi = 1.0f * float(M_PI) / float(1ul << (control-target));
	cuDoubleComplex phase = {cosf(phi), sinf(phi)};

	unsigned long tgt_bit = (1ul << target);
	unsigned long ctl_bit = (1ul << control);

	if ((i & tgt_bit) == 0) {
		fprintf(stderr, "Logic error: phase_shift_host() should not be called with '(i & tgt_bit) == 0'.\n");
		exit(1);
	}

	if ((i & ctl_bit) != 0) {
		v[i] = cuCmul(v[i], phase);
		// v[i^tgt_bit] stays unchanged
	}
}


// memproses bit i dan i^(1ul << target)
static void hadamard_host(unsigned long target, unsigned long i, cuDoubleComplex *v){
	unsigned long tgt_bit = (1ul << target);

	if ((i & tgt_bit) == 0) {
		fprintf(stderr, "Logic error: phase_shift_host() tidak seharusnya dipanggil dengan '(i & tgt_bit) == 0'.\n");
		exit(1);
	}

	unsigned long i_other = i^tgt_bit;
	cuDoubleComplex ai, aother;	// coefficient i and (i^tgt_bit)
	cuDoubleComplex v_i = v[i];
	cuDoubleComplex v_iother = v[i_other];

	cuDoubleComplex cuM_SQRT1_2 = make_cuDoubleComplex(M_SQRT1_2, 0);
	ai = cuCmul(cuM_SQRT1_2, cuCsub(v_iother, v_i));
	aother = cuCmul(cuM_SQRT1_2, cuCadd(v_iother, v_i));

	v[i] = ai;
	v[i_other] = aother;
}

static void qft_host_kernel(int tgt, unsigned long i, int width, cuDoubleComplex *v){
	unsigned long tgt_bit = (1ul << tgt);

	if ((i & tgt_bit) == 0) {
        fprintf(stderr, "Logic error: qft_host_kernel() should not be called with '(i & tgt_bit) == 0'.\n");
		exit(1);
	}

    // Note: Phase shifts (with the same target) commute.
    //for (int ctl=width-1; ctl>tgt; ctl--) {
	for (int ctl = tgt + 1; ctl < width; ctl++) {
		phase_shift_host(ctl, tgt, i, v);
	}
	hadamard_host(tgt, i, v);
}


// For each qubit, process each state (pair of coefficients) separately.
// This is structured so that it could be implemented in a GPU kernel.
// From here, the next step would be to perform coefficient grouping to achieve
// memory coalescing.
void qft_host(int width, cuDoubleComplex *v){
	unsigned long N = (1ul << width);

    // For each qubit...
	for (int tgt = width - 1; tgt >= 0; tgt--) {
        // For each state...
		for (unsigned long i = 0; i < N; i++) {
			unsigned long tgt_bit = (1ul << tgt);
			if ((i & tgt_bit) != 0) {
				qft_host_kernel(tgt, i, width, v);
			}
		}
	}
}