Skip to content

Commit 2d8b9d8

Browse files
vsoftcovsoftco
committed
Update
Signed-off-by: Vlad Gheorghiu <vsoftco@gmail.com>
1 parent 991b22c commit 2d8b9d8

7 files changed

Lines changed: 369 additions & 38 deletions

File tree

‎CHANGES.md‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,9 @@
22

33
- Significant performance improvements (10x to 1000x) for qubit state
44
manipulation and qubit circuit execution through extensive optimizations of
5-
`qpp::apply()`, `qpp::applyCTRL()`, and `qpp::measure_seq()`.
5+
`qpp::apply()`, `qpp::applyCTRL()`, `qpp::applyCTRL_fan()`,
6+
`qpp::measure_seq()`, `qpp::ptrace()`, `qpp::ptranspose()`, and
7+
`qpp::syspermute()`.
68
- **pyqpp** is now feature complete, providing Python users with
79
- Added `pyqpp.QEngine.get_stats_to_JSON()` to display the measurement
810
statistics in JSON format from **pyqpp**

‎benchmarks/src/apply_ctrl_fan_psi_bench.cpp‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
#define CATCH_CONFIG_RUNNER
2+
3+
#include <string>
4+
5+
#include <catch2/catch_all.hpp>
6+
7+
#include <qpp/internal/kernels/qubit/apply_ctrl_fan.hpp>
8+
#include <qpp/qpp.hpp>
9+
10+
#include "openmp_utils.hpp"
11+
12+
namespace {
13+
qpp::idx nq = 10; // default number of qubits if none provided at runtime
14+
} // namespace
15+
16+
int main(int argc, char* argv[]) {
17+
Catch::Session session;
18+
19+
auto cli = session.cli() |
20+
Catch::Clara::Opt(nq, "qubits")["--nq"]("Number of qubits");
21+
#ifdef QPP_OPENMP
22+
cli |= Catch::Clara::Opt(cli_threads, "threads")["--threads"](
23+
"Number of OpenMP threads/cores (ignored if the environment variable "
24+
"OMP_NUM_THREADS is set)");
25+
#endif // QPP_OPENMP
26+
session.cli(cli);
27+
28+
int returnCode = session.applyCommandLine(argc, argv);
29+
if (returnCode != 0) {
30+
return returnCode;
31+
}
32+
33+
if (session.config().showHelp()) {
34+
return 0;
35+
}
36+
37+
#ifdef QPP_OPENMP
38+
set_openmp_threads();
39+
#endif // QPP_OPENMP
40+
41+
return session.run();
42+
}
43+
44+
TEST_CASE("qpp::applyCTRL_fan() state vector benchmark",
45+
"[benchmark][applyCTRL_fan-state-vector]") {
46+
// Setup (NOT measured)
47+
REQUIRE(nq > 1);
48+
qpp::idx D = qpp::internal::safe_pow<qpp::idx>(2, nq);
49+
qpp::ket psi = qpp::randket(D); // D x 1 random state vector
50+
// Test worst-case scenario
51+
qpp::idx m = nq / 2;
52+
std::vector<qpp::idx> ctrl(m);
53+
std::vector<qpp::idx> target(m);
54+
std::vector<qpp::idx> shift(nq, 0);
55+
for (qpp::idx i = 0; i < nq / 2; ++i) {
56+
ctrl[i] = i;
57+
target[i] = m + i;
58+
}
59+
60+
// Benchmarked portion (executed repeatedly)
61+
BENCHMARK("CTRL-fan (psi) nq=" + std::to_string(nq)) {
62+
// CRITICAL: Return the result so the compiler doesn't optimize the
63+
// calculation away.
64+
return qpp::applyCTRL_fan(psi, qpp::gt.H, ctrl, target);
65+
};
66+
67+
// Benchmarked portion (executed repeatedly)
68+
BENCHMARK("CTRL-fan qubit (psi) nq=" + std::to_string(nq)) {
69+
// CRITICAL: Return the result so the compiler doesn't optimize the
70+
// calculation away.
71+
return qpp::internal::kernels::qubit::apply_ctrl_fan_psi(
72+
psi, qpp::gt.H, ctrl, target, shift, nq);
73+
};
74+
}

