-
Notifications
You must be signed in to change notification settings - Fork 6
/
heisenberg.cu
245 lines (186 loc) · 7.44 KB
/
heisenberg.cu
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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
#include "hamiltonian.h"
__device__ float HOffBondXHeisenberg(const int si, const int bra, const float JJ)
{
float valH;
//int S0, S1;
//int T0, T1;
valH = JJ*0.5; //contribution from the J part of the Hamiltonian
return valH;
}
__device__ float HOffBondYHeisenberg(const int si, const int bra, const float JJ)
{
float valH;
//int S0, S1;
//int T0, T1;
valH = JJ*0.5; //contribution from the J part of the Hamiltonian
return valH;
}
__device__ float HDiagPartHeisenberg1D(const int bra, int latticeSize, int3* d_Bond, const float JJ)
{
int S0b,S1b ; //spins (bra
int T0,T1; //site
//int P0, P1, P2, P3; //sites for plaquette (Q)
//int s0p, s1p, s2p, s3p;
float valH = 0.f;
for (int Ti=0; Ti<latticeSize; Ti++)
{
//***HEISENBERG PART
T0 = (d_Bond[Ti]).x; //lower left spin
S0b = (bra>>T0)&1;
//if (T0 != Ti) cout<<"Square error 3\n";
T1 = (d_Bond[Ti]).y; //first bond
S1b = (bra>>T1)&1; //unpack bra
valH += JJ*(S0b-0.5)*(S1b-0.5);
}//T0
//cout<<bra<<" "<<valH<<endl;
return valH;
}//HdiagPart
__device__ float HDiagPartHeisenberg2D(const int bra, int latticeSize, int3* d_Bond, const float JJ)
{
int S0b,S1b ; //spins (bra
int T0,T1; //site
//int P0, P1, P2, P3; //sites for plaquette (Q)
//int s0p, s1p, s2p, s3p;
float valH = 0.f;
for (int Ti=0; Ti<latticeSize; Ti++)
{
//***HEISENBERG PART
T0 = (d_Bond[Ti]).x; //lower left spin
S0b = (bra>>T0)&1;
//if (T0 != Ti) cout<<"Square error 3\n";
T1 = (d_Bond[Ti]).y; //first bond
S1b = (bra>>T1)&1; //unpack bra
valH += JJ*(S0b-0.5)*(S1b-0.5);
T1 = (d_Bond[Ti]).z; //second bond
S1b = (bra>>T1)&1; //unpack bra
valH += JJ*(S0b-0.5)*(S1b-0.5);
}//T0
//cout<<bra<<" "<<valH<<endl;
return valH;
}//HdiagPart
__global__ void FillDiagonalsHeisenberg(int* d_basis, f_hamiltonian H, int* d_Bond, parameters data)
{
int row = blockIdx.x*blockDim.x + threadIdx.x;
int latticeSize = data.nsite;
int site = threadIdx.x%(latticeSize);
int dim = H.sectorDim;
unsigned int tempi;
__shared__ int3 tempBond[32];
if (row < dim)
{
tempi = d_basis[row];
(tempBond[site]).x = d_Bond[site];
(tempBond[site]).y = d_Bond[latticeSize + site];
switch( data.dimension )
{
case 1 :
H.vals[row] = HDiagPartHeisenberg1D(tempi, latticeSize, tempBond, data.J1);
break;
case 2 :
(tempBond[site]).z = d_Bond[2*latticeSize + site];
H.vals[row] = HDiagPartHeisenberg2D(tempi, latticeSize, tempBond, data.J1);
break;
}
H.rows[row] = row;
H.cols[row] = row;
H.set[row] = 1;
}
else
{
H.rows[row] = 2*dim;
H.cols[row] = 2*dim;
H.set[row] = 0;
}
}
__global__ void FillSparseHeisenberg(int* d_basisPosition, int* d_basis, f_hamiltonian H, int* d_Bond, parameters data, int offset)
{
int latticeSize = data.nsite;
int dim = H.sectorDim;
int ii = (blockDim.x/(2*latticeSize))*blockIdx.x + threadIdx.x/(2*latticeSize) + offset*(blockDim.x/(2*latticeSize));
int T0 = threadIdx.x%(2*latticeSize);
#if __CUDA_ARCH__ < 200
const int arraySize = 512;
#elif __CUDA_ARCH__ >= 200
const int arraySize = 1024;
#else
#error Could not detect GPU architecture
#endif
__shared__ int3 tempBond[32];
int count;
__shared__ int tempPos[arraySize];
__shared__ float tempVal[arraySize];
//__shared__ uint tempi[arraySize];
unsigned int tempi;
__shared__ uint tempod[arraySize];
int stride = 2*data.dimension*latticeSize;
//int tempcount;
int site = T0%(latticeSize);
count = 0;
int rowTemp;
int start = (bool)(dim%arraySize) ? (dim/arraySize + 1)*arraySize : dim/arraySize;
int s;
int braSector;
//int si, sj;//sk,sl; //spin operators
//unsigned int tempi;// tempod; //tempj;
//cuDoubleComplex tempD;
bool compare;
if( ii < dim )
{
tempi = d_basis[ii];
if (T0 < 2*latticeSize)
{
//Putting bond info in shared memory
(tempBond[site]).x = d_Bond[site];
(tempBond[site]).y = d_Bond[latticeSize + site];
__syncthreads();
//Horizontal bond ---------------
s = (tempBond[site]).x;
tempod[threadIdx.x] = tempi;
braSector = (tempi & (1 << s)) >> s;
tempod[threadIdx.x] ^= (1<<s);
s = (tempBond[site]).y;
braSector ^= (tempi & (1 << s)) >> s;
tempod[threadIdx.x] ^= (1<<s);
//tempod[threadIdx.x] ^= (1<<si); //toggle bit
//tempod[threadIdx.x] ^= (1<<sj); //toggle bit
compare = (d_basisPosition[tempod[threadIdx.x]] != -1) && braSector;
compare &= (d_basisPosition[tempod[threadIdx.x]] > ii);
tempPos[threadIdx.x] = (compare) ? d_basisPosition[tempod[threadIdx.x]] : dim;
tempVal[threadIdx.x] = HOffBondXHeisenberg(site, tempi, data.J1);
count += (int)compare;
//tempcount = (T0/latticeSize);
rowTemp = (T0/latticeSize) ? ii : tempPos[threadIdx.x];
rowTemp = (compare) ? rowTemp : 2*dim;
H.vals[ ii*stride + 4*site + (T0/latticeSize) + start ] = tempVal[threadIdx.x]; // (T0/latticeSize) ? tempVal[threadIdx.x] : cuConj(tempVal[threadIdx.x]);
H.cols[ ii*stride + 4*site + (T0/latticeSize) + start ] = (T0/latticeSize) ? tempPos[threadIdx.x] : ii;
H.rows[ ii*stride + 4*site + (T0/latticeSize) + start ] = rowTemp;
H.set[ ii*stride + 4*site + (T0/latticeSize) + start ] = (int)compare;
if (data.dimension == 2)
{
(tempBond[site]).z = d_Bond[2*latticeSize + site];
//Vertical bond -----------------
s = (tempBond[site]).x;
tempod[threadIdx.x] = tempi;
braSector = (tempi & (1 << s)) >> s;
tempod[threadIdx.x] ^= (1<<s);
s = (tempBond[site]).z;
braSector ^= (tempi & (1 << s)) >> s;
tempod[threadIdx.x] ^= (1<<s);
//tempod[threadIdx.x] ^= (1<<si); //toggle bit
//tempod[threadIdx.x] ^= (1<<sj); //toggle bit
compare = (d_basisPosition[tempod[threadIdx.x]] != -1) && braSector;
compare &= (d_basisPosition[tempod[threadIdx.x]] > ii);
tempPos[threadIdx.x] = (compare) ? d_basisPosition[tempod[threadIdx.x]] : dim;
tempVal[threadIdx.x] = HOffBondYHeisenberg(site,tempi, data.J1);
count += (int)compare;
//tempcount = (T0/latticeSize);
rowTemp = (T0/latticeSize) ? ii : tempPos[threadIdx.x];
rowTemp = (compare) ? rowTemp : 2*dim;
H.vals[ ii*stride + 4*site + 2 + (T0/latticeSize) + start ] = tempVal[threadIdx.x]; // (T0/latticeSize) ? tempVal[threadIdx.x] : cuConj(tempVal[threadIdx.x]);
H.cols[ ii*stride + 4*site + 2 + (T0/latticeSize) + start ] = (T0/latticeSize) ? tempPos[threadIdx.x] : ii;
H.rows[ ii*stride + 4*site + 2 + (T0/latticeSize) + start ] = rowTemp;
H.set[ ii*stride + 4*site + 2 + (T0/latticeSize) + start ] = (int)compare;
}
}
}//end of ii
}//end of FillSparse