diff --git a/src/simulators/statevector/qubitvector.hpp b/src/simulators/statevector/qubitvector.hpp index 575983f58c..fab03580bd 100755 --- a/src/simulators/statevector/qubitvector.hpp +++ b/src/simulators/statevector/qubitvector.hpp @@ -33,6 +33,7 @@ using indexes_t = std::unique_ptr; using complex_t = std::complex; using cvector_t = std::vector; using rvector_t = std::vector; +template using areg_t = std::array; //============================================================================ // BIT MASKS and indexing @@ -154,9 +155,6 @@ class QubitVector { // Set all entries in the vector to 0. void zero(); - // Return as an int an N qubit bitstring for M-N qubit bit k with 0s inserted - // for N qubits at the locations specified by qubits_sorted. - // qubits_sorted must be sorted lowest to highest. Eg. {0, 1}. // index0 returns only the integer representation of a number of bits set // to zero inserted into an arbitrary bit string. // Eg: for qubits 0,2 in a state k = ba ( ba = 00 => k=0, etc). @@ -164,7 +162,15 @@ class QubitVector { // indexes0([1,3], k) -> int(0b0a) // Example: k = 77 = 1001101 , qubits_sorted = [1,4] // ==> output = 297 = 100101001 (with 0's put into places 1 and 4). - uint_t index0(const reg_t &qubits_sorted, const uint_t k) const; + template + uint_t index0(const list_t &qubits_sorted, const uint_t k) const; + + // Return a vector of of 2^N in ints, where each int corresponds to an N qubit + // bitstring for M-N qubit bits in state k, and the specified N qubits in states + // [0, ..., 2^N - 1] qubits_sorted must be sorted lowest to highest. Eg. {0, 1}. + // qubits specifies the location of the qubits in the retured strings. + template + areg_t<1ULL << N> indexes(const areg_t &qs, const areg_t &qubits_sorted, const uint_t k) const; // Return a std::unique_ptr to an array of of 2^N in ints // each int corresponds to an N qubit bitstring for M-N qubit bits in state k, @@ -419,20 +425,13 @@ class QubitVector { template void apply_lambda(Lambda&& func); - // Apply a single-qubit lambda function to all blocks of the statevector - // for the given qubit. The function signature should be: - // [&](const int_t k1, const int_t k2)->void - // where (k1, k2) are the 0 and 1 indexes for each qubit block - template - void apply_lambda(Lambda&& func, const uint_t qubit); - // Apply a N-qubit lambda function to all blocks of the statevector // for the given qubits. The function signature should be: // [&](const indexes_t &inds)->void // where inds are the 2 ** N indexes for each N-qubit block returned by // the dynamic_indexes function - template - void apply_lambda(Lambda&& func, const reg_t &qubits); + template + void apply_lambda(Lambda&& func, const list_t &qubits); //----------------------------------------------------------------------- // State matrix update with Lambda functions @@ -440,24 +439,14 @@ class QubitVector { // These functions loop through the indexes of the qubitvector data and // apply a lambda function taking a vector matrix argument to each entry. - // Apply a single-qubit matrix lambda function to all blocks of the - // statevector for the given qubit. The function signature should be: - // [&](const int_t k1, const int_t k2, const cvector_t &m)->void - // where (k1, k2) are the 0 and 1 indexes for each qubit block and - // m is a vectorized complex matrix. - template - void apply_matrix_lambda(Lambda&& func, - const uint_t qubit, - const cvector_t &mat); - // Apply a N-qubit matrix lambda function to all blocks of the statevector // for the given qubits. The function signature should be: // [&](const indexes_t &inds, const cvector_t &m)->void // where inds are the 2 ** N indexes for each N-qubit block returned by // the dynamic_indexes function and m is a vectorized complex matrix. - template + template void apply_matrix_lambda(Lambda&& func, - const reg_t &qubits, + const list_t &qubits, const cvector_t &mat); //----------------------------------------------------------------------- @@ -712,17 +701,34 @@ cvector_t QubitVector::vector() const { //------------------------------------------------------------------------------ template -uint_t QubitVector::index0(const reg_t& qubits_sorted, const uint_t k) const { +template +uint_t QubitVector::index0(const list_t &qubits_sorted, const uint_t k) const { uint_t lowbits, retval = k; - for (const auto& qubit : qubits_sorted) { - lowbits = retval & MASKS[qubit]; - retval >>= qubit; - retval <<= qubit + 1; + for (size_t j = 0; j < qubits_sorted.size(); j++) { + lowbits = retval & MASKS[qubits_sorted[j]]; + retval >>= qubits_sorted[j]; + retval <<= qubits_sorted[j] + 1; retval |= lowbits; } return retval; } +template +template +areg_t<1ULL << N> QubitVector::indexes(const areg_t &qs, + const areg_t &qubits_sorted, + const uint_t k) const { + areg_t<1ULL << N> ret; + ret[0] = index0(qubits_sorted, k); + for (size_t i = 0; i < N; i++) { + const auto n = 1ULL << i; + const auto bit = BITS[qs[i]]; + for (size_t j = 0; j < n; j++) + ret[n + j] = ret[j] | bit; + } + return ret; +} + template indexes_t QubitVector::indexes(const reg_t& qubits, const reg_t& qubits_sorted, @@ -926,42 +932,15 @@ void QubitVector::apply_lambda(Lambda&& func) { } } -// Single qubit -template -template -void QubitVector::apply_lambda(Lambda&& func, - const uint_t qubit) { - // Error checking - #ifdef DEBUG - check_qubit(qubit); - #endif - - const int_t END1 = data_size_; // end for k1 loop - const int_t END2 = BITS[qubit]; // end for k2 loop - const int_t STEP1 = BITS[qubit + 1]; // step for k1 loop -#pragma omp parallel if (num_qubits_ > omp_threshold_ && omp_threads_ > 1) num_threads(omp_threads_) - { -#ifdef _WIN32 - #pragma omp for -#else - #pragma omp for collapse(2) -#endif - for (int_t k1 = 0; k1 < END1; k1 += STEP1) - for (int_t k2 = 0; k2 < END2; k2++) { - std::forward(func)(k1, k2); - } - } -} - // Dynamic N-qubit template -template +template void QubitVector::apply_lambda(Lambda&& func, - const reg_t &qubits) { + const list_t &qubits) { // Error checking #ifdef DEBUG - for (const auto &qubit : qs) + for (const auto &qubit : qubits) check_qubit(qubit); #endif @@ -983,39 +962,10 @@ void QubitVector::apply_lambda(Lambda&& func, //------------------------------------------------------------------------------ // Matrix Lambda //------------------------------------------------------------------------------ - -// Single qubit -template -template -void QubitVector::apply_matrix_lambda(Lambda&& func, - const uint_t qubit, - const cvector_t &mat) { - // Error checking - #ifdef DEBUG - check_qubit(qubit); - #endif - - const int_t END1 = data_size_; // end for k1 loop - const int_t END2 = BITS[qubit]; // end for k2 loop - const int_t STEP1 = BITS[qubit + 1]; // step for k1 loop -#pragma omp parallel if (num_qubits_ > omp_threshold_ && omp_threads_ > 1) num_threads(omp_threads_) - { -#ifdef _WIN32 - #pragma omp for -#else - #pragma omp for collapse(2) -#endif - for (int_t k1 = 0; k1 < END1; k1 += STEP1) - for (int_t k2 = 0; k2 < END2; k2++) { - std::forward(func)(k1, k2, mat); - } - } -} - // Dynamic N-qubit template -template -void QubitVector::apply_matrix_lambda(Lambda&& func, const reg_t &qubits, const cvector_t &mat) { +template +void QubitVector::apply_matrix_lambda(Lambda&& func, const list_t &qubits, const cvector_t &mat) { // Error checking #ifdef DEBUG @@ -1317,35 +1267,35 @@ void QubitVector::apply_permutation_matrix(const reg_t& qubits, template void QubitVector::apply_x(const uint_t qubit) { // Lambda function for optimized Pauli-X gate - auto lambda = [&](const int_t k1, const int_t k2)->void { - const auto i0 = k1 | k2; - const auto i1 = i0 | BITS[qubit]; + auto lambda = [&](const areg_t<2> &inds)->void { + const auto i0 = inds[0]; + const auto i1 = inds[1]; std::swap(data_[i0], data_[i1]); }; - apply_lambda(lambda, qubit); + apply_lambda(lambda, (areg_t<1>){qubit}); } template void QubitVector::apply_y(const uint_t qubit) { // Lambda function for optimized Pauli-Y gate const complex_t I(0., 1.); - auto lambda = [&](const int_t k1, const int_t k2)->void { - const auto i0 = k1 | k2; - const auto i1 = i0 | BITS[qubit]; + auto lambda = [&](const areg_t<2> &inds)->void { + const auto i0 = inds[0]; + const auto i1 = inds[1]; const complex_t cache = data_[i0]; data_[i0] = -I * data_[i1]; // mat(0,1) data_[i1] = I * cache; // mat(1,0) }; - apply_lambda(lambda, qubit); + apply_lambda(lambda, (areg_t<1>){qubit}); } template void QubitVector::apply_z(const uint_t qubit) { // Lambda function for optimized Pauli-Z gate - auto lambda = [&](const int_t k1, const int_t k2)->void { - data_[k1 | k2 | BITS[qubit]] *= complex_t(-1.0, 0.0); + auto lambda = [&](const areg_t<2> &inds)->void { + data_[inds[1]] *= complex_t(-1.0, 0.0); }; - apply_lambda(lambda, qubit); + apply_lambda(lambda, (areg_t<1>){qubit}); } //------------------------------------------------------------------------------ @@ -1358,6 +1308,16 @@ void QubitVector::apply_mcx(const reg_t &qubits) { const size_t N = qubits.size(); const size_t pos0 = MASKS[N - 1]; const size_t pos1 = MASKS[N]; + + if (N == 2) { + // Lambda function for multi-controlled X gate + auto lambda = [&](const areg_t<4> &inds)->void { + std::swap(data_[inds[pos0]], data_[inds[pos1]]); + }; + apply_lambda(lambda, (areg_t<2>){qubits[0], qubits[1]}); + return; + } + // Lambda function for multi-controlled X gate auto lambda = [&](const indexes_t &inds)->void { std::swap(data_[inds[pos0]], data_[inds[pos1]]); @@ -1372,6 +1332,18 @@ void QubitVector::apply_mcy(const reg_t &qubits) { const size_t pos0 = MASKS[N - 1]; const size_t pos1 = MASKS[N]; const complex_t I(0., 1.); + + if (N == 2) { + // Lambda function for multi-controlled Y gate + auto lambda = [&](const areg_t<4> &inds)->void { + const complex_t cache = data_[inds[pos0]]; + data_[inds[pos0]] = -I * data_[inds[pos1]]; + data_[inds[pos1]] = I * cache; + }; + apply_lambda(lambda, (areg_t<2>){qubits[0], qubits[1]}); + return; + } + // Lambda function for multi-controlled Y gate auto lambda = [&](const indexes_t &inds)->void { const complex_t cache = data_[inds[pos0]]; @@ -1383,6 +1355,18 @@ void QubitVector::apply_mcy(const reg_t &qubits) { template void QubitVector::apply_mcz(const reg_t &qubits) { + const size_t N = qubits.size(); + + if (N == 2) { + // Lambda function for multi-controlled Z gate + auto lambda = [&](const areg_t<4> &inds)->void { + // Multiply last block index by -1 + data_[inds[MASKS[qubits.size()]]] *= -1.; + }; + apply_lambda(lambda, (areg_t<2>){qubits[0], qubits[1]}); + return; + } + // Lambda function for multi-controlled Z gate auto lambda = [&](const indexes_t &inds)->void { // Multiply last block index by -1 @@ -1398,6 +1382,19 @@ void QubitVector::apply_mcu(const reg_t &qubits, const size_t N = qubits.size(); const size_t pos0 = MASKS[N - 1]; const size_t pos1 = MASKS[N]; + + if (N == 2) { + // Lambda function for multi-controlled single-qubit gate + auto lambda = [&](const areg_t<4> &inds, + const cvector_t &_mat)->void { + const auto cache = data_[pos0]; + data_[pos0] = _mat[0] * data_[pos0] + _mat[2] * data_[pos1]; + data_[pos1] = _mat[1] * cache + _mat[3] * data_[pos1]; + }; + apply_matrix_lambda(lambda, (areg_t<2>){qubits[0], qubits[1]}, mat); + return; + } + // Lambda function for multi-controlled single-qubit gate auto lambda = [&](const indexes_t &inds, const cvector_t &_mat)->void { @@ -1440,109 +1437,108 @@ void QubitVector::apply_matrix(const uint_t qubit, } // Lambda function for single-qubit matrix multiplication - const int_t BIT = BITS[qubit]; - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto pos0 = k1 | k2; - const auto pos1 = pos0 | BIT; + const auto pos0 = inds[0]; + const auto pos1 = inds[1]; const auto cache = data_[pos0]; data_[pos0] = _mat[0] * data_[pos0] + _mat[2] * data_[pos1]; data_[pos1] = _mat[1] * cache + _mat[3] * data_[pos1]; }; - apply_matrix_lambda(lambda, qubit, mat); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, mat); } template void QubitVector::apply_diagonal_matrix(const uint_t qubit, const cvector_t& diag) { // TODO: This should be changed so it isn't checking doubles with == - const int_t BIT = BITS[qubit]; if (diag[0] == 1.0) { // [[1, 0], [0, z]] matrix if (diag[1] == 1.0) return; // Identity if (diag[1] == complex_t(0., -1.)) { // [[1, 0], [0, -i]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2 | BIT; + const auto k = inds[1]; double cache = data_[k].imag(); data_[k].imag(data_[k].real() * -1.); data_[k].real(cache); }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else if (diag[1] == complex_t(0., 1.)) { // [[1, 0], [0, i]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2 | BIT; + const auto k = inds[1]; double cache = data_[k].imag(); data_[k].imag(data_[k].real()); data_[k].real(cache * -1.); }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else if (diag[0] == 0.0) { // [[1, 0], [0, 0]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - data_[k1 | k2 | BIT] = 0.0; + data_[inds[1]] = 0.0; }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else { // general [[1, 0], [0, z]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2 | BIT; + const auto k = inds[1]; data_[k] *= _mat[1]; }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } } else if (diag[1] == 1.0) { // [[z, 0], [0, 1]] matrix if (diag[0] == complex_t(0., -1.)) { // [[-i, 0], [0, 1]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2; + const auto k = inds[1]; double cache = data_[k].imag(); data_[k].imag(data_[k].real() * -1.); data_[k].real(cache); }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else if (diag[0] == complex_t(0., 1.)) { // [[i, 0], [0, 1]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2; + const auto k = inds[1]; double cache = data_[k].imag(); data_[k].imag(data_[k].real()); data_[k].real(cache * -1.); }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else if (diag[0] == 0.0) { // [[0, 0], [0, 1]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - data_[k1 | k2] = 0.0; + data_[inds[0]] = 0.0; }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } else { // general [[z, 0], [0, 1]] - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2; + const auto k = inds[0]; data_[k] *= _mat[0]; }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } } else { // Lambda function for diagonal matrix multiplication - auto lambda = [&](const int_t k1, const int_t k2, + auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat)->void { - const auto k = k1 | k2; - data_[k] *= _mat[0]; - data_[k | BIT] *= _mat[1]; + const auto k0 = inds[0]; + const auto k1 = inds[1]; + data_[k0] *= _mat[0]; + data_[k1] *= _mat[1]; }; - apply_matrix_lambda(lambda, qubit, diag); + apply_matrix_lambda(lambda, (areg_t<1>) {qubit}, diag); } }