‎benchmarks/src/apply_ctrl_fan_rho_bench.cpp‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
#define CATCH_CONFIG_RUNNER
2+
3+
#include <string>
4+
5+
#include <catch2/catch_all.hpp>
6+
7+
#include <qpp/internal/kernels/qubit/apply_ctrl_fan.hpp>
8+
#include <qpp/qpp.hpp>
9+
10+
#include "openmp_utils.hpp"
11+
12+
namespace {
13+
qpp::idx nq = 10; // default number of qubits if none provided at runtime
14+
} // namespace
15+
16+
int main(int argc, char* argv[]) {
17+
Catch::Session session;
18+
19+
auto cli = session.cli() |
20+
Catch::Clara::Opt(nq, "qubits")["--nq"]("Number of qubits");
21+
#ifdef QPP_OPENMP
22+
cli |= Catch::Clara::Opt(cli_threads, "threads")["--threads"](
23+
"Number of OpenMP threads/cores (ignored if the environment variable "
24+
"OMP_NUM_THREADS is set)");
25+
#endif // QPP_OPENMP
26+
session.cli(cli);
27+
28+
int returnCode = session.applyCommandLine(argc, argv);
29+
if (returnCode != 0) {
30+
return returnCode;
31+
}
32+
33+
if (session.config().showHelp()) {
34+
return 0;
35+
}
36+
37+
#ifdef QPP_OPENMP
38+
set_openmp_threads();
39+
#endif // QPP_OPENMP
40+
41+
return session.run();
42+
}
43+
44+
TEST_CASE("qpp::applyCTRL_fan() density matrix benchmark",
45+
"[benchmark][applyCTRL_fan-density-matrix]") {
46+
// Setup (NOT measured)
47+
REQUIRE(nq > 1);
48+
qpp::idx D = qpp::internal::safe_pow<qpp::idx>(2, nq);
49+
qpp::cmat rho = qpp::randrho(D); // D x D random density matrix
50+
// Test worst-case scenario
51+
qpp::idx mid = nq / 2;
52+
std::vector<qpp::idx> ctrl(mid);
53+
std::vector<qpp::idx> target(mid);
54+
std::vector<qpp::idx> shift(nq, 0);
55+
for (qpp::idx i = 0; i < mid; ++i) {
56+
ctrl[i] = i;
57+
target[i] = mid + i;
58+
}
59+
60+
// Benchmarked portion (executed repeatedly)
61+
BENCHMARK("CTRL-fan (rho) nq=" + std::to_string(nq)) {
62+
// CRITICAL: Return the result so the compiler doesn't optimize the
63+
// calculation away.
64+
return qpp::applyCTRL_fan(rho, qpp::gt.H, ctrl, target);
65+
};
66+
67+
// Benchmarked portion (executed repeatedly)
68+
BENCHMARK("CTRL-fan qubit (rho) nq=" + std::to_string(nq)) {
69+
// CRITICAL: Return the result so the compiler doesn't optimize the
70+
// calculation away.
71+
return qpp::internal::kernels::qubit::apply_ctrl_fan_rho(
72+
rho, qpp::gt.H, ctrl, target, shift, nq);
73+
};
74+
}

‎include/qpp/classes/gates.hpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -664,7 +664,8 @@ class Gates final : public internal::Singleton<const Gates> // const Singleton
664664
if (elem >= d) {
665665
throw exception::OutOfRange("qpp::Gates::CTRL()", "shift");
666666
}
667-
elem = d - elem;
667+
// invert shift mod D
668+
elem = (d - elem) % d;
668669
}
669670
}
670671
// END EXCEPTION CHECKS

‎include/qpp/internal/kernels/qubit/apply_ctrl_fan.hpp‎

