From 5d9b223201e540d23e7654e4a9ce1bd19b3ff212 Mon Sep 17 00:00:00 2001 From: "Christopher J. Wood" Date: Fri, 26 Apr 2019 13:36:16 -0400 Subject: [PATCH] Add static indexing to qubitvector norm (#178) * add check for diagonal to 1-qubit norm --- src/simulators/statevector/qubitvector.hpp | 156 +++++++++++++++++---- 1 file changed, 130 insertions(+), 26 deletions(-) diff --git a/src/simulators/statevector/qubitvector.hpp b/src/simulators/statevector/qubitvector.hpp index 47ef010d5f..0de4f7c863 100755 --- a/src/simulators/statevector/qubitvector.hpp +++ b/src/simulators/statevector/qubitvector.hpp @@ -1920,51 +1920,148 @@ template double QubitVector::norm(const reg_t &qubits, const cvector_t &mat) const { const uint_t N = qubits.size(); - const uint_t DIM = BITS[N]; + // Error checking #ifdef DEBUG check_vector(mat, 2 * N); #endif - // Lambda function for N-qubit matrix norm - auto lambda = [&](const indexes_t &inds, const cvector_t &_mat, - double &val_re, double &val_im)->void { - (void)val_im; // unused - for (size_t i = 0; i < DIM; i++) { - complex_t vi = 0; - for (size_t j = 0; j < DIM; j++) - vi += _mat[i + DIM * j] * data_[inds[j]]; - val_re += std::real(vi * std::conj(vi)); + // Static array optimized lambda functions + switch (N) { + case 1: + return norm(qubits[0], mat); + case 2: { + // Lambda function for 2-qubit matrix norm + auto lambda = [&](const areg_t<4> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 4; i++) { + complex_t vi = 0; + for (size_t j = 0; j < 4; j++) + vi += _mat[i + 4 * j] * data_[inds[j]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<2> qubits_arr = {{qubits[0], qubits[1]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); } - }; - // Use the lambda function - return std::real(apply_reduction_lambda(lambda, qubits, mat)); + case 3: { + // Lambda function for 3-qubit matrix norm + auto lambda = [&](const areg_t<8> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 8; i++) { + complex_t vi = 0; + for (size_t j = 0; j < 8; j++) + vi += _mat[i + 8 * j] * data_[inds[j]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<3> qubits_arr = {{qubits[0], qubits[1], qubits[2]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); + } + case 4: { + // Lambda function for 4-qubit matrix norm + auto lambda = [&](const areg_t<16> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 16; i++) { + complex_t vi = 0; + for (size_t j = 0; j < 16; j++) + vi += _mat[i + 16 * j] * data_[inds[j]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<4> qubits_arr = {{qubits[0], qubits[1], qubits[2], qubits[3]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); + } + default: { + // Lambda function for N-qubit matrix norm + const uint_t DIM = BITS[N]; + auto lambda = [&](const indexes_t &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < DIM; i++) { + complex_t vi = 0; + for (size_t j = 0; j < DIM; j++) + vi += _mat[i + DIM * j] * data_[inds[j]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + // Use the lambda function + return std::real(apply_reduction_lambda(lambda, qubits, mat)); + } + } // end switch } template double QubitVector::norm_diagonal(const reg_t &qubits, const cvector_t &mat) const { const uint_t N = qubits.size(); - const uint_t DIM = BITS[N]; // Error checking #ifdef DEBUG check_vector(mat, N); #endif - // Lambda function for N-qubit matrix norm - auto lambda = [&](const indexes_t &inds, - const cvector_t &_mat, - double &val_re, - double &val_im)->void { - (void)val_im; // unused - for (size_t i = 0; i < DIM; i++) { - const auto vi = _mat[i] * data_[inds[i]]; - val_re += std::real(vi * std::conj(vi)); + // Static array optimized lambda functions + switch (N) { + case 1: + return norm_diagonal(qubits[0], mat); + case 2: { + // Lambda function for 2-qubit matrix norm + auto lambda = [&](const areg_t<4> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 4; i++) { + const auto vi = _mat[i] * data_[inds[i]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<2> qubits_arr = {{qubits[0], qubits[1]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); } - }; - // Use the lambda function - return std::real(apply_reduction_lambda(lambda, qubits, mat)); + case 3: { + // Lambda function for 3-qubit matrix norm + auto lambda = [&](const areg_t<8> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 8; i++) { + const auto vi = _mat[i] * data_[inds[i]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<3> qubits_arr = {{qubits[0], qubits[1], qubits[2]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); + } + case 4: { + // Lambda function for 4-qubit matrix norm + auto lambda = [&](const areg_t<16> &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < 16; i++) { + const auto vi = _mat[i] * data_[inds[i]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + areg_t<4> qubits_arr = {{qubits[0], qubits[1], qubits[2], qubits[3]}}; + return std::real(apply_reduction_lambda(lambda, qubits_arr, mat)); + } + default: { + // Lambda function for N-qubit matrix norm + const uint_t DIM = BITS[N]; + auto lambda = [&](const indexes_t &inds, const cvector_t &_mat, + double &val_re, double &val_im)->void { + (void)val_im; // unused + for (size_t i = 0; i < DIM; i++) { + const auto vi = _mat[i] * data_[inds[i]]; + val_re += std::real(vi * std::conj(vi)); + } + }; + // Use the lambda function + return std::real(apply_reduction_lambda(lambda, qubits, mat)); + } + } // end switch } //------------------------------------------------------------------------------ @@ -1976,6 +2073,13 @@ double QubitVector::norm(const uint_t qubit, const cvector_t &mat) const #ifdef DEBUG check_vector(mat, 2); #endif + + // Check if input matrix is diagonal, and if so use diagonal function. + if (mat[1] == 0.0 && mat[2] == 0.0) { + const cvector_t diag = {{mat[0], mat[3]}}; + return norm_diagonal(qubit, diag); + } + // Lambda function for norm reduction to real value. auto lambda = [&](const areg_t<2> &inds, const cvector_t &_mat,