The Quantum Exact Simulation Toolkit v4.2.0
Loading...
Searching...
No Matches
evolve.cpp
1/** @file
2 * Testing utilities which evolve a reference state
3 * (qvector or qmatrix) under the action of a
4 * reference operation. These are slow, serial,
5 * un-optimised, defensively-designed routines.
6 *
7 * @author Tyson Jones
8 */
9
10#include "qvector.hpp"
11#include "qmatrix.hpp"
12#include "macros.hpp"
13#include "linalg.hpp"
14
15#include <tuple>
16#include <vector>
17#include <algorithm>
18
19using std::vector;
20
21
22
23/*
24 * OPERATOR MATRICES
25 */
26
27
28qmatrix getSwapMatrix(int qb1, int qb2, int numQb) {
29 DEMAND( numQb > 1 );
30 DEMAND( (qb1 >= 0 && qb1 < numQb) );
31 DEMAND( (qb2 >= 0 && qb2 < numQb) );
32
33 if (qb1 == qb2)
34 return getIdentityMatrix(getPow2(numQb));
35
36 if (qb1 > qb2)
37 std::swap(qb1, qb2);
38
39 qmatrix out;
40
41 // qubits are either adjacent
42 if (qb2 == qb1 + 1) {
43 out = qmatrix{{1,0,0,0},{0,0,1,0},{0,1,0,0},{0,0,0,1}};
44
45 // or distant
46 } else {
47 int block = getPow2(qb2 - qb1);
48 out = getZeroMatrix(block*2);
49 qmatrix iden = getIdentityMatrix(block/2);
50
51 // Lemma 3.1 of arxiv.org/pdf/1711.09765.pdf
52 qmatrix p0{{1,0},{0,0}};
53 qmatrix l0{{0,1},{0,0}};
54 qmatrix l1{{0,0},{1,0}};
55 qmatrix p1{{0,0},{0,1}};
56
57 // notating a^(n+1) = identity(getPow2(n)) (otimes) a, we construct the matrix
58 // [ p0^(N) l1^N ]
59 // [ l0^(N) p1^N ]
60 // where N = qb2 - qb1 */
61 setSubMatrix(out, getKroneckerProduct(iden, p0), 0, 0);
62 setSubMatrix(out, getKroneckerProduct(iden, l0), block, 0);
63 setSubMatrix(out, getKroneckerProduct(iden, l1), 0, block);
64 setSubMatrix(out, getKroneckerProduct(iden, p1), block, block);
65 }
66
67 // pad swap with outer identities
68 if (qb1 > 0)
69 out = getKroneckerProduct(out, getIdentityMatrix(getPow2(qb1)));
70
71 if (qb2 < numQb-1)
72 out = getKroneckerProduct(getIdentityMatrix(getPow2(numQb-qb2-1)), out);
73
74 return out;
75}
76
77
78auto getSwapAndUnswapMatrices(vector<int> ctrls, vector<int> targs, size_t numQubits) {
79 DEMAND( numQubits >= ctrls.size() + targs.size() );
80
81 // matrices which swap targs+ctrls to be contiguous from 0
82 qmatrix swaps = getIdentityMatrix(getPow2(numQubits));
83 qmatrix unswaps = getIdentityMatrix(getPow2(numQubits));
84
85 // swap targs to {0, ..., ntargs - 1}
86 for (size_t i=0; i<targs.size(); i++) {
87
88 if (i == (size_t) targs[i])
89 continue;
90
91 qmatrix m = getSwapMatrix(i, targs[i], numQubits);
92 swaps = m * swaps;
93 unswaps = unswaps * m;
94
95 std::replace(ctrls.begin(), ctrls.end(), (int) i, targs[i]);
96 std::replace(targs.begin(), targs.end(), (int) i, targs[i]);
97 }
98
99 // swap ctrls to {ntargs, ..., ntargs + nctrls - 1}
100 for (size_t i=0; i<ctrls.size(); i++) {
101
102 size_t j = i + targs.size();
103 if (j == (size_t) ctrls[i])
104 continue;
105
106 qmatrix m = getSwapMatrix(j, ctrls[i], numQubits);
107 swaps = m * swaps;
108 unswaps = unswaps * m;
109
110 std::replace(ctrls.begin(), ctrls.end(), (int) j, ctrls[i]);
111 }
112
113 return std::tuple{swaps, unswaps};
114}
115
116
117qmatrix getFullStateOperator(vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix, size_t numQubits) {
118 DEMAND( numQubits >= ctrls.size() + targs.size() );
119 DEMAND( getPow2(targs.size()) == (qindex) matrix.size() );
120 DEMAND( (ctrlStates.empty() || ctrlStates.size() == ctrls.size()) );
121
122 // construct controlled-(matrix) upon lowest order qubits
123 qmatrix full = getControlledMatrix(matrix, ctrls.size());
124
125 // left-pad 'full' to be numQubits large
126 if (numQubits > ctrls.size() + targs.size()) {
127 size_t pad = getPow2(numQubits - ctrls.size() - targs.size());
128 full = getKroneckerProduct(getIdentityMatrix(pad), full);
129 }
130
131 // apply swaps to retarget 'full' to given ctrls and targs
132 auto [swaps, unswaps] = getSwapAndUnswapMatrices(ctrls, targs, numQubits);
133 qmatrix out = unswaps * full * swaps;
134
135 // apply NOT to all zero-controlled qubits (recurses just once)
136 qmatrix matrX = {{0,1},{1,0}};
137 for (size_t i=0; i<ctrlStates.size(); i++) {
138
139 if (ctrlStates[i] == 1)
140 continue;
141
142 qmatrix fullX = getFullStateOperator({}, {}, {ctrls[i]}, matrX, numQubits);
143 out = fullX * out * fullX;
144 }
145
146 return out;
147}
148
149
150
151/*
152 * EVOLUTION
153 */
154
155
156// overloads with no targs (given full operator)
157
158void applyReferenceOperator(qvector& state, qmatrix matrix) {
159 DEMAND( state.size() == matrix.size() );
160
161 state = matrix * state;
162}
163void applyReferenceOperator(qmatrix& state, qmatrix matrix) {
164 DEMAND( state.size() == matrix.size() );
165
166 state = matrix * state * getConjugateTranspose(matrix);
167}
168
169void leftapplyReferenceOperator(qvector& state, qmatrix matrix) {
170 DEMAND( state.size() == matrix.size() );
171
172 // for statevectors, multiplying is the same as applying
173 applyReferenceOperator(state, matrix);
174}
175
176void leftapplyReferenceOperator(qmatrix& state, qmatrix matrix) {
177 DEMAND( state.size() == matrix.size() );
178
179 // we left-multiply upon density matrices only
180 state = matrix * state;
181}
182
183void rightapplyReferenceOperator(qmatrix& state, qmatrix matrix) {
184 DEMAND( state.size() == matrix.size() );
185
186 // we right-multiply upon density matrices only
187 state = state * matrix;
188}
189
190
191// overloads with ctrls, states and targs (given sub-operator)
192
193void applyReferenceOperator(qvector& state, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix) {
194
195 qmatrix fullOp = getFullStateOperator(ctrls, ctrlStates, targs, matrix, getLog2(state.size()));
196 applyReferenceOperator(state, fullOp);
197}
198
199void applyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix) {
200
201 qmatrix fullOp = getFullStateOperator(ctrls, ctrlStates, targs, matrix, getLog2(state.size()));
202 applyReferenceOperator(state, fullOp);
203}
204
205void leftapplyReferenceOperator(qvector& state, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix) {
206
207 applyReferenceOperator(state, ctrls, ctrlStates, targs, matrix);
208}
209
210void leftapplyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix) {
211
212 qmatrix left = getFullStateOperator(ctrls, ctrlStates, targs, matrix, getLog2(state.size()));
213 leftapplyReferenceOperator(state, left);
214}
215
216void rightapplyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> ctrlStates, vector<int> targs, qmatrix matrix) {
217
218 qmatrix left = getFullStateOperator(ctrls, ctrlStates, targs, matrix, getLog2(state.size()));
219 rightapplyReferenceOperator(state, left);
220}
221
222
223// overloads with only ctrls and targs
224
225void applyReferenceOperator(qvector& state, vector<int> ctrls, vector<int> targs, qmatrix matrix) {
226
227 applyReferenceOperator(state, ctrls, {}, targs, matrix);
228}
229void applyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> targs, qmatrix matrix) {
230
231 applyReferenceOperator(state, ctrls, {}, targs, matrix);
232}
233void leftapplyReferenceOperator(qvector& state, vector<int> ctrls, vector<int> targs, qmatrix matrix) {
234
235 leftapplyReferenceOperator(state, ctrls, {}, targs, matrix);
236}
237void leftapplyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> targs, qmatrix matrix) {
238
239 leftapplyReferenceOperator(state, ctrls, {}, targs, matrix);
240}
241void rightapplyReferenceOperator(qmatrix& state, vector<int> ctrls, vector<int> targs, qmatrix matrix) {
242
243 rightapplyReferenceOperator(state, ctrls, {}, targs, matrix);
244}
245
246
247// overloads with only targs
248
249void applyReferenceOperator(qvector& state, vector<int> targs, qmatrix matrix) {
250
251 applyReferenceOperator(state, {}, {}, targs, matrix);
252}
253void applyReferenceOperator(qmatrix& state, vector<int> targs, qmatrix matrix) {
254
255 applyReferenceOperator(state, {}, {}, targs, matrix);
256}
257void leftapplyReferenceOperator(qvector& state, vector<int> targs, qmatrix matrix) {
258
259 leftapplyReferenceOperator(state, {}, {}, targs, matrix);
260}
261void leftapplyReferenceOperator(qmatrix& state, vector<int> targs, qmatrix matrix) {
262
263 leftapplyReferenceOperator(state, {}, {}, targs, matrix);
264}
265void rightapplyReferenceOperator(qmatrix& state, vector<int> targs, qmatrix matrix) {
266
267 rightapplyReferenceOperator(state, {}, {}, targs, matrix);
268}
269
270
271// overloads with only targs and kraus operators
272
273void applyReferenceOperator(qmatrix& state, vector<int> targs, vector<qmatrix> matrices) {
274
275 qmatrix in = state;
276 qmatrix out = getZeroMatrix(state.size());
277
278 for (auto& matrix : matrices) {
279 state = in;
280 applyReferenceOperator(state, targs, matrix);
281 out += state;
282 }
283
284 state = out;
285}
qmatrix getSwapMatrix(int qb1, int qb2, int numQb)
Definition evolve.cpp:28
qmatrix getKroneckerProduct(qmatrix a, qmatrix b)
Definition linalg.cpp:523
qmatrix getConjugateTranspose(qmatrix m)
Definition linalg.cpp:291
qmatrix getIdentityMatrix(size_t dim)
Definition qmatrix.cpp:30
void setSubMatrix(qmatrix &dest, qmatrix sub, size_t r, size_t c)
Definition qmatrix.cpp:203
qmatrix getZeroMatrix(size_t dim)
Definition qmatrix.cpp:18