Lines changed: 14 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,7 @@ template <typename Derived1, typename Derived2>
6262
apply_ctrl_fan_psi(const Eigen::MatrixBase<Derived1>& state,
6363
const Eigen::MatrixBase<Derived2>& A,
6464
const std::vector<idx>& ctrl, const std::vector<idx>& target,
65-
idx n,
66-
std::optional<std::vector<idx>> shift = std::nullopt) {
65+
const std::vector<idx>& shift, idx n) {
6766

6867
const expr_t<Derived1>& rstate = state.derived();
6968
const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
@@ -80,14 +79,12 @@ apply_ctrl_fan_psi(const Eigen::MatrixBase<Derived1>& state,
8079
assert(rA.rows() == 2 && "Aggressively optimized for qubits (d=2) only");
8180

8281
idx D = 1ULL << n; // Total dimension 2^n
83-
assert(rstate.rows() == D && "State vector dimension mismatch with 2^n");
82+
assert(static_cast<idx>(rstate.rows()) == D &&
83+
"State vector dimension mismatch with 2^n");
8484

85-
if (shift.has_value()) {
86-
assert(shift.value().size() == ctrl.size() &&
87-
"Shift size must match ctrl size");
88-
for (idx s : shift.value()) {
89-
assert(s < 2 && "Shift values must be < 2 for qubits");
90-
}
85+
assert(shift.size() == ctrl.size() && "Shift size must match ctrl size");
86+
for (idx s : shift) {
87+
assert(s < 2 && "Shift values must be < 2 for qubits");
9188
}
9289

9390
if (D == 1) {
@@ -102,9 +99,7 @@ apply_ctrl_fan_psi(const Eigen::MatrixBase<Derived1>& state,
10299
for (size_t i = 0; i < ctrl.size(); ++i) {
103100
idx bit = 1ULL << (n - 1 - ctrl[i]); // MSB convention index
104101
ctrl_mask |= bit;
105-
106-
idx s = shift.has_value() ? shift.value()[i] : 0;
107-
if (s == 0) {
102+
if (shift[i] == 0) {
108103
ctrl_val |= bit; // Default to triggering on |1>
109104
}
110105
}
@@ -169,8 +164,7 @@ template <typename Derived1, typename Derived2>
169164
apply_ctrl_fan_rho(const Eigen::MatrixBase<Derived1>& state,
170165
const Eigen::MatrixBase<Derived2>& A,
171166
const std::vector<idx>& ctrl, const std::vector<idx>& target,
172-
idx n,
173-
std::optional<std::vector<idx>> shift = std::nullopt) {
167+
const std::vector<idx>& shift, idx n) {
174168

175169
const expr_t<Derived1>& rstate = state.derived();
176170
const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
@@ -187,15 +181,13 @@ apply_ctrl_fan_rho(const Eigen::MatrixBase<Derived1>& state,
187181
assert(rA.rows() == 2 && "Aggressively optimized for qubits (d=2) only");
188182

189183
idx D = 1ULL << n; // Total dimension 2^n
190-
assert(rstate.rows() == D && rstate.cols() == D &&
184+
assert(static_cast<idx>(rstate.rows()) == D &&
185+
static_cast<idx>(rstate.cols()) == D &&
191186
"Density matrix dimension mismatch with 2^n");
192187

193-
if (shift.has_value()) {
194-
assert(shift.value().size() == ctrl.size() &&
195-
"Shift size must match ctrl size");
196-
for (idx s : shift.value()) {
197-
assert(s < 2 && "Shift values must be < 2 for qubits");
198-
}
188+
assert(shift.size() == ctrl.size() && "Shift size must match ctrl size");
189+
for (idx s : shift) {
190+
assert(s < 2 && "Shift values must be < 2 for qubits");
199191
}
200192

201193
if (D == 1) {
@@ -210,9 +202,7 @@ apply_ctrl_fan_rho(const Eigen::MatrixBase<Derived1>& state,
210202
for (size_t i = 0; i < ctrl.size(); ++i) {
211203
idx bit = 1ULL << (n - 1 - ctrl[i]); // MSB convention index
212204
ctrl_mask |= bit;
213-
214-
idx s = shift.has_value() ? shift.value()[i] : 0;
215-
if (s == 0) {
205+
if (shift[i] == 0) {
216206
ctrl_val |= bit; // Default to triggering on |1>
217207
}
218208
}

‎include/qpp/operations.hpp‎

Lines changed: 41 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@
5252
#include "qpp/classes/gates.hpp"
5353
#include "qpp/internal/kernels/qubit/apply.hpp"
5454
#include "qpp/internal/kernels/qubit/apply_ctrl.hpp"
55+
#include "qpp/internal/kernels/qubit/apply_ctrl_fan.hpp"
5556
#include "qpp/internal/kernels/qubit/ptrace.hpp"
5657
#include "qpp/internal/kernels/qubit/ptranspose.hpp"
5758
#include "qpp/internal/kernels/qubit/syspermute.hpp"
@@ -730,7 +731,8 @@ applyCTRL(const Eigen::MatrixBase<Derived1>& state,
730731
if (elem >= d) {
731732
throw exception::OutOfRange("qpp::applyCTRL()", "shift");
732733
}
733-
elem = d - elem;
734+
// invert shift mod D
735+
elem = (d - elem) % d;
734736
}
735737
}
736738
// END EXCEPTION CHECKS
@@ -806,7 +808,7 @@ applyCTRL(const Eigen::MatrixBase<Derived1>& state,
806808
D);
807809
// project on complement
808810
dyn_col_vect<typename Derived1::Scalar> chopped_psi_bar =
809-
state / d - chopped_psi;
811+
(state / d) - chopped_psi;
810812
// apply gates on target if non-zero
811813
if (chopped_psi !=
812814
dyn_col_vect<typename Derived1::Scalar>::Zero(D)) {
@@ -1084,14 +1086,42 @@ applyCTRL_fan(const Eigen::MatrixBase<Derived1>& state,
10841086
if (elem >= d) {
10851087
throw exception::OutOfRange("qpp::applyCTRL_fan()", "shift");
10861088
}
1087-
elem = d - elem;
1089+
// invert shift mod D
1090+
elem = (d - elem) % d;
10881091
}
10891092
}
10901093
// END EXCEPTION CHECKS
10911094
if (!shift.has_value()) {
10921095
shift = std::vector<idx>(ctrl.size(), 0);
10931096
}
10941097

1098+
// qubit optimizations
1099+
#ifdef QPP_QUBIT_OPTIMIZATIONS
1100+
if (internal::all_qubits(dims)) {
1101+
idx n = dims.size();
1102+
// ket
1103+
if (internal::check_cvector(rstate)) {
1104+
// check that dims match state vector
1105+
if (!internal::check_dims_match_cvect(dims, rstate)) {
1106+
throw exception::DimsMismatchCvector("qpp::applyCTRL_fan()",
1107+
"dims/state");
1108+
}
1109+
return internal::kernels::qubit::apply_ctrl_fan_psi(
1110+
state, A, ctrl, target, shift.value(), n);
1111+
1112+
}
1113+
// density matrix
1114+
else {
1115+
// check that dims match density matrix
1116+
if (!internal::check_dims_match_mat(dims, rstate)) {
1117+
throw exception::DimsMismatchCvector("qpp::applyCTRL_fan()",
1118+
"dims/state");
1119+
}
1120+
return internal::kernels::qubit::apply_ctrl_fan_rho(
1121+
state, A, ctrl, target, shift.value(), n);
1122+
}
1123+
}
1124+
#endif // QPP_QUBIT_OPTIMIZATIONS
10951125
// construct the table of A^k
10961126
std::vector<dyn_mat<typename Derived1::Scalar>> Ak;
10971127
for (idx k = 0; k < d; ++k) {
@@ -1124,7 +1154,7 @@ applyCTRL_fan(const Eigen::MatrixBase<Derived1>& state,
11241154
D);
11251155
// project on complement
11261156
dyn_col_vect<typename Derived1::Scalar> chopped_psi_bar =
1127-
state / d - chopped_psi;
1157+
(state / d) - chopped_psi;
11281158
// apply gates on target if non-zero
11291159
if (chopped_psi !=
11301160
dyn_col_vect<typename Derived1::Scalar>::Zero(D)) {
@@ -1198,7 +1228,7 @@ applyCTRL_fan(const Eigen::MatrixBase<Derived1>& state,
11981228
#pragma omp critical
11991229
#endif // QPP_OPENMP
12001230
{
1201-
result += kron(psi_i, phi_i_bra);
1231+
result += psi_i * phi_i_bra;
12021232
}
12031233
} // end for(i)
12041234

@@ -1300,7 +1330,7 @@ applyCTRL_fan(const Eigen::MatrixBase<Derived1>& state,
13001330
#pragma omp parallel for
13011331
#endif // QPP_OPENMP
13021332
for (idx a = 0; a < Din; ++a) {
1303-
MES(a * Din + a, 0) = 1;
1333+
MES((a * Din) + a, 0) = 1;
13041334
}
13051335

13061336
cmat Omega = MES * adjoint(MES);
@@ -1436,8 +1466,8 @@ inline std::vector<cmat> choi2kraus(const cmat& A) {
14361466
for (idx b = 0; b < Dout; ++b) {
14371467
for (idx m = 0; m < Din; ++m) {
14381468
for (idx n = 0; n < Din; ++n) {
1439-
result(a * Dout + b, m * Din + n) =
1440-
A(m * Dout + a, n * Dout + b);
1469+
result((a * Dout) + b, (m * Din) + n) =
1470+
A((m * Dout) + a, (n * Dout) + b);
14411471
}
14421472
}
14431473
}
@@ -1510,8 +1540,8 @@ inline cmat choi2super(const cmat& A) {
15101540
for (idx b = 0; b < Dout; ++b) {
15111541
for (idx m = 0; m < Din; ++m) {
15121542
for (idx n = 0; n < Din; ++n) {
1513-
result(m * Dout + a, n * Dout + b) =
1514-
A(a * Dout + b, m * Din + n);
1543+
result((m * Dout) + a, (n * Dout) + b) =
1544+
A((a * Dout) + b, (m * Din) + n);
15151545
}
15161546
}
15171547
}
@@ -1581,7 +1611,7 @@ inline cmat choi2super(const cmat& A) {
15811611
for (idx b = 0; b < Dout; ++b) {
15821612
// compute result(ab,mn)=<a|E(|m><n)|b>
15831613
B(b) = 1;
1584-
result(a * Dout + b, m * Din + n) =
1614+
result((a * Dout) + b, (m * Din) + n) =
15851615
static_cast<cmat>(A * EMN * B).value();
15861616
B(b) = 0;
15871617
}

0 commit comments

Comments
 (0)