Sequential Quantum Gate Decomposer  v1.9.6
Powerful decomposition of general unitarias into one- and two-qubit gates gates
apply_kernel_to_state_vector_input_AVX.cpp
Go to the documentation of this file.
1 /*
2 Created on Fri Jun 26 14:13:26 2020
3 Copyright 2020 Peter Rakyta, Ph.D.
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 
17 @author: Peter Rakyta, Ph.D.
18 */
24 #include "apply_kernel_to_state_vector_input_AVX.h"
25 #include "apply_large_kernel_to_input_AVX.h"
26 #include <immintrin.h>
27 #include "tbb/tbb.h"
28 
29 void apply_2qbit_kernel_to_state_vector_input_AVX_impl(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
30 void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
31 void apply_2qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
32 void apply_3qbit_kernel_to_state_vector_input_AVX_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
33 void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
34 void apply_3qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
35 void apply_4qbit_kernel_to_state_vector_input_AVX_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
36 void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
37 void apply_4qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
38 void apply_5qbit_kernel_to_state_vector_input_AVX_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
39 void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
40 void apply_5qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size);
41 void apply_2qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
42 void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
43 void apply_2qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
44 void apply_3qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
45 void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
46 void apply_3qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
47 void apply_4qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
48 void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
49 void apply_4qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
50 void apply_5qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
51 void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
52 void apply_5qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size);
53 
54 
55 
56 
66 void
67 apply_kernel_to_state_vector_input_AVX(Matrix& u3_1qbit, Matrix& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
68 
69  unsigned int bitmask_low = (1 << target_qbit) - 1;
70  unsigned int bitmask_high = ~bitmask_low;
71 
72  int control_qbit_step_index = (1<<control_qbit);
73 
74  if ( control_qbit == 0 ) {
75 
76  for (int idx=0; idx<matrix_size/2; idx++ ) {
77 
78  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
79  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
80 
81  // the index pair with target qubit state 1
82  int current_idx_pair = current_idx | (1<<target_qbit);
83 
84  if ( current_idx & control_qbit_step_index ) {
85 
86  QGD_Complex16 element = input[current_idx];
87  QGD_Complex16 element_pair = input[current_idx_pair];
88 
89 
90  QGD_Complex16&& tmp1 = mult(u3_1qbit[0], element);
91  QGD_Complex16&& tmp2 = mult(u3_1qbit[1], element_pair);
92 
93  input[current_idx].real = tmp1.real + tmp2.real;
94  input[current_idx].imag = tmp1.imag + tmp2.imag;
95 
96  QGD_Complex16&& tmp3 = mult(u3_1qbit[2], element);
97  QGD_Complex16&& tmp4 = mult(u3_1qbit[3], element_pair);
98 
99  input[current_idx_pair].real = tmp3.real + tmp4.real;
100  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
101 
102 
103 
104  }
105  else if (deriv) {
106  // when calculating derivatives, the constant element should be zeros
107  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
108  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
109  }
110  else {
111  // leave the state as it is
112  continue;
113  }
114 
115 
116  }
117 
118  }
119  else if ( target_qbit == 0 ) {
120 
121 /*
122 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
123 
124 under MIT License
125 
126 Copyright (c) 2018 Qulacs Authors
127 
128 Permission is hereby granted, free of charge, to any person obtaining a copy
129 of this software and associated documentation files (the "Software"), to deal
130 in the Software without restriction, including without limitation the rights
131 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
132 copies of the Software, and to permit persons to whom the Software is
133 furnished to do so, subject to the following conditions:
134 
135 The above copyright notice and this permission notice shall be included in all
136 copies or substantial portions of the Software.
137 */
138 
139  __m256d mv00 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
140  __m256d mv01 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
141  __m256d mv20 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
142  __m256d mv21 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
143 
144  for (int idx=0; idx<matrix_size/2; idx++ ) {
145 
146  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
147  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
148 
149  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
150 
151 
152  double *ptr = (double*)input.get_data() + 2 * current_idx;
153  __m256d data = _mm256_loadu_pd(ptr);
154 
155  __m256d data_u0 = _mm256_mul_pd(data, mv00);
156  __m256d data_u1 = _mm256_mul_pd(data, mv01);
157  __m256d data_u2 = _mm256_hadd_pd(data_u0, data_u1);
158  data_u2 = _mm256_permute4x64_pd(data_u2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
159 
160  __m256d data_d0 = _mm256_mul_pd(data, mv20);
161  __m256d data_d1 = _mm256_mul_pd(data, mv21);
162  __m256d data_d2 = _mm256_hadd_pd(data_d0, data_d1);
163  data_d2 = _mm256_permute4x64_pd(data_d2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
164 
165  __m256d data_r = _mm256_hadd_pd(data_u2, data_d2);
166 
167  data_r = _mm256_permute4x64_pd(data_r, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
168  _mm256_storeu_pd(ptr, data_r);
169 
170  }
171  else if (deriv) {
172  // when calculating derivatives, the constant element should be zeros
173  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex16));
174  }
175  else {
176  // leave the state as it is
177  continue;
178  }
179 
180 
181  }
182 
183  }
184  else {
185 
186 
187 /*
188 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
189 
190 under MIT License
191 
192 Copyright (c) 2018 Qulacs Authors
193 
194 Permission is hereby granted, free of charge, to any person obtaining a copy
195 of this software and associated documentation files (the "Software"), to deal
196 in the Software without restriction, including without limitation the rights
197 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
198 copies of the Software, and to permit persons to whom the Software is
199 furnished to do so, subject to the following conditions:
200 
201 The above copyright notice and this permission notice shall be included in all
202 copies or substantial portions of the Software.
203 */
204 
205  __m256d mv00 = _mm256_set_pd(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
206  __m256d mv01 = _mm256_set_pd( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
207  __m256d mv10 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
208  __m256d mv11 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
209  __m256d mv20 = _mm256_set_pd(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
210  __m256d mv21 = _mm256_set_pd( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
211  __m256d mv30 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
212  __m256d mv31 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
213 
214 
215  for (int idx=0; idx<matrix_size/2; idx=idx+2 ) {
216 
217  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
218  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
219 
220  // the index pair with target qubit state 1
221  int current_idx_pair = current_idx | (1<<target_qbit);
222 
223  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
224 
225 
226  double* element = (double*)input.get_data() + 2 * current_idx;
227  double* element_pair = (double*)input.get_data() + 2 * current_idx_pair;
228 
229 
230  __m256d data0 = _mm256_loadu_pd(element);
231  __m256d data1 = _mm256_loadu_pd(element_pair);
232 
233  __m256d data_u2 = _mm256_mul_pd(data0, mv00);
234  __m256d data_u3 = _mm256_mul_pd(data1, mv10);
235  __m256d data_u4 = _mm256_mul_pd(data0, mv01);
236  __m256d data_u5 = _mm256_mul_pd(data1, mv11);
237 
238  __m256d data_u6 = _mm256_hadd_pd(data_u2, data_u4);
239  __m256d data_u7 = _mm256_hadd_pd(data_u3, data_u5);
240 
241  __m256d data_d2 = _mm256_mul_pd(data0, mv20);
242  __m256d data_d3 = _mm256_mul_pd(data1, mv30);
243  __m256d data_d4 = _mm256_mul_pd(data0, mv21);
244  __m256d data_d5 = _mm256_mul_pd(data1, mv31);
245 
246  __m256d data_d6 = _mm256_hadd_pd(data_d2, data_d4);
247  __m256d data_d7 = _mm256_hadd_pd(data_d3, data_d5);
248 
249  __m256d data_r0 = _mm256_add_pd(data_u6, data_u7);
250  __m256d data_r1 = _mm256_add_pd(data_d6, data_d7);
251 
252  _mm256_storeu_pd(element, data_r0);
253  _mm256_storeu_pd(element_pair, data_r1);
254 
255 
256  }
257  else if (deriv) {
258  // when calculating derivatives, the constant element should be zeros
259  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
260  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
261  }
262  else {
263  // leave the state as it is
264  continue;
265  }
266 
267 
268  //std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
269 
270 
271 
272  }
273 
274 
275 
276  } // else
277 
278 }
279 
280 
281 
282 void
283 apply_kernel_to_state_vector_input_AVX32(Matrix_float& u3_1qbit, Matrix_float& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
284 
285  unsigned int bitmask_low = (1 << target_qbit) - 1;
286  unsigned int bitmask_high = ~bitmask_low;
287 
288  int control_qbit_step_index = (1<<control_qbit);
289 
290  if ( control_qbit == 0 ) {
291 
292  for (int idx=0; idx<matrix_size/2; idx++ ) {
293 
294  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
295  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
296 
297  // the index pair with target qubit state 1
298  int current_idx_pair = current_idx | (1<<target_qbit);
299 
300  if ( current_idx & control_qbit_step_index ) {
301 
302  QGD_Complex8 element = input[current_idx];
303  QGD_Complex8 element_pair = input[current_idx_pair];
304 
305 
306  QGD_Complex8&& tmp1 = mult(u3_1qbit[0], element);
307  QGD_Complex8&& tmp2 = mult(u3_1qbit[1], element_pair);
308 
309  input[current_idx].real = tmp1.real + tmp2.real;
310  input[current_idx].imag = tmp1.imag + tmp2.imag;
311 
312  QGD_Complex8&& tmp3 = mult(u3_1qbit[2], element);
313  QGD_Complex8&& tmp4 = mult(u3_1qbit[3], element_pair);
314 
315  input[current_idx_pair].real = tmp3.real + tmp4.real;
316  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
317 
318 
319 
320  }
321  else if (deriv) {
322  // when calculating derivatives, the constant element should be zeros
323  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
324  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
325  }
326  else {
327  // leave the state as it is
328  continue;
329  }
330 
331 
332  }
333 
334  }
335  else if ( target_qbit == 0 ) {
336 
337 /*
338 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
339 
340 under MIT License
341 
342 Copyright (c) 2018 Qulacs Authors
343 
344 Permission is hereby granted, free of charge, to any person obtaining a copy
345 of this software and associated documentation files (the "Software"), to deal
346 in the Software without restriction, including without limitation the rights
347 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
348 copies of the Software, and to permit persons to whom the Software is
349 furnished to do so, subject to the following conditions:
350 
351 The above copyright notice and this permission notice shall be included in all
352 copies or substantial portions of the Software.
353 */
354 
355  __m256 mv00 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
356  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
357  __m256 mv01 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
358  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
359  __m256 mv20 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
360  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
361  __m256 mv21 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
362  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
363 
364  for (int idx=0; idx<matrix_size/2; idx++ ) {
365 
366  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
367  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
368 
369  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
370 
371 
372  float *ptr = (float*)input.get_data() + 2 * current_idx;
373  __m128 data_low = _mm_loadu_ps(ptr);
374  __m256 data = _mm256_castps128_ps256(data_low);
375  data = _mm256_insertf128_ps(data, data_low, 1);
376 
377  __m256 data_u0 = _mm256_mul_ps(data, mv00);
378  __m256 data_u1 = _mm256_mul_ps(data, mv01);
379  __m256 data_u2 = _mm256_hadd_ps(data_u0, data_u1);
380  data_u2 = _mm256_permutevar8x32_ps(data_u2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
381 
382  __m256 data_d0 = _mm256_mul_ps(data, mv20);
383  __m256 data_d1 = _mm256_mul_ps(data, mv21);
384  __m256 data_d2 = _mm256_hadd_ps(data_d0, data_d1);
385  data_d2 = _mm256_permutevar8x32_ps(data_d2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
386 
387  __m256 data_r = _mm256_hadd_ps(data_u2, data_d2);
388 
389  data_r = _mm256_permutevar8x32_ps(data_r, _mm256_setr_epi32(0, 4, 2, 6, 1, 5, 3, 7));
390  _mm_storeu_ps(ptr, _mm256_castps256_ps128(data_r));
391 
392  }
393  else if (deriv) {
394  // when calculating derivatives, the constant element should be zeros
395  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex8));
396  }
397  else {
398  // leave the state as it is
399  continue;
400  }
401 
402 
403  }
404 
405  }
406  else {
407 
408 
409 /*
410 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
411 
412 under MIT License
413 
414 Copyright (c) 2018 Qulacs Authors
415 
416 Permission is hereby granted, free of charge, to any person obtaining a copy
417 of this software and associated documentation files (the "Software"), to deal
418 in the Software without restriction, including without limitation the rights
419 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
420 copies of the Software, and to permit persons to whom the Software is
421 furnished to do so, subject to the following conditions:
422 
423 The above copyright notice and this permission notice shall be included in all
424 copies or substantial portions of the Software.
425 */
426 
427  __m256 mv00 = _mm256_set_ps(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
428  -u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
429  __m256 mv01 = _mm256_set_ps( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
430  u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
431  __m256 mv10 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real,
432  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
433  __m256 mv11 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag,
434  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
435  __m256 mv20 = _mm256_set_ps(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
436  -u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
437  __m256 mv21 = _mm256_set_ps( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
438  u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
439  __m256 mv30 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real,
440  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
441  __m256 mv31 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag,
442  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
443 
444 
445  if ( target_qbit == 1 ) {
446  // Stride=2 complex32=4 floats; element_pair=element+4 floats.
447  // One 256-bit load covers both; broadcast each 128-bit half to both lanes.
448  for (int idx=0; idx<matrix_size/2; idx=idx+2 ) {
449 
450  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
451  int current_idx_pair = current_idx | (1<<target_qbit);
452 
453  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
454 
455  float* element = (float*)input.get_data() + 2 * current_idx;
456  //float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
457 
458  __m256 data_both = _mm256_loadu_ps(element);
459  __m256 data0 = _mm256_set_m128(_mm256_castps256_ps128(data_both), _mm256_castps256_ps128(data_both));
460  __m256 data1 = _mm256_set_m128(_mm256_extractf128_ps(data_both, 1), _mm256_extractf128_ps(data_both, 1));
461 
462  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
463  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
464  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
465  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
466 
467  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
468  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
469 
470  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
471  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
472  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
473  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
474 
475  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
476  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
477 
478  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
479  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
480 
481  // hadd_ps gives [Re(r0),Re(r1),Im(r0),Im(r1)] per lane; shuffle lower lane
482  // to [Re(r0),Im(r0),Re(r1),Im(r1)]; combined 256-bit store covers both.
483  __m128 r0 = _mm256_castps256_ps128(data_r0);
484  r0 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 1, 2, 0));
485  __m128 r1 = _mm256_castps256_ps128(data_r1);
486  r1 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 1, 2, 0));
487  _mm256_storeu_ps(element, _mm256_set_m128(r1, r0));
488 
489  }
490  else if (deriv) {
491  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
492  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
493  }
494  else { continue; }
495 
496  }
497  }
498  else {
499  // target_qbit >= 2: stride >= 4 complex32 = 8 floats; full 256-bit loads, no overlap.
500  // idx+=4: 4 complex32 per side per iteration — matches f64 register-width scaling.
501  for (int idx=0; idx<matrix_size/2; idx=idx+4 ) {
502 
503  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
504  int current_idx_pair = current_idx | (1<<target_qbit);
505 
506  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
507 
508  float* element = (float*)input.get_data() + 2 * current_idx;
509  float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
510 
511  __m256 data0 = _mm256_loadu_ps(element);
512  __m256 data1 = _mm256_loadu_ps(element_pair);
513 
514  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
515  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
516  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
517  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
518 
519  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
520  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
521 
522  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
523  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
524  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
525  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
526 
527  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
528  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
529 
530  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
531  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
532 
533  // hadd_ps within each 128-bit lane gives [Re(r0),Re(r1),Im(r0),Im(r1)];
534  // _mm256_shuffle_ps deinterleaves both lanes: [Re(r0),Im(r0),Re(r1),Im(r1)].
535  data_r0 = _mm256_shuffle_ps(data_r0, data_r0, _MM_SHUFFLE(3, 1, 2, 0));
536  data_r1 = _mm256_shuffle_ps(data_r1, data_r1, _MM_SHUFFLE(3, 1, 2, 0));
537  _mm256_storeu_ps(element, data_r0);
538  _mm256_storeu_ps(element_pair, data_r1);
539 
540  }
541  else if (deriv) {
542  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
543  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
544  }
545  else { continue; }
546 
547  }
548  }
549 
550 
551 
552  } // else
553 
554 }
555 
556 void
557 apply_kernel_to_state_vector_input_parallel_OpenMP_AVX32(Matrix_float& u3_1qbit, Matrix_float& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
558 
559  unsigned int bitmask_low = (1 << target_qbit) - 1;
560  unsigned int bitmask_high = ~bitmask_low;
561 
562  int control_qbit_step_index = (1<<control_qbit);
563 
564  if ( control_qbit == 0 ) {
565 #pragma omp parallel for
566  for (int idx=0; idx<matrix_size/2; idx++ ) {
567 
568  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
569  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
570 
571  // the index pair with target qubit state 1
572  int current_idx_pair = current_idx | (1<<target_qbit);
573 
574  if ( current_idx & control_qbit_step_index ) {
575 
576  QGD_Complex8 element = input[current_idx];
577  QGD_Complex8 element_pair = input[current_idx_pair];
578 
579 
580  QGD_Complex8&& tmp1 = mult(u3_1qbit[0], element);
581  QGD_Complex8&& tmp2 = mult(u3_1qbit[1], element_pair);
582 
583  input[current_idx].real = tmp1.real + tmp2.real;
584  input[current_idx].imag = tmp1.imag + tmp2.imag;
585 
586  QGD_Complex8&& tmp3 = mult(u3_1qbit[2], element);
587  QGD_Complex8&& tmp4 = mult(u3_1qbit[3], element_pair);
588 
589  input[current_idx_pair].real = tmp3.real + tmp4.real;
590  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
591 
592 
593 
594  }
595  else if (deriv) {
596  // when calculating derivatives, the constant element should be zeros
597  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
598  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
599  }
600  else {
601  // leave the state as it is
602  continue;
603  }
604 
605 
606  }
607 
608  }
609  else if ( target_qbit == 0 ) {
610 
611 /*
612 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
613 
614 under MIT License
615 
616 Copyright (c) 2018 Qulacs Authors
617 
618 Permission is hereby granted, free of charge, to any person obtaining a copy
619 of this software and associated documentation files (the "Software"), to deal
620 in the Software without restriction, including without limitation the rights
621 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
622 copies of the Software, and to permit persons to whom the Software is
623 furnished to do so, subject to the following conditions:
624 
625 The above copyright notice and this permission notice shall be included in all
626 copies or substantial portions of the Software.
627 */
628 
629  __m256 mv00 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
630  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
631  __m256 mv01 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
632  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
633  __m256 mv20 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
634  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
635  __m256 mv21 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
636  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
637 
638 #pragma omp parallel for
639  for (int idx=0; idx<matrix_size/2; idx++ ) {
640 
641  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
642  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
643 
644  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
645 
646 
647  float *ptr = (float*)input.get_data() + 2 * current_idx;
648  __m128 data_low = _mm_loadu_ps(ptr);
649  __m256 data = _mm256_castps128_ps256(data_low);
650  data = _mm256_insertf128_ps(data, data_low, 1);
651 
652  __m256 data_u0 = _mm256_mul_ps(data, mv00);
653  __m256 data_u1 = _mm256_mul_ps(data, mv01);
654  __m256 data_u2 = _mm256_hadd_ps(data_u0, data_u1);
655  data_u2 = _mm256_permutevar8x32_ps(data_u2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
656 
657  __m256 data_d0 = _mm256_mul_ps(data, mv20);
658  __m256 data_d1 = _mm256_mul_ps(data, mv21);
659  __m256 data_d2 = _mm256_hadd_ps(data_d0, data_d1);
660  data_d2 = _mm256_permutevar8x32_ps(data_d2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
661 
662  __m256 data_r = _mm256_hadd_ps(data_u2, data_d2);
663 
664  data_r = _mm256_permutevar8x32_ps(data_r, _mm256_setr_epi32(0, 4, 2, 6, 1, 5, 3, 7));
665  _mm_storeu_ps(ptr, _mm256_castps256_ps128(data_r));
666 
667  }
668  else if (deriv) {
669  // when calculating derivatives, the constant element should be zeros
670  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex8));
671  }
672  else {
673  // leave the state as it is
674  continue;
675  }
676 
677 
678  }
679 
680  }
681  else {
682 
683 
684 /*
685 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
686 
687 under MIT License
688 
689 Copyright (c) 2018 Qulacs Authors
690 
691 Permission is hereby granted, free of charge, to any person obtaining a copy
692 of this software and associated documentation files (the "Software"), to deal
693 in the Software without restriction, including without limitation the rights
694 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
695 copies of the Software, and to permit persons to whom the Software is
696 furnished to do so, subject to the following conditions:
697 
698 The above copyright notice and this permission notice shall be included in all
699 copies or substantial portions of the Software.
700 */
701 
702  __m256 mv00 = _mm256_set_ps(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
703  -u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
704  __m256 mv01 = _mm256_set_ps( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
705  u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
706  __m256 mv10 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real,
707  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
708  __m256 mv11 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag,
709  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
710  __m256 mv20 = _mm256_set_ps(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
711  -u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
712  __m256 mv21 = _mm256_set_ps( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
713  u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
714  __m256 mv30 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real,
715  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
716  __m256 mv31 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag,
717  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
718 
719  if ( target_qbit == 1 ) {
720  // Stride=2 complex32=4 floats; element_pair=element+4 floats.
721  // One 256-bit load covers both; broadcast each 128-bit half to both lanes.
722 #pragma omp parallel for
723  for (int idx=0; idx<matrix_size/2; idx=idx+2 ) {
724 
725  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
726  int current_idx_pair = current_idx | (1<<target_qbit);
727 
728  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
729 
730  float* element = (float*)input.get_data() + 2 * current_idx;
731  //float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
732 
733  __m256 data_both = _mm256_loadu_ps(element);
734  __m256 data0 = _mm256_set_m128(_mm256_castps256_ps128(data_both), _mm256_castps256_ps128(data_both));
735  __m256 data1 = _mm256_set_m128(_mm256_extractf128_ps(data_both, 1), _mm256_extractf128_ps(data_both, 1));
736 
737  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
738  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
739  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
740  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
741 
742  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
743  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
744 
745  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
746  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
747  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
748  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
749 
750  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
751  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
752 
753  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
754  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
755 
756  // hadd_ps gives [Re(r0),Re(r1),Im(r0),Im(r1)] per lane; shuffle lower lane
757  // to [Re(r0),Im(r0),Re(r1),Im(r1)]; combined 256-bit store covers both.
758  __m128 r0 = _mm256_castps256_ps128(data_r0);
759  r0 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 1, 2, 0));
760  __m128 r1 = _mm256_castps256_ps128(data_r1);
761  r1 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 1, 2, 0));
762  _mm256_storeu_ps(element, _mm256_set_m128(r1, r0));
763 
764  }
765  else if (deriv) {
766  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
767  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
768  }
769  else { continue; }
770 
771  }
772  }
773  else {
774  // target_qbit >= 2: stride >= 4 complex32 = 8 floats; full 256-bit loads, no overlap.
775  // idx+=4: 4 complex32 per side per iteration — matches f64 register-width scaling.
776 #pragma omp parallel for
777  for (int idx=0; idx<matrix_size/2; idx=idx+4 ) {
778 
779  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
780  int current_idx_pair = current_idx | (1<<target_qbit);
781 
782  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
783 
784  float* element = (float*)input.get_data() + 2 * current_idx;
785  float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
786 
787  __m256 data0 = _mm256_loadu_ps(element);
788  __m256 data1 = _mm256_loadu_ps(element_pair);
789 
790  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
791  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
792  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
793  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
794 
795  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
796  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
797 
798  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
799  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
800  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
801  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
802 
803  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
804  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
805 
806  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
807  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
808 
809  // hadd_ps within each 128-bit lane gives [Re(r0),Re(r1),Im(r0),Im(r1)];
810  // _mm256_shuffle_ps deinterleaves both lanes: [Re(r0),Im(r0),Re(r1),Im(r1)].
811  data_r0 = _mm256_shuffle_ps(data_r0, data_r0, _MM_SHUFFLE(3, 1, 2, 0));
812  data_r1 = _mm256_shuffle_ps(data_r1, data_r1, _MM_SHUFFLE(3, 1, 2, 0));
813  _mm256_storeu_ps(element, data_r0);
814  _mm256_storeu_ps(element_pair, data_r1);
815 
816  }
817  else if (deriv) {
818  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
819  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
820  }
821  else { continue; }
822 
823  }
824  }
825 
826 
827 
828  } // else
829 
830 }
831 
832 void
833 apply_kernel_to_state_vector_input_parallel_AVX32(Matrix_float& u3_1qbit, Matrix_float& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
834 
835 
836  int grain_size = 64;
837 
838  unsigned int bitmask_low = (1 << target_qbit) - 1;
839  unsigned int bitmask_high = ~bitmask_low;
840 
841  int control_qbit_step_index = (1<<control_qbit);
842 
843  tbb::affinity_partitioner aff_p;
844 
845  if ( control_qbit == 0 ) {
846  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
847 
848  for (int idx=r.begin(); idx<r.end(); idx++) {
849 
850 
851  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
852  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
853 
854  // the index pair with target qubit state 1
855  int current_idx_pair = current_idx | (1<<target_qbit);
856 
857  if ( current_idx & control_qbit_step_index ) {
858 
859  QGD_Complex8 element = input[current_idx];
860  QGD_Complex8 element_pair = input[current_idx_pair];
861 
862 
863  QGD_Complex8&& tmp1 = mult(u3_1qbit[0], element);
864  QGD_Complex8&& tmp2 = mult(u3_1qbit[1], element_pair);
865 
866  input[current_idx].real = tmp1.real + tmp2.real;
867  input[current_idx].imag = tmp1.imag + tmp2.imag;
868 
869  QGD_Complex8&& tmp3 = mult(u3_1qbit[2], element);
870  QGD_Complex8&& tmp4 = mult(u3_1qbit[3], element_pair);
871 
872  input[current_idx_pair].real = tmp3.real + tmp4.real;
873  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
874 
875  }
876  else if (deriv) {
877  // when calculating derivatives, the constant element should be zeros
878  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
879  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
880  }
881  else {
882  // leave the state as it is
883  continue;
884  }
885 
886 
887  //std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
888 
889 
890  }
891  }, aff_p);
892 
893 
894  }
895  else if ( target_qbit == 0 ) {
896 
897 /*
898 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
899 
900 under MIT License
901 
902 Copyright (c) 2018 Qulacs Authors
903 
904 Permission is hereby granted, free of charge, to any person obtaining a copy
905 of this software and associated documentation files (the "Software"), to deal
906 in the Software without restriction, including without limitation the rights
907 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
908 copies of the Software, and to permit persons to whom the Software is
909 furnished to do so, subject to the following conditions:
910 
911 The above copyright notice and this permission notice shall be included in all
912 copies or substantial portions of the Software.
913 */
914 
915  __m256 mv00 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
916  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
917  __m256 mv01 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
918  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
919  __m256 mv20 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
920  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
921  __m256 mv21 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
922  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
923 
924  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
925 
926  for (int idx=r.begin(); idx<r.end(); idx++) {
927 
928  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
929  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
930 
931  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
932 
933 
934  float *ptr = (float*)input.get_data() + 2 * current_idx;
935  __m128 data_low = _mm_loadu_ps(ptr);
936  __m256 data = _mm256_castps128_ps256(data_low);
937  data = _mm256_insertf128_ps(data, data_low, 1);
938 
939  __m256 data_u0 = _mm256_mul_ps(data, mv00);
940  __m256 data_u1 = _mm256_mul_ps(data, mv01);
941  __m256 data_u2 = _mm256_hadd_ps(data_u0, data_u1);
942  data_u2 = _mm256_permutevar8x32_ps(data_u2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
943 
944  __m256 data_d0 = _mm256_mul_ps(data, mv20);
945  __m256 data_d1 = _mm256_mul_ps(data, mv21);
946  __m256 data_d2 = _mm256_hadd_ps(data_d0, data_d1);
947  data_d2 = _mm256_permutevar8x32_ps(data_d2, _mm256_setr_epi32(0, 1, 4, 5, 2, 3, 6, 7));
948 
949  __m256 data_r = _mm256_hadd_ps(data_u2, data_d2);
950 
951  data_r = _mm256_permutevar8x32_ps(data_r, _mm256_setr_epi32(0, 4, 2, 6, 1, 5, 3, 7));
952  _mm_storeu_ps(ptr, _mm256_castps256_ps128(data_r));
953 
954  }
955  else if (deriv) {
956  // when calculating derivatives, the constant element should be zeros
957  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex8));
958  }
959  else {
960  // leave the state as it is
961  continue;
962  }
963 
964 
965  }
966  }, aff_p);
967 
968  }
969  else {
970 
971 /*
972 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
973 
974 under MIT License
975 
976 Copyright (c) 2018 Qulacs Authors
977 
978 Permission is hereby granted, free of charge, to any person obtaining a copy
979 of this software and associated documentation files (the "Software"), to deal
980 in the Software without restriction, including without limitation the rights
981 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
982 copies of the Software, and to permit persons to whom the Software is
983 furnished to do so, subject to the following conditions:
984 
985 The above copyright notice and this permission notice shall be included in all
986 copies or substantial portions of the Software.
987 */
988  __m256 mv00 = _mm256_set_ps(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real,
989  -u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
990  __m256 mv01 = _mm256_set_ps( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag,
991  u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
992  __m256 mv10 = _mm256_set_ps(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real,
993  -u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
994  __m256 mv11 = _mm256_set_ps( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag,
995  u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
996  __m256 mv20 = _mm256_set_ps(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real,
997  -u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
998  __m256 mv21 = _mm256_set_ps( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag,
999  u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
1000  __m256 mv30 = _mm256_set_ps(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real,
1001  -u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
1002  __m256 mv31 = _mm256_set_ps( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag,
1003  u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
1004 
1005 
1006  if ( target_qbit == 1 ) {
1007  // Stride=2 complex32=4 floats; element_pair=element+4 floats.
1008  // One 256-bit load covers both; broadcast each 128-bit half to both lanes.
1009  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
1010 
1011  for (int idx=r.begin(); idx<r.end(); idx=idx+2) {
1012 
1013  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1014  int current_idx_pair = current_idx | (1<<target_qbit);
1015 
1016  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1017 
1018  float* element = (float*)input.get_data() + 2 * current_idx;
1019  //float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
1020 
1021  __m256 data_both = _mm256_loadu_ps(element);
1022  __m256 data0 = _mm256_set_m128(_mm256_castps256_ps128(data_both), _mm256_castps256_ps128(data_both));
1023  __m256 data1 = _mm256_set_m128(_mm256_extractf128_ps(data_both, 1), _mm256_extractf128_ps(data_both, 1));
1024 
1025  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
1026  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
1027  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
1028  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
1029 
1030  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
1031  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
1032 
1033  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
1034  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
1035  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
1036  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
1037 
1038  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
1039  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
1040 
1041  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
1042  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
1043 
1044  // hadd_ps gives [Re(r0),Re(r1),Im(r0),Im(r1)] per lane; shuffle lower lane
1045  // to [Re(r0),Im(r0),Re(r1),Im(r1)]; combined 256-bit store covers both.
1046  __m128 r0 = _mm256_castps256_ps128(data_r0);
1047  r0 = _mm_shuffle_ps(r0, r0, _MM_SHUFFLE(3, 1, 2, 0));
1048  __m128 r1 = _mm256_castps256_ps128(data_r1);
1049  r1 = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(3, 1, 2, 0));
1050  _mm256_storeu_ps(element, _mm256_set_m128(r1, r0));
1051 
1052  }
1053  else if (deriv) {
1054  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
1055  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
1056  }
1057  else { continue; }
1058 
1059  }
1060  }, aff_p);
1061  }
1062  else {
1063  // target_qbit >= 2: stride >= 4 complex32 = 8 floats; full 256-bit loads, no overlap.
1064  // idx+=4: 4 complex32 per side per iteration — matches f64 register-width scaling.
1065  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
1066 
1067  for (int idx=r.begin(); idx<r.end(); idx=idx+4) {
1068 
1069  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1070  int current_idx_pair = current_idx | (1<<target_qbit);
1071 
1072  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1073 
1074  float* element = (float*)input.get_data() + 2 * current_idx;
1075  float* element_pair = (float*)input.get_data() + 2 * current_idx_pair;
1076 
1077  __m256 data0 = _mm256_loadu_ps(element);
1078  __m256 data1 = _mm256_loadu_ps(element_pair);
1079 
1080  __m256 data_u2 = _mm256_mul_ps(data0, mv00);
1081  __m256 data_u3 = _mm256_mul_ps(data1, mv10);
1082  __m256 data_u4 = _mm256_mul_ps(data0, mv01);
1083  __m256 data_u5 = _mm256_mul_ps(data1, mv11);
1084 
1085  __m256 data_u6 = _mm256_hadd_ps(data_u2, data_u4);
1086  __m256 data_u7 = _mm256_hadd_ps(data_u3, data_u5);
1087 
1088  __m256 data_d2 = _mm256_mul_ps(data0, mv20);
1089  __m256 data_d3 = _mm256_mul_ps(data1, mv30);
1090  __m256 data_d4 = _mm256_mul_ps(data0, mv21);
1091  __m256 data_d5 = _mm256_mul_ps(data1, mv31);
1092 
1093  __m256 data_d6 = _mm256_hadd_ps(data_d2, data_d4);
1094  __m256 data_d7 = _mm256_hadd_ps(data_d3, data_d5);
1095 
1096  __m256 data_r0 = _mm256_add_ps(data_u6, data_u7);
1097  __m256 data_r1 = _mm256_add_ps(data_d6, data_d7);
1098 
1099  // hadd_ps within each 128-bit lane gives [Re(r0),Re(r1),Im(r0),Im(r1)];
1100  // _mm256_shuffle_ps deinterleaves both lanes: [Re(r0),Im(r0),Re(r1),Im(r1)].
1101  data_r0 = _mm256_shuffle_ps(data_r0, data_r0, _MM_SHUFFLE(3, 1, 2, 0));
1102  data_r1 = _mm256_shuffle_ps(data_r1, data_r1, _MM_SHUFFLE(3, 1, 2, 0));
1103  _mm256_storeu_ps(element, data_r0);
1104  _mm256_storeu_ps(element_pair, data_r1);
1105 
1106  }
1107  else if (deriv) {
1108  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex8));
1109  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex8));
1110  }
1111  else { continue; }
1112 
1113  }
1114  }, aff_p);
1115  }
1116 
1117 
1118  } // else
1119 
1120 
1121 
1122 
1123 }
1124 
1125 void apply_2qbit_kernel_to_state_vector_input_AVX(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1126  apply_2qbit_kernel_to_state_vector_input_AVX_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1127 }
1128 
1129 void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1130  apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1131 }
1132 
1133 void apply_2qbit_kernel_to_state_vector_input_AVX_TBB(Matrix& two_qbit_unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1134  apply_2qbit_kernel_to_state_vector_input_AVX_TBB_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1135 }
1136 
1137 void apply_3qbit_kernel_to_state_vector_input_AVX(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1138  apply_3qbit_kernel_to_state_vector_input_AVX_impl(unitary, input, std::move(involved_qbits), matrix_size);
1139 }
1140 
1141 void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1142  apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(unitary, input, std::move(involved_qbits), matrix_size);
1143 }
1144 
1145 void apply_3qbit_kernel_to_state_vector_input_AVX_TBB(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1146  apply_3qbit_kernel_to_state_vector_input_AVX_TBB_impl(unitary, input, std::move(involved_qbits), matrix_size);
1147 }
1148 
1149 void apply_4qbit_kernel_to_state_vector_input_AVX(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1150  apply_4qbit_kernel_to_state_vector_input_AVX_impl(unitary, input, std::move(involved_qbits), matrix_size);
1151 }
1152 
1153 void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1154  apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(unitary, input, std::move(involved_qbits), matrix_size);
1155 }
1156 
1157 void apply_4qbit_kernel_to_state_vector_input_AVX_TBB(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1158  apply_4qbit_kernel_to_state_vector_input_AVX_TBB_impl(unitary, input, std::move(involved_qbits), matrix_size);
1159 }
1160 
1161 void apply_5qbit_kernel_to_state_vector_input_AVX(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1162  apply_5qbit_kernel_to_state_vector_input_AVX_impl(unitary, input, std::move(involved_qbits), matrix_size);
1163 }
1164 
1165 void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1166  apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(unitary, input, std::move(involved_qbits), matrix_size);
1167 }
1168 
1169 void apply_5qbit_kernel_to_state_vector_input_AVX_TBB(Matrix& unitary, Matrix& input, std::vector<int> involved_qbits, const int& matrix_size) {
1170  apply_5qbit_kernel_to_state_vector_input_AVX_TBB_impl(unitary, input, std::move(involved_qbits), matrix_size);
1171 }
1172 
1173 void apply_2qbit_kernel_to_state_vector_input_AVX32(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1174  apply_2qbit_kernel_to_state_vector_input_AVX32_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1175 }
1176 
1177 void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1178  apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1179 }
1180 
1181 void apply_2qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float& two_qbit_unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1182  apply_2qbit_kernel_to_state_vector_input_AVX_TBB32_impl(two_qbit_unitary, input, std::move(involved_qbits), matrix_size);
1183 }
1184 
1185 void apply_3qbit_kernel_to_state_vector_input_AVX32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1186  apply_3qbit_kernel_to_state_vector_input_AVX32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1187 }
1188 
1189 void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1190  apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1191 }
1192 
1193 void apply_3qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1194  apply_3qbit_kernel_to_state_vector_input_AVX_TBB32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1195 }
1196 
1197 void apply_4qbit_kernel_to_state_vector_input_AVX32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1198  apply_4qbit_kernel_to_state_vector_input_AVX32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1199 }
1200 
1201 void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1202  apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1203 }
1204 
1205 void apply_4qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1206  apply_4qbit_kernel_to_state_vector_input_AVX_TBB32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1207 }
1208 
1209 void apply_5qbit_kernel_to_state_vector_input_AVX32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1210  apply_5qbit_kernel_to_state_vector_input_AVX32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1211 }
1212 
1213 void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1214  apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1215 }
1216 
1217 void apply_5qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float& unitary, Matrix_float& input, std::vector<int> involved_qbits, const int& matrix_size) {
1218  apply_5qbit_kernel_to_state_vector_input_AVX_TBB32_impl(unitary, input, std::move(involved_qbits), matrix_size);
1219 }
1220 
1230 void
1231 apply_kernel_to_state_vector_input_parallel_OpenMP_AVX(Matrix& u3_1qbit, Matrix& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
1232 
1233  unsigned int bitmask_low = (1 << target_qbit) - 1;
1234  unsigned int bitmask_high = ~bitmask_low;
1235 
1236  int control_qbit_step_index = (1<<control_qbit);
1237 
1238  if ( control_qbit == 0 ) {
1239 #pragma omp parallel for
1240  for (int idx=0; idx<matrix_size/2; idx++ ) {
1241 
1242  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1243  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1244 
1245  // the index pair with target qubit state 1
1246  int current_idx_pair = current_idx | (1<<target_qbit);
1247 
1248  if ( current_idx & control_qbit_step_index ) {
1249 
1250  QGD_Complex16 element = input[current_idx];
1251  QGD_Complex16 element_pair = input[current_idx_pair];
1252 
1253 
1254  QGD_Complex16&& tmp1 = mult(u3_1qbit[0], element);
1255  QGD_Complex16&& tmp2 = mult(u3_1qbit[1], element_pair);
1256 
1257  input[current_idx].real = tmp1.real + tmp2.real;
1258  input[current_idx].imag = tmp1.imag + tmp2.imag;
1259 
1260  QGD_Complex16&& tmp3 = mult(u3_1qbit[2], element);
1261  QGD_Complex16&& tmp4 = mult(u3_1qbit[3], element_pair);
1262 
1263  input[current_idx_pair].real = tmp3.real + tmp4.real;
1264  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
1265 
1266 
1267 
1268  }
1269  else if (deriv) {
1270  // when calculating derivatives, the constant element should be zeros
1271  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
1272  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
1273  }
1274  else {
1275  // leave the state as it is
1276  continue;
1277  }
1278 
1279 
1280  }
1281 
1282  }
1283  else if ( target_qbit == 0 ) {
1284 
1285 /*
1286 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
1287 
1288 under MIT License
1289 
1290 Copyright (c) 2018 Qulacs Authors
1291 
1292 Permission is hereby granted, free of charge, to any person obtaining a copy
1293 of this software and associated documentation files (the "Software"), to deal
1294 in the Software without restriction, including without limitation the rights
1295 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
1296 copies of the Software, and to permit persons to whom the Software is
1297 furnished to do so, subject to the following conditions:
1298 
1299 The above copyright notice and this permission notice shall be included in all
1300 copies or substantial portions of the Software.
1301 */
1302 
1303  __m256d mv00 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
1304  __m256d mv01 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
1305  __m256d mv20 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
1306  __m256d mv21 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
1307 
1308 #pragma omp parallel for
1309  for (int idx=0; idx<matrix_size/2; idx++ ) {
1310 
1311  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1312  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1313 
1314  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1315 
1316 
1317  double *ptr = (double*)input.get_data() + 2 * current_idx;
1318  __m256d data = _mm256_loadu_pd(ptr);
1319 
1320  __m256d data_u0 = _mm256_mul_pd(data, mv00);
1321  __m256d data_u1 = _mm256_mul_pd(data, mv01);
1322  __m256d data_u2 = _mm256_hadd_pd(data_u0, data_u1);
1323  data_u2 = _mm256_permute4x64_pd(data_u2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1324 
1325  __m256d data_d0 = _mm256_mul_pd(data, mv20);
1326  __m256d data_d1 = _mm256_mul_pd(data, mv21);
1327  __m256d data_d2 = _mm256_hadd_pd(data_d0, data_d1);
1328  data_d2 = _mm256_permute4x64_pd(data_d2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1329 
1330  __m256d data_r = _mm256_hadd_pd(data_u2, data_d2);
1331 
1332  data_r = _mm256_permute4x64_pd(data_r, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1333  _mm256_storeu_pd(ptr, data_r);
1334 
1335  }
1336  else if (deriv) {
1337  // when calculating derivatives, the constant element should be zeros
1338  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex16));
1339  }
1340  else {
1341  // leave the state as it is
1342  continue;
1343  }
1344 
1345 
1346  }
1347 
1348  }
1349  else {
1350 
1351 
1352 /*
1353 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
1354 
1355 under MIT License
1356 
1357 Copyright (c) 2018 Qulacs Authors
1358 
1359 Permission is hereby granted, free of charge, to any person obtaining a copy
1360 of this software and associated documentation files (the "Software"), to deal
1361 in the Software without restriction, including without limitation the rights
1362 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
1363 copies of the Software, and to permit persons to whom the Software is
1364 furnished to do so, subject to the following conditions:
1365 
1366 The above copyright notice and this permission notice shall be included in all
1367 copies or substantial portions of the Software.
1368 */
1369 
1370  __m256d mv00 = _mm256_set_pd(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
1371  __m256d mv01 = _mm256_set_pd( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
1372  __m256d mv10 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
1373  __m256d mv11 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
1374  __m256d mv20 = _mm256_set_pd(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
1375  __m256d mv21 = _mm256_set_pd( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
1376  __m256d mv30 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
1377  __m256d mv31 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
1378 
1379 #pragma omp parallel for
1380  for (int idx=0; idx<matrix_size/2; idx=idx+2 ) {
1381 
1382  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1383  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1384 
1385  // the index pair with target qubit state 1
1386  int current_idx_pair = current_idx | (1<<target_qbit);
1387 
1388  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1389 
1390 
1391  double* element = (double*)input.get_data() + 2 * current_idx;
1392  double* element_pair = (double*)input.get_data() + 2 * current_idx_pair;
1393 
1394 
1395  __m256d data0 = _mm256_loadu_pd(element);
1396  __m256d data1 = _mm256_loadu_pd(element_pair);
1397 
1398  __m256d data_u2 = _mm256_mul_pd(data0, mv00);
1399  __m256d data_u3 = _mm256_mul_pd(data1, mv10);
1400  __m256d data_u4 = _mm256_mul_pd(data0, mv01);
1401  __m256d data_u5 = _mm256_mul_pd(data1, mv11);
1402 
1403  __m256d data_u6 = _mm256_hadd_pd(data_u2, data_u4);
1404  __m256d data_u7 = _mm256_hadd_pd(data_u3, data_u5);
1405 
1406  __m256d data_d2 = _mm256_mul_pd(data0, mv20);
1407  __m256d data_d3 = _mm256_mul_pd(data1, mv30);
1408  __m256d data_d4 = _mm256_mul_pd(data0, mv21);
1409  __m256d data_d5 = _mm256_mul_pd(data1, mv31);
1410 
1411  __m256d data_d6 = _mm256_hadd_pd(data_d2, data_d4);
1412  __m256d data_d7 = _mm256_hadd_pd(data_d3, data_d5);
1413 
1414  __m256d data_r0 = _mm256_add_pd(data_u6, data_u7);
1415  __m256d data_r1 = _mm256_add_pd(data_d6, data_d7);
1416 
1417  _mm256_storeu_pd(element, data_r0);
1418  _mm256_storeu_pd(element_pair, data_r1);
1419 
1420 
1421  }
1422  else if (deriv) {
1423  // when calculating derivatives, the constant element should be zeros
1424  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
1425  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
1426  }
1427  else {
1428  // leave the state as it is
1429  continue;
1430  }
1431 
1432 
1433  //std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
1434 
1435 
1436 
1437  }
1438 
1439 
1440 
1441  } // else
1442 
1443 }
1444 
1445 
1446 
1447 
1448 
1458 void
1459 apply_kernel_to_state_vector_input_parallel_AVX(Matrix& u3_1qbit, Matrix& input, const bool& deriv, const int& target_qbit, const int& control_qbit, const int& matrix_size) {
1460 
1461 
1462  int grain_size = 64;
1463 
1464  unsigned int bitmask_low = (1 << target_qbit) - 1;
1465  unsigned int bitmask_high = ~bitmask_low;
1466 
1467  int control_qbit_step_index = (1<<control_qbit);
1468 
1469  tbb::affinity_partitioner aff_p;
1470 
1471  if ( control_qbit == 0 ) {
1472  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
1473 
1474  for (int idx=r.begin(); idx<r.end(); idx++) {
1475 
1476 
1477  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1478  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1479 
1480  // the index pair with target qubit state 1
1481  int current_idx_pair = current_idx | (1<<target_qbit);
1482 
1483  if ( current_idx & control_qbit_step_index ) {
1484 
1485  QGD_Complex16 element = input[current_idx];
1486  QGD_Complex16 element_pair = input[current_idx_pair];
1487 
1488 
1489  QGD_Complex16&& tmp1 = mult(u3_1qbit[0], element);
1490  QGD_Complex16&& tmp2 = mult(u3_1qbit[1], element_pair);
1491 
1492  input[current_idx].real = tmp1.real + tmp2.real;
1493  input[current_idx].imag = tmp1.imag + tmp2.imag;
1494 
1495  QGD_Complex16&& tmp3 = mult(u3_1qbit[2], element);
1496  QGD_Complex16&& tmp4 = mult(u3_1qbit[3], element_pair);
1497 
1498  input[current_idx_pair].real = tmp3.real + tmp4.real;
1499  input[current_idx_pair].imag = tmp3.imag + tmp4.imag;
1500 
1501  }
1502  else if (deriv) {
1503  // when calculating derivatives, the constant element should be zeros
1504  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
1505  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
1506  }
1507  else {
1508  // leave the state as it is
1509  continue;
1510  }
1511 
1512 
1513  //std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
1514 
1515 
1516  }
1517  }, aff_p);
1518 
1519 
1520  }
1521  else if ( target_qbit == 0 ) {
1522 
1523 /*
1524 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
1525 
1526 under MIT License
1527 
1528 Copyright (c) 2018 Qulacs Authors
1529 
1530 Permission is hereby granted, free of charge, to any person obtaining a copy
1531 of this software and associated documentation files (the "Software"), to deal
1532 in the Software without restriction, including without limitation the rights
1533 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
1534 copies of the Software, and to permit persons to whom the Software is
1535 furnished to do so, subject to the following conditions:
1536 
1537 The above copyright notice and this permission notice shall be included in all
1538 copies or substantial portions of the Software.
1539 */
1540 
1541  __m256d mv00 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
1542  __m256d mv01 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
1543  __m256d mv20 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
1544  __m256d mv21 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
1545 
1546  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
1547 
1548  for (int idx=r.begin(); idx<r.end(); idx++) {
1549 
1550  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1551  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1552 
1553  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1554 
1555 
1556  double *ptr = (double*)input.get_data() + 2 * current_idx;
1557  __m256d data = _mm256_loadu_pd(ptr);
1558 
1559  __m256d data_u0 = _mm256_mul_pd(data, mv00);
1560  __m256d data_u1 = _mm256_mul_pd(data, mv01);
1561  __m256d data_u2 = _mm256_hadd_pd(data_u0, data_u1);
1562  data_u2 = _mm256_permute4x64_pd(data_u2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1563 
1564  __m256d data_d0 = _mm256_mul_pd(data, mv20);
1565  __m256d data_d1 = _mm256_mul_pd(data, mv21);
1566  __m256d data_d2 = _mm256_hadd_pd(data_d0, data_d1);
1567  data_d2 = _mm256_permute4x64_pd(data_d2, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1568 
1569  __m256d data_r = _mm256_hadd_pd(data_u2, data_d2);
1570 
1571  data_r = _mm256_permute4x64_pd(data_r, 216); // (3210) -> (3120) : 1*0 + 4*2 + 16*1 + 64*3 = 216
1572  _mm256_storeu_pd(ptr, data_r);
1573 
1574  }
1575  else if (deriv) {
1576  // when calculating derivatives, the constant element should be zeros
1577  memset(input.get_data() + current_idx, 0, input.cols * 2 *sizeof(QGD_Complex16));
1578  }
1579  else {
1580  // leave the state as it is
1581  continue;
1582  }
1583 
1584 
1585  }
1586  }, aff_p);
1587 
1588  }
1589  else {
1590 
1591 /*
1592 AVX kernel developed according to https://github.com/qulacs/qulacs/blob/main/src/csim/update_ops_matrix_dense_single.cpp
1593 
1594 under MIT License
1595 
1596 Copyright (c) 2018 Qulacs Authors
1597 
1598 Permission is hereby granted, free of charge, to any person obtaining a copy
1599 of this software and associated documentation files (the "Software"), to deal
1600 in the Software without restriction, including without limitation the rights
1601 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
1602 copies of the Software, and to permit persons to whom the Software is
1603 furnished to do so, subject to the following conditions:
1604 
1605 The above copyright notice and this permission notice shall be included in all
1606 copies or substantial portions of the Software.
1607 */
1608  __m256d mv00 = _mm256_set_pd(-u3_1qbit[0].imag, u3_1qbit[0].real, -u3_1qbit[0].imag, u3_1qbit[0].real);
1609  __m256d mv01 = _mm256_set_pd( u3_1qbit[0].real, u3_1qbit[0].imag, u3_1qbit[0].real, u3_1qbit[0].imag);
1610  __m256d mv10 = _mm256_set_pd(-u3_1qbit[1].imag, u3_1qbit[1].real, -u3_1qbit[1].imag, u3_1qbit[1].real);
1611  __m256d mv11 = _mm256_set_pd( u3_1qbit[1].real, u3_1qbit[1].imag, u3_1qbit[1].real, u3_1qbit[1].imag);
1612  __m256d mv20 = _mm256_set_pd(-u3_1qbit[2].imag, u3_1qbit[2].real, -u3_1qbit[2].imag, u3_1qbit[2].real);
1613  __m256d mv21 = _mm256_set_pd( u3_1qbit[2].real, u3_1qbit[2].imag, u3_1qbit[2].real, u3_1qbit[2].imag);
1614  __m256d mv30 = _mm256_set_pd(-u3_1qbit[3].imag, u3_1qbit[3].real, -u3_1qbit[3].imag, u3_1qbit[3].real);
1615  __m256d mv31 = _mm256_set_pd( u3_1qbit[3].real, u3_1qbit[3].imag, u3_1qbit[3].real, u3_1qbit[3].imag);
1616 
1617 
1618  tbb::parallel_for( tbb::blocked_range<int>(0,matrix_size/2,grain_size), [&](tbb::blocked_range<int> r) {
1619 
1620  for (int idx=r.begin(); idx<r.end(); idx=idx+2) {
1621  // generate index by inserting state 0 into the place of the target qbit while pushing high bits left by one
1622  int current_idx = ((idx & bitmask_high) << 1) | (idx & bitmask_low);
1623 
1624  // the index pair with target qubit state 1
1625  int current_idx_pair = current_idx | (1<<target_qbit);
1626 
1627  if (control_qbit < 0 || (current_idx & control_qbit_step_index) ) {
1628 
1629 
1630  double* element = (double*)input.get_data() + 2 * current_idx;
1631  double* element_pair = (double*)input.get_data() + 2 * current_idx_pair;
1632 
1633 
1634  __m256d data0 = _mm256_loadu_pd(element);
1635  __m256d data1 = _mm256_loadu_pd(element_pair);
1636 
1637  __m256d data_u2 = _mm256_mul_pd(data0, mv00);
1638  __m256d data_u3 = _mm256_mul_pd(data1, mv10);
1639  __m256d data_u4 = _mm256_mul_pd(data0, mv01);
1640  __m256d data_u5 = _mm256_mul_pd(data1, mv11);
1641 
1642  __m256d data_u6 = _mm256_hadd_pd(data_u2, data_u4);
1643  __m256d data_u7 = _mm256_hadd_pd(data_u3, data_u5);
1644 
1645  __m256d data_d2 = _mm256_mul_pd(data0, mv20);
1646  __m256d data_d3 = _mm256_mul_pd(data1, mv30);
1647  __m256d data_d4 = _mm256_mul_pd(data0, mv21);
1648  __m256d data_d5 = _mm256_mul_pd(data1, mv31);
1649 
1650  __m256d data_d6 = _mm256_hadd_pd(data_d2, data_d4);
1651  __m256d data_d7 = _mm256_hadd_pd(data_d3, data_d5);
1652 
1653  __m256d data_r0 = _mm256_add_pd(data_u6, data_u7);
1654  __m256d data_r1 = _mm256_add_pd(data_d6, data_d7);
1655 
1656  _mm256_storeu_pd(element, data_r0);
1657  _mm256_storeu_pd(element_pair, data_r1);
1658 
1659 
1660  }
1661  else if (deriv) {
1662  // when calculating derivatives, the constant element should be zeros
1663  memset(input.get_data() + current_idx, 0, input.cols * sizeof(QGD_Complex16));
1664  memset(input.get_data() + current_idx_pair, 0, input.cols * sizeof(QGD_Complex16));
1665  }
1666  else {
1667  // leave the state as it is
1668  continue;
1669  }
1670 
1671 
1672  //std::cout << current_idx_target << " " << current_idx_target_pair << std::endl;
1673 
1674 
1675 
1676  }
1677  }, aff_p);
1678 
1679 
1680 
1681  } // else
1682 
1683 
1684 
1685 
1686 }
1687 
1688 
1689 
1690 
1691 
1692 
1693 
apply_2qbit_kernel_to_state_vector_input_AVX
void apply_2qbit_kernel_to_state_vector_input_AVX(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1125
apply_kernel_to_state_vector_input_AVX.h
apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32
void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1213
apply_3qbit_kernel_to_state_vector_input_AVX32
void apply_3qbit_kernel_to_state_vector_input_AVX32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1185
apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32
void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1177
tests.decomposition.test_optmization_problem_combined.real
real
Definition: test_optmization_problem_combined.py:67
apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl
void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5992
apply_3qbit_kernel_to_state_vector_input_AVX_TBB_impl
void apply_3qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply three-qubit gate kernel to a state vector using AVX optimization and TBB parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:2268
apply_5qbit_kernel_to_state_vector_input_AVX
void apply_5qbit_kernel_to_state_vector_input_AVX(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1161
apply_3qbit_kernel_to_state_vector_input_AVX_TBB
void apply_3qbit_kernel_to_state_vector_input_AVX_TBB(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1145
apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP
void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1141
apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP
void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1153
apply_5qbit_kernel_to_state_vector_input_AVX_TBB_impl
void apply_5qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply five-qubit gate kernel to a state vector using AVX optimization and TBB parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:3179
apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP_impl
void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply five-qubit gate kernel to a state vector using AVX optimization and OpenMP parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:2842
apply_4qbit_kernel_to_state_vector_input_AVX
void apply_4qbit_kernel_to_state_vector_input_AVX(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1149
apply_4qbit_kernel_to_state_vector_input_AVX_TBB32
void apply_4qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1205
QGD_Complex8::real
float real
real part
Definition: QGDTypes.h:47
apply_4qbit_kernel_to_state_vector_input_AVX32
void apply_4qbit_kernel_to_state_vector_input_AVX32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1197
QGD_Complex8
Structure type representing single-precision complex numbers.
Definition: QGDTypes.h:46
apply_2qbit_kernel_to_state_vector_input_AVX_TBB_impl
void apply_2qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply two-qubit gate kernel to a state vector using AVX optimization and TBB parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:1825
apply_kernel_to_state_vector_input_parallel_OpenMP_AVX
void apply_kernel_to_state_vector_input_parallel_OpenMP_AVX(Matrix &u3_1qbit, Matrix &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
Parallel AVX kernel on a state vector (parallelized with OpenMP)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1231
apply_kernel_to_state_vector_input_AVX
void apply_kernel_to_state_vector_input_AVX(Matrix &u3_1qbit, Matrix &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
AVX kernel on a state vector.
Definition: apply_kernel_to_state_vector_input_AVX.cpp:67
apply_2qbit_kernel_to_state_vector_input_AVX32
void apply_2qbit_kernel_to_state_vector_input_AVX32(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1173
example.data
data
load the unitary from file
Definition: example.py:51
apply_kernel_to_state_vector_input_AVX32
void apply_kernel_to_state_vector_input_AVX32(Matrix_float &u3_1qbit, Matrix_float &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:283
mult
QGD_Complex16 mult(QGD_Complex16 &a, QGD_Complex16 &b)
Call to calculate the product of two complex scalars.
Definition: common.cpp:298
apply_kernel_to_state_vector_input_parallel_AVX32
void apply_kernel_to_state_vector_input_parallel_AVX32(Matrix_float &u3_1qbit, Matrix_float &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:833
QGD_Complex8::imag
float imag
imaginary part
Definition: QGDTypes.h:48
apply_kernel_to_state_vector_input_parallel_AVX
void apply_kernel_to_state_vector_input_parallel_AVX(Matrix &u3_1qbit, Matrix &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
Parallel AVX kernel on a state vector (parallelized with Intel TBB)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1459
apply_4qbit_kernel_to_state_vector_input_AVX_TBB
void apply_4qbit_kernel_to_state_vector_input_AVX_TBB(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1157
apply_5qbit_kernel_to_state_vector_input_AVX_impl
void apply_5qbit_kernel_to_state_vector_input_AVX_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply five-qubit gate kernel to a state vector using AVX optimization.
Definition: apply_large_kernel_to_input_AVX.cpp:1216
matrix_base::get_data
scalar * get_data() const
Call to get the pointer to the stored data.
Definition: matrix_base.hpp:285
apply_2qbit_kernel_to_state_vector_input_AVX32_impl
void apply_2qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5948
apply_4qbit_kernel_to_state_vector_input_AVX_impl
void apply_4qbit_kernel_to_state_vector_input_AVX_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply four-qubit gate kernel to a state vector using AVX optimization.
Definition: apply_large_kernel_to_input_AVX.cpp:1013
apply_2qbit_kernel_to_state_vector_input_AVX_impl
void apply_2qbit_kernel_to_state_vector_input_AVX_impl(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply two-qubit gate kernel to a state vector using AVX optimization.
Definition: apply_large_kernel_to_input_AVX.cpp:516
apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl
void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5952
apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP_impl
void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply four-qubit gate kernel to a state vector using AVX optimization and OpenMP parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:2443
apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32
void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1201
matrix_base::cols
int cols
The number of columns.
Definition: matrix_base.hpp:44
example.matrix_size
matrix_size
[load Umtx]
Definition: example.py:58
apply_4qbit_kernel_to_state_vector_input_AVX32_impl
void apply_4qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5976
apply_5qbit_kernel_to_state_vector_input_AVX_TBB
void apply_5qbit_kernel_to_state_vector_input_AVX_TBB(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1169
apply_5qbit_kernel_to_state_vector_input_AVX32_impl
void apply_5qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5988
apply_3qbit_kernel_to_state_vector_input_AVX32_impl
void apply_3qbit_kernel_to_state_vector_input_AVX32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5964
apply_5qbit_kernel_to_state_vector_input_AVX_TBB32
void apply_5qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1217
QC_sim_benchmark.target_qbit
target_qbit
Definition: QC_sim_benchmark.py:86
example_get_circuit_unitary.unitary
unitary
Definition: example_get_circuit_unitary.py:105
QGD_Complex16
Structure type representing complex numbers in the SQUANDER package.
Definition: QGDTypes.h:38
apply_2qbit_kernel_to_state_vector_input_AVX_TBB32_impl
void apply_2qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5956
apply_kernel_to_state_vector_input_parallel_OpenMP_AVX32
void apply_kernel_to_state_vector_input_parallel_OpenMP_AVX32(Matrix_float &u3_1qbit, Matrix_float &input, const bool &deriv, const int &target_qbit, const int &control_qbit, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:557
Matrix
Double-precision complex matrix (float64).
Definition: matrix.h:38
apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl
void apply_4qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5980
apply_5qbit_kernel_to_state_vector_input_AVX_TBB32_impl
void apply_5qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5996
apply_large_kernel_to_input_AVX.h
Header file for AVX-optimized implementations for applying multi-qubit gate kernels to quantum state ...
Matrix_float
Single-precision complex matrix (float32).
Definition: matrix_float.h:41
apply_2qbit_kernel_to_state_vector_input_AVX_TBB32
void apply_2qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float &two_qbit_unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1181
apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP
void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1129
apply_2qbit_kernel_to_state_vector_input_AVX_TBB
void apply_2qbit_kernel_to_state_vector_input_AVX_TBB(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1133
apply_3qbit_kernel_to_state_vector_input_AVX_TBB32_impl
void apply_3qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5972
apply_4qbit_kernel_to_state_vector_input_AVX_TBB_impl
void apply_4qbit_kernel_to_state_vector_input_AVX_TBB_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply four-qubit gate kernel to a state vector using AVX optimization and TBB parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:2640
QGD_Complex16::real
double real
the real part of a complex number
Definition: QGDTypes.h:40
apply_3qbit_kernel_to_state_vector_input_AVX_impl
void apply_3qbit_kernel_to_state_vector_input_AVX_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply a 3-qubit quantum gate (unitary matrix) to a state vector using AVX intrinsics.
Definition: apply_large_kernel_to_input_AVX.cpp:849
apply_3qbit_kernel_to_state_vector_input_AVX_TBB32
void apply_3qbit_kernel_to_state_vector_input_AVX_TBB32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1193
apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32
void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1189
apply_4qbit_kernel_to_state_vector_input_AVX_TBB32_impl
void apply_4qbit_kernel_to_state_vector_input_AVX_TBB32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5984
apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP_impl
void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply three-qubit gate kernel to a state vector using AVX optimization and OpenMP parallelization...
Definition: apply_large_kernel_to_input_AVX.cpp:2095
apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP_impl
void apply_2qbit_kernel_to_state_vector_input_AVX_OpenMP_impl(Matrix &two_qbit_unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Apply two-qubit gate kernel to a state vector using AVX optimization and OpenMP parallelization.
Definition: apply_large_kernel_to_input_AVX.cpp:1553
apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP
void apply_5qbit_kernel_to_state_vector_input_AVX_OpenMP(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1165
apply_3qbit_kernel_to_state_vector_input_AVX
void apply_3qbit_kernel_to_state_vector_input_AVX(Matrix &unitary, Matrix &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1137
apply_5qbit_kernel_to_state_vector_input_AVX32
void apply_5qbit_kernel_to_state_vector_input_AVX32(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_kernel_to_state_vector_input_AVX.cpp:1209
QGD_Complex16::imag
double imag
the imaginary part of a complex number
Definition: QGDTypes.h:42
QC_sim_benchmark.control_qbit
control_qbit
Definition: QC_sim_benchmark.py:86
apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl
void apply_3qbit_kernel_to_state_vector_input_AVX_OpenMP32_impl(Matrix_float &unitary, Matrix_float &input, std::vector< int > involved_qbits, const int &matrix_size)
Definition: apply_large_kernel_to_input_AVX.cpp:5968
Generated on Thu Jul 9 2026 07:59:23 for Sequential Quantum Gate Decomposer by   doxygen 1.8.13