XACC
operations.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2020 Vlad Gheorghiu.
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef OPERATIONS_H_
33 #define OPERATIONS_H_
34 
35 namespace qpp {
36 
55 template <typename Derived1, typename Derived2>
56 dyn_mat<typename Derived1::Scalar>
57 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
58  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
59  const std::vector<idx>& target, const std::vector<idx>& dims,
60  std::vector<idx> shift = {}) {
61  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
62  state.derived();
63  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
64 
65  // EXCEPTION CHECKS
66 
67  // check types
68  if (!std::is_same<typename Derived1::Scalar,
69  typename Derived2::Scalar>::value)
70  throw exception::TypeMismatch("qpp::applyCTRL()");
71 
72  // check zero sizes
73  if (!internal::check_nonzero_size(rA))
74  throw exception::ZeroSize("qpp::applyCTRL()");
75 
76  // check zero sizes
77  if (!internal::check_nonzero_size(rstate))
78  throw exception::ZeroSize("qpp::applyCTRL()");
79 
80  // check zero sizes
81  if (!internal::check_nonzero_size(target))
82  throw exception::ZeroSize("qpp::applyCTRL()");
83 
84  // check square matrix for the gate
85  if (!internal::check_square_mat(rA))
86  throw exception::MatrixNotSquare("qpp::applyCTRL()");
87 
88  // check valid state and matching dimensions
89  if (internal::check_cvector(rstate)) {
90  if (!internal::check_dims_match_cvect(dims, state))
91  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
92  } else if (internal::check_square_mat(rstate)) {
93  if (!internal::check_dims_match_mat(dims, state))
94  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
95  } else
96  throw exception::MatrixNotSquareNorCvector("qpp::applyCTRL()");
97 
98  // check that ctrl subsystem is valid w.r.t. dims
99  if (!internal::check_subsys_match_dims(ctrl, dims))
100  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
101 
102  // check that all control subsystems have the same dimension
103  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
104  for (idx i = 1; i < ctrl.size(); ++i)
105  if (dims[ctrl[i]] != d)
106  throw exception::DimsNotEqual("qpp::applyCTRL()");
107 
108  // check that dimension is valid
109  if (!internal::check_dims(dims))
110  throw exception::DimsInvalid("qpp::applyCTRL()");
111 
112  // check that target is valid w.r.t. dims
113  if (!internal::check_subsys_match_dims(target, dims))
114  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
115 
116  // check that gate matches the dimensions of the target
117  std::vector<idx> target_dims(target.size());
118  for (idx i = 0; i < target.size(); ++i)
119  target_dims[i] = dims[target[i]];
120  if (!internal::check_dims_match_mat(target_dims, rA))
121  throw exception::MatrixMismatchSubsys("qpp::applyCTRL()");
122 
123  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
124  ctrlgate.insert(std::end(ctrlgate), std::begin(target), std::end(target));
125  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
126 
127  // check that ctrl + gate subsystem is valid
128  // with respect to local dimensions
129  if (!internal::check_subsys_match_dims(ctrlgate, dims))
130  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
131 
132  // check shift
133  if (!shift.empty() && (shift.size() != ctrl.size()))
134  throw exception::SizeMismatch("qpp::applyCTRL()");
135  if (!shift.empty())
136  for (auto&& elem : shift)
137  if (elem >= d)
138  throw exception::OutOfRange("qpp::applyCTRL()");
139  // END EXCEPTION CHECKS
140 
141  if (shift.empty())
142  shift = std::vector<idx>(ctrl.size(), 0);
143 
144  // construct the table of A^i and (A^dagger)^i
145  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
146  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
147  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i) {
148  Ai.emplace_back(powm(rA, i));
149  Aidagger.emplace_back(powm(adjoint(rA), i));
150  }
151 
152  idx D = static_cast<idx>(rstate.rows()); // total dimension
153  idx n = dims.size(); // total number of subsystems
154  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
155  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
156  idx targetsize = target.size(); // number of subsystems of the target
157  // dimension of ctrl subsystem
158  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
159  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
160 
161  idx Cdims[maxn]; // local dimensions
162  idx CdimsA[maxn]; // local dimensions
163  idx CdimsCTRL[maxn]; // local dimensions
164  idx CdimsCTRLA_bar[maxn]; // local dimensions
165 
166  // compute the complementary subsystem of ctrlgate w.r.t. dims
167  std::vector<idx> ctrlgate_bar = complement(ctrlgate, n);
168  // number of subsystems that are complementary to the ctrl+gate
169  idx ctrlgate_barsize = ctrlgate_bar.size();
170 
171  idx DCTRLA_bar = 1; // dimension of the rest
172  for (idx i = 0; i < ctrlgate_barsize; ++i)
173  DCTRLA_bar *= dims[ctrlgate_bar[i]];
174 
175  for (idx k = 0; k < n; ++k)
176  Cdims[k] = dims[k];
177  for (idx k = 0; k < targetsize; ++k)
178  CdimsA[k] = dims[target[k]];
179  for (idx k = 0; k < ctrlsize; ++k)
180  CdimsCTRL[k] = d;
181  for (idx k = 0; k < ctrlgate_barsize; ++k)
182  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
183 
184  // worker, computes the coefficient and the index for the ket case
185  // used in #pragma omp parallel for collapse
186  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
187  ->std::pair<typename Derived1::Scalar, idx> {
188  idx indx = 0;
189  typename Derived1::Scalar coeff = 0;
190 
191  idx Cmidx[maxn]; // the total multi-index
192  idx CmidxA[maxn]; // the gate part multi-index
193  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
194 
195  // compute the index
196 
197  // set the CTRL part
198  for (idx k = 0; k < ctrlsize; ++k) {
199  Cmidx[ctrl[k]] = (i_ + d - shift[k]) % d;
200  }
201 
202  // set the rest
203  internal::n2multiidx(r_, n - ctrlgatesize, CdimsCTRLA_bar,
204  CmidxCTRLA_bar);
205  for (idx k = 0; k < n - ctrlgatesize; ++k) {
206  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
207  }
208 
209  // set the A part
210  internal::n2multiidx(m_, targetsize, CdimsA, CmidxA);
211  for (idx k = 0; k < targetsize; ++k) {
212  Cmidx[target[k]] = CmidxA[k];
213  }
214 
215  // we now got the total index
216  indx = internal::multiidx2n(Cmidx, n, Cdims);
217 
218  // compute the coefficient
219  for (idx n_ = 0; n_ < DA; ++n_) {
220  internal::n2multiidx(n_, targetsize, CdimsA, CmidxA);
221  for (idx k = 0; k < targetsize; ++k) {
222  Cmidx[target[k]] = CmidxA[k];
223  }
224  coeff +=
225  Ai[i_](m_, n_) * rstate(internal::multiidx2n(Cmidx, n, Cdims));
226  }
227 
228  return std::make_pair(coeff, indx);
229  }; /* end coeff_idx_ket */
230 
231  // worker, computes the coefficient and the index
232  // for the density matrix case
233  // used in #pragma omp parallel for collapse
234  auto coeff_idx_rho = [&](idx i1_, idx m1_, idx r1_, idx i2_, idx m2_,
235  idx r2_) noexcept
236  ->std::tuple<typename Derived1::Scalar, idx, idx> {
237  idx idxrow = 0;
238  idx idxcol = 0;
239  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
240 
241  idx Cmidxrow[maxn]; // the total row multi-index
242  idx Cmidxcol[maxn]; // the total col multi-index
243  idx CmidxArow[maxn]; // the gate part row multi-index
244  idx CmidxAcol[maxn]; // the gate part col multi-index
245  idx CmidxCTRLrow[maxn]; // the control row multi-index
246  idx CmidxCTRLcol[maxn]; // the control col multi-index
247  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
248  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
249 
250  // compute the ket/bra indexes
251 
252  // set the CTRL part
253  internal::n2multiidx(i1_, ctrlsize, CdimsCTRL, CmidxCTRLrow);
254  internal::n2multiidx(i2_, ctrlsize, CdimsCTRL, CmidxCTRLcol);
255 
256  for (idx k = 0; k < ctrlsize; ++k) {
257  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
258  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
259  }
260 
261  // set the rest
262  internal::n2multiidx(r1_, n - ctrlgatesize, CdimsCTRLA_bar,
263  CmidxCTRLA_barrow);
264  internal::n2multiidx(r2_, n - ctrlgatesize, CdimsCTRLA_bar,
265  CmidxCTRLA_barcol);
266  for (idx k = 0; k < n - ctrlgatesize; ++k) {
267  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
268  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
269  }
270 
271  // set the A part
272  internal::n2multiidx(m1_, targetsize, CdimsA, CmidxArow);
273  internal::n2multiidx(m2_, targetsize, CdimsA, CmidxAcol);
274  for (idx k = 0; k < target.size(); ++k) {
275  Cmidxrow[target[k]] = CmidxArow[k];
276  Cmidxcol[target[k]] = CmidxAcol[k];
277  }
278 
279  // we now got the total row/col indexes
280  idxrow = internal::multiidx2n(Cmidxrow, n, Cdims);
281  idxcol = internal::multiidx2n(Cmidxcol, n, Cdims);
282 
283  // check whether all CTRL row and col multi indexes are equal
284  bool all_ctrl_rows_equal = true;
285  bool all_ctrl_cols_equal = true;
286 
287  idx first_ctrl_row, first_ctrl_col;
288  if (ctrlsize > 0) {
289  first_ctrl_row = (CmidxCTRLrow[0] + shift[0]) % d;
290  first_ctrl_col = (CmidxCTRLcol[0] + shift[0]) % d;
291  } else {
292  first_ctrl_row = first_ctrl_col = 1;
293  }
294 
295  for (idx k = 1; k < ctrlsize; ++k) {
296  if ((CmidxCTRLrow[k] + shift[k]) % d != first_ctrl_row) {
297  all_ctrl_rows_equal = false;
298  break;
299  }
300  }
301  for (idx k = 1; k < ctrlsize; ++k) {
302  if ((CmidxCTRLcol[k] + shift[k]) % d != first_ctrl_col) {
303  all_ctrl_cols_equal = false;
304  break;
305  }
306  }
307 
308  // at least one control activated, compute the coefficient
309  for (idx n1_ = 0; n1_ < DA; ++n1_) {
310  internal::n2multiidx(n1_, targetsize, CdimsA, CmidxArow);
311  for (idx k = 0; k < targetsize; ++k) {
312  Cmidxrow[target[k]] = CmidxArow[k];
313  }
314  idx idxrowtmp = internal::multiidx2n(Cmidxrow, n, Cdims);
315 
316  if (all_ctrl_rows_equal) {
317  lhs = Ai[first_ctrl_row](m1_, n1_);
318  } else {
319  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
320  }
321 
322  for (idx n2_ = 0; n2_ < DA; ++n2_) {
323  internal::n2multiidx(n2_, targetsize, CdimsA, CmidxAcol);
324  for (idx k = 0; k < targetsize; ++k) {
325  Cmidxcol[target[k]] = CmidxAcol[k];
326  }
327 
328  if (all_ctrl_cols_equal) {
329  rhs = Aidagger[first_ctrl_col](n2_, m2_);
330  } else {
331  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
332  }
333 
334  idx idxcoltmp = internal::multiidx2n(Cmidxcol, n, Cdims);
335 
336  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
337  }
338  }
339 
340  return std::make_tuple(coeff, idxrow, idxcol);
341  }; /* end coeff_idx_rho */
342 
343  //************ ket ************//
344  if (internal::check_cvector(rstate)) // we have a ket
345  {
346  // check that dims match state vector
347  if (!internal::check_dims_match_cvect(dims, rstate))
348  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
349  if (D == 1)
350  return rstate;
351 
352  dyn_mat<typename Derived1::Scalar> result = rstate;
353 
354 #ifdef WITH_OPENMP_
355 #pragma omp parallel for collapse(2)
356 #endif // WITH_OPENMP_
357  for (idx m = 0; m < DA; ++m)
358  for (idx r = 0; r < DCTRLA_bar; ++r) {
359  if (ctrlsize == 0) // no control
360  {
361  result(coeff_idx_ket(1, m, r).second) =
362  coeff_idx_ket(1, m, r).first;
363  } else
364  for (idx i = 0; i < d; ++i) {
365  result(coeff_idx_ket(i, m, r).second) =
366  coeff_idx_ket(i, m, r).first;
367  }
368  }
369 
370  return result;
371  }
372  //************ density matrix ************//
373  else // we have a density operator
374  {
375  // check that dims match state matrix
376  if (!internal::check_dims_match_mat(dims, rstate))
377  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
378 
379  if (D == 1)
380  return rstate;
381 
382  dyn_mat<typename Derived1::Scalar> result = rstate;
383 
384 #ifdef WITH_OPENMP_
385 #pragma omp parallel for collapse(4)
386 #endif // WITH_OPENMP_
387  for (idx m1 = 0; m1 < DA; ++m1)
388  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
389  for (idx m2 = 0; m2 < DA; ++m2)
390  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
391  if (ctrlsize == 0) // no control
392  {
393  auto coeff_idxes =
394  coeff_idx_rho(1, m1, r1, 1, m2, r2);
395  result(std::get<1>(coeff_idxes),
396  std::get<2>(coeff_idxes)) =
397  std::get<0>(coeff_idxes);
398  } else {
399  for (idx i1 = 0; i1 < Dctrl; ++i1)
400  for (idx i2 = 0; i2 < Dctrl; ++i2) {
401  auto coeff_idxes =
402  coeff_idx_rho(i1, m1, r1, i2, m2, r2);
403  result(std::get<1>(coeff_idxes),
404  std::get<2>(coeff_idxes)) =
405  std::get<0>(coeff_idxes);
406  }
407  }
408 
409  return result;
410  }
411 }
412 
431 template <typename Derived1, typename Derived2>
432 dyn_mat<typename Derived1::Scalar>
433 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
434  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
435  const std::vector<idx>& target, idx d = 2,
436  const std::vector<idx>& shift = {}) {
437  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
438  state.derived();
439  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
440 
441  // EXCEPTION CHECKS
442 
443  // check zero size
444  if (!internal::check_nonzero_size(rstate))
445  throw exception::ZeroSize("qpp::applyCTRL()");
446 
447  // check valid dims
448  if (d < 2)
449  throw exception::DimsInvalid("qpp::applyCTRL()");
450  // END EXCEPTION CHECKS
451 
452  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
453  std::vector<idx> dims(n, d); // local dimensions vector
454 
455  return applyCTRL(rstate, rA, ctrl, target, dims, shift);
456 }
457 
470 template <typename Derived1, typename Derived2>
471 dyn_mat<typename Derived1::Scalar>
472 apply(const Eigen::MatrixBase<Derived1>& state,
473  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& target,
474  const std::vector<idx>& dims) {
475  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
476  state.derived();
477  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
478 
479  // EXCEPTION CHECKS
480 
481  // check types
482  if (!std::is_same<typename Derived1::Scalar,
483  typename Derived2::Scalar>::value)
484  throw exception::TypeMismatch("qpp::apply()");
485 
486  // check zero sizes
487  if (!internal::check_nonzero_size(rA))
488  throw exception::ZeroSize("qpp::apply()");
489 
490  // check zero sizes
491  if (!internal::check_nonzero_size(rstate))
492  throw exception::ZeroSize("qpp::apply()");
493 
494  // check zero sizes
495  if (!internal::check_nonzero_size(target))
496  throw exception::ZeroSize("qpp::apply()");
497 
498  // check square matrix for the gate
499  if (!internal::check_square_mat(rA))
500  throw exception::MatrixNotSquare("qpp::apply()");
501 
502  // check that dimension is valid
503  if (!internal::check_dims(dims))
504  throw exception::DimsInvalid("qpp::apply()");
505 
506  // check that target is valid w.r.t. dims
507  if (!internal::check_subsys_match_dims(target, dims))
508  throw exception::SubsysMismatchDims("qpp::apply()");
509 
510  // check valid state and matching dimensions
511  if (internal::check_cvector(rstate)) {
512  if (!internal::check_dims_match_cvect(dims, state))
513  throw exception::DimsMismatchCvector("qpp::apply()");
514  } else if (internal::check_square_mat(rstate)) {
515  if (!internal::check_dims_match_mat(dims, state))
516  throw exception::DimsMismatchMatrix("qpp::apply()");
517  } else
518  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
519 
520  // check that gate matches the dimensions of the target
521  std::vector<idx> subsys_dims(target.size());
522  for (idx i = 0; i < target.size(); ++i)
523  subsys_dims[i] = dims[target[i]];
524  if (!internal::check_dims_match_mat(subsys_dims, rA))
525  throw exception::MatrixMismatchSubsys("qpp::apply()");
526  // END EXCEPTION CHECKS
527 
528  //************ ket ************//
529  if (internal::check_cvector(rstate)) // we have a ket
530  return applyCTRL(rstate, rA, {}, target, dims);
531  //************ density matrix ************//
532  else // we have a density operator
533  return applyCTRL(rstate, rA, {}, target, dims);
534 }
535 
548 template <typename Derived1, typename Derived2>
549 dyn_mat<typename Derived1::Scalar>
550 apply(const Eigen::MatrixBase<Derived1>& state,
551  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& target,
552  idx d = 2) {
553  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
554  state.derived();
555  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
556 
557  // EXCEPTION CHECKS
558 
559  // check zero size
560  if (!internal::check_nonzero_size(rstate))
561  throw exception::ZeroSize("qpp::apply()");
562 
563  // check valid dims
564  if (d < 2)
565  throw exception::DimsInvalid("qpp::apply()");
566  // END EXCEPTION CHECKS
567 
568  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
569  std::vector<idx> dims(n, d); // local dimensions vector
570 
571  return apply(rstate, rA, target, dims);
572 }
573 
582 template <typename Derived>
583 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
584  const cmat& rA = A.derived();
585 
586  // EXCEPTION CHECKS
587 
588  if (!internal::check_nonzero_size(rA))
589  throw exception::ZeroSize("qpp::apply()");
590  if (!internal::check_square_mat(rA))
591  throw exception::MatrixNotSquare("qpp::apply()");
592  if (Ks.empty())
593  throw exception::ZeroSize("qpp::apply()");
594  if (!internal::check_square_mat(Ks[0]))
595  throw exception::MatrixNotSquare("qpp::apply()");
596  if (Ks[0].rows() != rA.rows())
597  throw exception::DimsMismatchMatrix("qpp::apply()");
598  for (auto&& elem : Ks)
599  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
600  throw exception::DimsNotEqual("qpp::apply()");
601  // END EXCEPTION CHECKS
602 
603  cmat result = cmat::Zero(rA.rows(), rA.rows());
604 
605 #ifdef WITH_OPENMP_
606 #pragma omp parallel for
607 #endif // WITH_OPENMP_
608  for (idx i = 0; i < Ks.size(); ++i) {
609 #ifdef WITH_OPENMP_
610 #pragma omp critical
611 #endif // WITH_OPENMP_
612  { result += Ks[i] * rA * adjoint(Ks[i]); }
613  }
614 
615  return result;
616 }
617 
628 template <typename Derived>
629 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
630  const std::vector<idx>& target, const std::vector<idx>& dims) {
631  const cmat& rA = A.derived();
632 
633  // EXCEPTION CHECKS
634 
635  // check zero sizes
636  if (!internal::check_nonzero_size(rA))
637  throw exception::ZeroSize("qpp::apply()");
638 
639  // check zero sizes
640  if (!internal::check_nonzero_size(target))
641  throw exception::ZeroSize("qpp::apply()");
642 
643  // check square matrix for the A
644  if (!internal::check_square_mat(rA))
645  throw exception::MatrixNotSquare("qpp::apply()");
646 
647  // check that dimension is valid
648  if (!internal::check_dims(dims))
649  throw exception::DimsInvalid("qpp::apply()");
650 
651  // check that target is valid w.r.t. dims
652  if (!internal::check_subsys_match_dims(target, dims))
653  throw exception::SubsysMismatchDims("qpp::apply()");
654 
655  // check valid state and matching dimensions
656  if (internal::check_cvector(rA)) {
657  if (!internal::check_dims_match_cvect(dims, rA))
658  throw exception::DimsMismatchCvector("qpp::apply()");
659  } else if (internal::check_square_mat(rA)) {
660  if (!internal::check_dims_match_mat(dims, rA))
661  throw exception::DimsMismatchMatrix("qpp::apply()");
662  } else
663  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
664 
665  std::vector<idx> subsys_dims(target.size());
666  for (idx i = 0; i < target.size(); ++i)
667  subsys_dims[i] = dims[target[i]];
668 
669  // check the Kraus operators
670  if (Ks.empty())
671  throw exception::ZeroSize("qpp::apply()");
672  if (!internal::check_square_mat(Ks[0]))
673  throw exception::MatrixNotSquare("qpp::apply()");
674  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
675  throw exception::MatrixMismatchSubsys("qpp::apply()");
676  for (auto&& elem : Ks)
677  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
678  throw exception::DimsNotEqual("qpp::apply()");
679  // END EXCEPTION CHECKS
680 
681  cmat result = cmat::Zero(rA.rows(), rA.rows());
682 
683  for (idx i = 0; i < Ks.size(); ++i)
684  result += apply(rA, Ks[i], target, dims);
685 
686  return result;
687 }
688 
699 template <typename Derived>
700 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
701  const std::vector<idx>& target, idx d = 2) {
702  const cmat& rA = A.derived();
703 
704  // EXCEPTION CHECKS
705 
706  // check zero sizes
707  if (!internal::check_nonzero_size(rA))
708  throw exception::ZeroSize("qpp::apply()");
709 
710  // check valid dims
711  if (d < 2)
712  throw exception::DimsInvalid("qpp::apply()");
713  // END EXCEPTION CHECKS
714 
715  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
716  std::vector<idx> dims(n, d); // local dimensions vector
717 
718  return apply(rA, Ks, target, dims);
719 }
720 
732 inline cmat kraus2super(const std::vector<cmat>& Ks) {
733  // EXCEPTION CHECKS
734 
735  if (Ks.empty())
736  throw exception::ZeroSize("qpp::kraus2super()");
737  if (!internal::check_nonzero_size(Ks[0]))
738  throw exception::ZeroSize("qpp::kraus2super()");
739  if (!internal::check_square_mat(Ks[0]))
740  throw exception::MatrixNotSquare("qpp::kraus2super()");
741  for (auto&& elem : Ks)
742  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
743  throw exception::DimsNotEqual("qpp::kraus2super()");
744  // END EXCEPTION CHECKS
745 
746  idx D = static_cast<idx>(Ks[0].rows());
747 
748  cmat result(D * D, D * D);
749  cmat MN = cmat::Zero(D, D);
750  bra A = bra::Zero(D);
751  ket B = ket::Zero(D);
752  cmat EMN = cmat::Zero(D, D);
753 
754 #ifdef WITH_OPENMP_
755 #pragma omp parallel for collapse(2)
756 #endif // WITH_OPENMP_
757  for (idx m = 0; m < D; ++m) {
758  for (idx n = 0; n < D; ++n) {
759 #ifdef WITH_OPENMP_
760 #pragma omp critical
761 #endif // WITH_OPENMP_
762  {
763  // compute E(|m><n|)
764  MN(m, n) = 1;
765  for (idx i = 0; i < Ks.size(); ++i)
766  EMN += Ks[i] * MN * adjoint(Ks[i]);
767  MN(m, n) = 0;
768 
769  for (idx a = 0; a < D; ++a) {
770  A(a) = 1;
771  for (idx b = 0; b < D; ++b) {
772  // compute result(ab,mn)=<a|E(|m><n)|b>
773  B(b) = 1;
774  result(a * D + b, m * D + n) =
775  static_cast<cmat>(A * EMN * B).value();
776  B(b) = 0;
777  }
778  A(a) = 0;
779  }
780  EMN = cmat::Zero(D, D);
781  }
782  }
783  }
784 
785  return result;
786 }
787 
803 inline cmat kraus2choi(const std::vector<cmat>& Ks) {
804  // EXCEPTION CHECKS
805 
806  if (Ks.empty())
807  throw exception::ZeroSize("qpp::kraus2choi()");
808  if (!internal::check_nonzero_size(Ks[0]))
809  throw exception::ZeroSize("qpp::kraus2choi()");
810  if (!internal::check_square_mat(Ks[0]))
811  throw exception::MatrixNotSquare("qpp::kraus2choi()");
812  for (auto&& elem : Ks)
813  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
814  throw exception::DimsNotEqual("qpp::kraus2choi()");
815  // END EXCEPTION CHECKS
816 
817  idx D = static_cast<idx>(Ks[0].rows());
818 
819  // construct the D x D \sum |jj> vector
820  // (un-normalized maximally entangled state)
821  cmat MES = cmat::Zero(D * D, 1);
822  for (idx a = 0; a < D; ++a)
823  MES(a * D + a) = 1;
824 
825  cmat Omega = MES * adjoint(MES);
826 
827  cmat result = cmat::Zero(D * D, D * D);
828 
829 #ifdef WITH_OPENMP_
830 #pragma omp parallel for
831 #endif // WITH_OPENMP_
832  for (idx i = 0; i < Ks.size(); ++i) {
833 #ifdef WITH_OPENMP_
834 #pragma omp critical
835 #endif // WITH_OPENMP_
836  {
837  result += kron(cmat::Identity(D, D), Ks[i]) * Omega *
838  adjoint(kron(cmat::Identity(D, D), Ks[i]));
839  }
840  }
841 
842  return result;
843 }
844 
858 inline std::vector<cmat> choi2kraus(const cmat& A) {
859  // EXCEPTION CHECKS
860 
861  if (!internal::check_nonzero_size(A))
862  throw exception::ZeroSize("qpp::choi2kraus()");
863  if (!internal::check_square_mat(A))
864  throw exception::MatrixNotSquare("qpp::choi2kraus()");
865  idx D = internal::get_dim_subsys(A.rows(), 2);
866  // check equal dimensions
867  if (D * D != static_cast<idx>(A.rows()))
868  throw exception::DimsInvalid("qpp::choi2kraus()");
869  // END EXCEPTION CHECKS
870 
871  dmat ev = hevals(A);
872  cmat evec = hevects(A);
873  std::vector<cmat> result;
874 
875  for (idx i = 0; i < D * D; ++i) {
876  if (std::abs(ev(i)) > 0)
877  result.emplace_back(std::sqrt(std::abs(ev(i))) *
878  reshape(evec.col(i), D, D));
879  }
880 
881  return result;
882 }
883 
891 inline cmat choi2super(const cmat& A) {
892  // EXCEPTION CHECKS
893 
894  if (!internal::check_nonzero_size(A))
895  throw exception::ZeroSize("qpp::choi2super()");
896  if (!internal::check_square_mat(A))
897  throw exception::MatrixNotSquare("qpp::choi2super()");
898  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
899  // check equal dimensions
900  if (D * D != static_cast<idx>(A.rows()))
901  throw exception::DimsInvalid("qpp::choi2super()");
902  // END EXCEPTION CHECKS
903 
904  cmat result(D * D, D * D);
905 
906 #ifdef WITH_OPENMP_
907 #pragma omp parallel for collapse(4)
908 #endif // WITH_OPENMP_
909  for (idx a = 0; a < D; ++a)
910  for (idx b = 0; b < D; ++b)
911  for (idx m = 0; m < D; ++m)
912  for (idx n = 0; n < D; ++n)
913  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
914 
915  return result;
916 }
917 
925 inline cmat super2choi(const cmat& A) {
926  // EXCEPTION CHECKS
927 
928  if (!internal::check_nonzero_size(A))
929  throw exception::ZeroSize("qpp::super2choi()");
930  if (!internal::check_square_mat(A))
931  throw exception::MatrixNotSquare("qpp::super2choi()");
932  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
933  // check equal dimensions
934  if (D * D != static_cast<idx>(A.rows()))
935  throw exception::DimsInvalid("qpp::super2choi()");
936  // END EXCEPTION CHECKS
937 
938  cmat result(D * D, D * D);
939 
940 #ifdef WITH_OPENMP_
941 #pragma omp parallel for collapse(4)
942 #endif // WITH_OPENMP_
943  for (idx a = 0; a < D; ++a)
944  for (idx b = 0; b < D; ++b)
945  for (idx m = 0; m < D; ++m)
946  for (idx n = 0; n < D; ++n)
947  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
948 
949  return result;
950 }
951 
965 template <typename Derived>
966 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
967  const std::vector<idx>& dims) {
968  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
969 
970  // EXCEPTION CHECKS
971 
972  // check zero-size
973  if (!internal::check_nonzero_size(rA))
974  throw exception::ZeroSize("qpp::ptrace1()");
975 
976  // check that dims is a valid dimension vector
977  if (!internal::check_dims(dims))
978  throw exception::DimsInvalid("qpp::ptrace1()");
979 
980  // check dims has only 2 elements
981  if (dims.size() != 2)
982  throw exception::NotBipartite("qpp::ptrace1()");
983  // END EXCEPTION CHECKS
984 
985  idx DA = dims[0];
986  idx DB = dims[1];
987 
988  dyn_mat<typename Derived::Scalar> result =
989  dyn_mat<typename Derived::Scalar>::Zero(DB, DB);
990 
991  //************ ket ************//
992  if (internal::check_cvector(rA)) // we have a ket
993  {
994  // check that dims match the dimension of A
995  if (!internal::check_dims_match_cvect(dims, rA))
996  throw exception::DimsMismatchCvector("qpp::ptrace1()");
997 
998  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
999  typename Derived::Scalar sum = 0;
1000  for (idx m = 0; m < DA; ++m)
1001  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
1002 
1003  return sum;
1004  }; /* end worker */
1005 
1006 #ifdef WITH_OPENMP_
1007 #pragma omp parallel for collapse(2)
1008 #endif // WITH_OPENMP_
1009  // column major order for speed
1010  for (idx j = 0; j < DB; ++j)
1011  for (idx i = 0; i < DB; ++i)
1012  result(i, j) = worker(i, j);
1013  }
1014  //************ density matrix ************//
1015  else if (internal::check_square_mat(rA)) // we have a density operator
1016  {
1017  // check that dims match the dimension of A
1018  if (!internal::check_dims_match_mat(dims, rA))
1019  throw exception::DimsMismatchMatrix("qpp::ptrace1()");
1020 
1021  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
1022  typename Derived::Scalar sum = 0;
1023  for (idx m = 0; m < DA; ++m)
1024  sum += rA(m * DB + i, m * DB + j);
1025 
1026  return sum;
1027  }; /* end worker */
1028 
1029 #ifdef WITH_OPENMP_
1030 #pragma omp parallel for collapse(2)
1031 #endif // WITH_OPENMP_
1032  // column major order for speed
1033  for (idx j = 0; j < DB; ++j)
1034  for (idx i = 0; i < DB; ++i)
1035  result(i, j) = worker(i, j);
1036  }
1037  //************ Exception: not ket nor density matrix ************//
1038  else
1039  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1040 
1041  return result;
1042 }
1043 
1057 template <typename Derived>
1058 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1059  idx d = 2) {
1060  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1061 
1062  // EXCEPTION CHECKS
1063 
1064  // check zero size
1065  if (!internal::check_nonzero_size(rA))
1066  throw exception::ZeroSize("qpp::ptrace1()");
1067 
1068  // check valid dims
1069  if (d == 0)
1070  throw exception::DimsInvalid("qpp::ptrace1()");
1071  // END EXCEPTION CHECKS
1072 
1073  std::vector<idx> dims(2, d); // local dimensions vector
1074 
1075  return ptrace1(rA, dims);
1076 }
1077 
1091 template <typename Derived>
1092 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1093  const std::vector<idx>& dims) {
1094  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1095 
1096  // EXCEPTION CHECKS
1097 
1098  // check zero-size
1099  if (!internal::check_nonzero_size(rA))
1100  throw exception::ZeroSize("qpp::ptrace2()");
1101 
1102  // check that dims is a valid dimension vector
1103  if (!internal::check_dims(dims))
1104  throw exception::DimsInvalid("qpp::ptrace2()");
1105 
1106  // check dims has only 2 elements
1107  if (dims.size() != 2)
1108  throw exception::NotBipartite("qpp::ptrace2()");
1109  // END EXCEPTION CHECKS
1110 
1111  idx DA = dims[0];
1112  idx DB = dims[1];
1113 
1114  dyn_mat<typename Derived::Scalar> result =
1115  dyn_mat<typename Derived::Scalar>::Zero(DA, DA);
1116 
1117  //************ ket ************//
1118  if (internal::check_cvector(rA)) // we have a ket
1119  {
1120  // check that dims match the dimension of A
1121  if (!internal::check_dims_match_cvect(dims, rA))
1122  throw exception::DimsMismatchCvector("qpp::ptrace2()");
1123 
1124  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
1125  typename Derived::Scalar sum = 0;
1126  for (idx m = 0; m < DB; ++m)
1127  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1128 
1129  return sum;
1130  }; /* end worker */
1131 
1132 #ifdef WITH_OPENMP_
1133 #pragma omp parallel for collapse(2)
1134 #endif // WITH_OPENMP_
1135  // column major order for speed
1136  for (idx j = 0; j < DA; ++j)
1137  for (idx i = 0; i < DA; ++i)
1138  result(i, j) = worker(i, j);
1139  }
1140  //************ density matrix ************//
1141  else if (internal::check_square_mat(rA)) // we have a density operator
1142  {
1143  // check that dims match the dimension of A
1144  if (!internal::check_dims_match_mat(dims, rA))
1145  throw exception::DimsMismatchMatrix("qpp::ptrace2()");
1146 
1147 #ifdef WITH_OPENMP_
1148 #pragma omp parallel for collapse(2)
1149 #endif // WITH_OPENMP_
1150  // column major order for speed
1151  for (idx j = 0; j < DA; ++j)
1152  for (idx i = 0; i < DA; ++i)
1153  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1154  }
1155  //************ Exception: not ket nor density matrix ************//
1156  else
1157  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1158 
1159  return result;
1160 }
1161 
1175 template <typename Derived>
1176 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1177  idx d = 2) {
1178  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1179 
1180  // EXCEPTION CHECKS
1181 
1182  // check zero size
1183  if (!internal::check_nonzero_size(rA))
1184  throw exception::ZeroSize("qpp::ptrace2()");
1185 
1186  // check valid dims
1187  if (d == 0)
1188  throw exception::DimsInvalid("qpp::ptrace2()");
1189  // END EXCEPTION CHECKS
1190 
1191  std::vector<idx> dims(2, d); // local dimensions vector
1192 
1193  return ptrace2(rA, dims);
1194 }
1195 
1210 template <typename Derived>
1211 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1212  const std::vector<idx>& target,
1213  const std::vector<idx>& dims) {
1214  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1215 
1216  // EXCEPTION CHECKS
1217 
1218  // check zero-size
1219  if (!internal::check_nonzero_size(rA))
1220  throw exception::ZeroSize("qpp::ptrace()");
1221 
1222  // check that dims is a valid dimension vector
1223  if (!internal::check_dims(dims))
1224  throw exception::DimsInvalid("qpp::ptrace()");
1225 
1226  // check that target is valid w.r.t. dims
1227  if (!internal::check_subsys_match_dims(target, dims))
1228  throw exception::SubsysMismatchDims("qpp::ptrace()");
1229 
1230  // check valid state and matching dimensions
1231  if (internal::check_cvector(rA)) {
1232  if (!internal::check_dims_match_cvect(dims, rA))
1233  throw exception::DimsMismatchCvector("qpp::ptrace()");
1234  } else if (internal::check_square_mat(rA)) {
1235  if (!internal::check_dims_match_mat(dims, rA))
1236  throw exception::DimsMismatchMatrix("qpp::ptrace()");
1237  } else
1238  throw exception::MatrixNotSquareNorCvector("qpp::ptrace()");
1239  // END EXCEPTION CHECKS
1240 
1241  idx D = static_cast<idx>(rA.rows());
1242  idx n = dims.size();
1243  idx n_subsys = target.size();
1244  idx n_subsys_bar = n - n_subsys;
1245  idx Dsubsys = 1;
1246  for (idx i = 0; i < n_subsys; ++i)
1247  Dsubsys *= dims[target[i]];
1248  idx Dsubsys_bar = D / Dsubsys;
1249 
1250  idx Cdims[maxn];
1251  idx Csubsys[maxn];
1252  idx Cdimssubsys[maxn];
1253  idx Csubsys_bar[maxn];
1254  idx Cdimssubsys_bar[maxn];
1255 
1256  idx Cmidxcolsubsys_bar[maxn];
1257 
1258  std::vector<idx> subsys_bar = complement(target, n);
1259  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1260  std::begin(Csubsys_bar));
1261 
1262  for (idx i = 0; i < n; ++i) {
1263  Cdims[i] = dims[i];
1264  }
1265  for (idx i = 0; i < n_subsys; ++i) {
1266  Csubsys[i] = target[i];
1267  Cdimssubsys[i] = dims[target[i]];
1268  }
1269  for (idx i = 0; i < n_subsys_bar; ++i) {
1270  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1271  }
1272 
1273  dyn_mat<typename Derived::Scalar> result =
1274  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1275 
1276  //************ ket ************//
1277  if (internal::check_cvector(rA)) // we have a ket
1278  {
1279  if (target.size() == dims.size()) {
1280  result(0, 0) = (adjoint(rA) * rA).value();
1281  return result;
1282  }
1283 
1284  if (target.empty())
1285  return rA * adjoint(rA);
1286 
1287  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1288  // use static allocation for speed!
1289 
1290  idx Cmidxrow[maxn];
1291  idx Cmidxcol[maxn];
1292  idx Cmidxrowsubsys_bar[maxn];
1293  idx Cmidxsubsys[maxn];
1294 
1295  /* get the row multi-indexes of the complement */
1296  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1297  Cmidxrowsubsys_bar);
1298  /* write them in the global row/col multi-indexes */
1299  for (idx k = 0; k < n_subsys_bar; ++k) {
1300  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1301  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1302  }
1303  typename Derived::Scalar sm = 0;
1304  for (idx a = 0; a < Dsubsys; ++a) {
1305  // get the multi-index over which we do the summation
1306  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1307  // write it into the global row/col multi-indexes
1308  for (idx k = 0; k < n_subsys; ++k)
1309  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1310  Cmidxsubsys[k];
1311 
1312  // now do the sum
1313  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims)) *
1314  std::conj(rA(internal::multiidx2n(Cmidxcol, n, Cdims)));
1315  }
1316 
1317  return sm;
1318  }; /* end worker */
1319 
1320  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1321  {
1322  // compute the column multi-indexes of the complement
1323  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1324  Cmidxcolsubsys_bar);
1325 #ifdef WITH_OPENMP_
1326 #pragma omp parallel for
1327 #endif // WITH_OPENMP_
1328  for (idx i = 0; i < Dsubsys_bar; ++i) {
1329  result(i, j) = worker(i);
1330  }
1331  }
1332  }
1333  //************ density matrix ************//
1334  else // we have a density operator
1335  {
1336  if (target.size() == dims.size()) {
1337  result(0, 0) = rA.trace();
1338  return result;
1339  }
1340 
1341  if (target.empty())
1342  return rA;
1343 
1344  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1345  // use static allocation for speed!
1346 
1347  idx Cmidxrow[maxn];
1348  idx Cmidxcol[maxn];
1349  idx Cmidxrowsubsys_bar[maxn];
1350  idx Cmidxsubsys[maxn];
1351 
1352  /* get the row/col multi-indexes of the complement */
1353  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1354  Cmidxrowsubsys_bar);
1355  /* write them in the global row/col multi-indexes */
1356  for (idx k = 0; k < n_subsys_bar; ++k) {
1357  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1358  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1359  }
1360  typename Derived::Scalar sm = 0;
1361  for (idx a = 0; a < Dsubsys; ++a) {
1362  // get the multi-index over which we do the summation
1363  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1364  // write it into the global row/col multi-indexes
1365  for (idx k = 0; k < n_subsys; ++k)
1366  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1367  Cmidxsubsys[k];
1368 
1369  // now do the sum
1370  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims),
1371  internal::multiidx2n(Cmidxcol, n, Cdims));
1372  }
1373 
1374  return sm;
1375  }; /* end worker */
1376 
1377  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1378  {
1379  // compute the column multi-indexes of the complement
1380  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1381  Cmidxcolsubsys_bar);
1382 #ifdef WITH_OPENMP_
1383 #pragma omp parallel for
1384 #endif // WITH_OPENMP_
1385  for (idx i = 0; i < Dsubsys_bar; ++i) {
1386  result(i, j) = worker(i);
1387  }
1388  }
1389  }
1390 
1391  return result;
1392 }
1393 
1408 template <typename Derived>
1409 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1410  const std::vector<idx>& target,
1411  idx d = 2) {
1412  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1413 
1414  // EXCEPTION CHECKS
1415 
1416  // check zero size
1417  if (!internal::check_nonzero_size(rA))
1418  throw exception::ZeroSize("qpp::ptrace()");
1419 
1420  // check valid dims
1421  if (d < 2)
1422  throw exception::DimsInvalid("qpp::ptrace()");
1423  // END EXCEPTION CHECKS
1424 
1425  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1426  std::vector<idx> dims(n, d); // local dimensions vector
1427 
1428  return ptrace(rA, target, dims);
1429 }
1430 
1444 template <typename Derived>
1445 dyn_mat<typename Derived::Scalar>
1446 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& target,
1447  const std::vector<idx>& dims) {
1448  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1449 
1450  // EXCEPTION CHECKS
1451 
1452  // check zero-size
1453  if (!internal::check_nonzero_size(rA))
1454  throw exception::ZeroSize("qpp::ptranspose()");
1455 
1456  // check that dims is a valid dimension vector
1457  if (!internal::check_dims(dims))
1458  throw exception::DimsInvalid("qpp::ptranspose()");
1459 
1460  // check that target is valid w.r.t. dims
1461  if (!internal::check_subsys_match_dims(target, dims))
1462  throw exception::SubsysMismatchDims("qpp::ptranspose()");
1463 
1464  // check valid state and matching dimensions
1465  if (internal::check_cvector(rA)) {
1466  if (!internal::check_dims_match_cvect(dims, rA))
1467  throw exception::DimsMismatchCvector("qpp::ptranspose()");
1468  } else if (internal::check_square_mat(rA)) {
1469  if (!internal::check_dims_match_mat(dims, rA))
1470  throw exception::DimsMismatchMatrix("qpp::ptranspose()");
1471  } else
1472  throw exception::MatrixNotSquareNorCvector("qpp::ptranspose()");
1473  // END EXCEPTION CHECKS
1474 
1475  idx D = static_cast<idx>(rA.rows());
1476  idx n = dims.size();
1477  idx n_subsys = target.size();
1478  idx Cdims[maxn];
1479  idx Cmidxcol[maxn];
1480  idx Csubsys[maxn];
1481 
1482  // copy dims in Cdims and target in Csubsys
1483  for (idx i = 0; i < n; ++i)
1484  Cdims[i] = dims[i];
1485  for (idx i = 0; i < n_subsys; ++i)
1486  Csubsys[i] = target[i];
1487 
1488  dyn_mat<typename Derived::Scalar> result(D, D);
1489 
1490  //************ ket ************//
1491  if (internal::check_cvector(rA)) // we have a ket
1492  {
1493  if (target.size() == dims.size())
1494  return (rA * adjoint(rA)).transpose();
1495 
1496  if (target.empty())
1497  return rA * adjoint(rA);
1498 
1499  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1500  // use static allocation for speed!
1501  idx midxcoltmp[maxn];
1502  idx midxrow[maxn];
1503 
1504  for (idx k = 0; k < n; ++k)
1505  midxcoltmp[k] = Cmidxcol[k];
1506 
1507  /* compute the row multi-index */
1508  internal::n2multiidx(i, n, Cdims, midxrow);
1509 
1510  for (idx k = 0; k < n_subsys; ++k)
1511  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1512 
1513  /* writes the result */
1514  return rA(internal::multiidx2n(midxrow, n, Cdims)) *
1515  std::conj(rA(internal::multiidx2n(midxcoltmp, n, Cdims)));
1516  }; /* end worker */
1517 
1518  for (idx j = 0; j < D; ++j) {
1519  // compute the column multi-index
1520  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1521 
1522 #ifdef WITH_OPENMP_
1523 #pragma omp parallel for
1524 #endif // WITH_OPENMP_
1525  for (idx i = 0; i < D; ++i)
1526  result(i, j) = worker(i);
1527  }
1528  }
1529  //************ density matrix ************//
1530  else // we have a density operator
1531  {
1532  if (target.size() == dims.size())
1533  return rA.transpose();
1534 
1535  if (target.empty())
1536  return rA;
1537 
1538  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1539  // use static allocation for speed!
1540  idx midxcoltmp[maxn];
1541  idx midxrow[maxn];
1542 
1543  for (idx k = 0; k < n; ++k)
1544  midxcoltmp[k] = Cmidxcol[k];
1545 
1546  /* compute the row multi-index */
1547  internal::n2multiidx(i, n, Cdims, midxrow);
1548 
1549  for (idx k = 0; k < n_subsys; ++k)
1550  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1551 
1552  /* writes the result */
1553  return rA(internal::multiidx2n(midxrow, n, Cdims),
1554  internal::multiidx2n(midxcoltmp, n, Cdims));
1555  }; /* end worker */
1556 
1557  for (idx j = 0; j < D; ++j) {
1558  // compute the column multi-index
1559  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1560 
1561 #ifdef WITH_OPENMP_
1562 #pragma omp parallel for
1563 #endif // WITH_OPENMP_
1564  for (idx i = 0; i < D; ++i)
1565  result(i, j) = worker(i);
1566  }
1567  }
1568 
1569  return result;
1570 }
1571 
1585 template <typename Derived>
1586 dyn_mat<typename Derived::Scalar>
1587 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& target,
1588  idx d = 2) {
1589  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1590 
1591  // EXCEPTION CHECKS
1592 
1593  // check zero size
1594  if (!internal::check_nonzero_size(rA))
1595  throw exception::ZeroSize("qpp::ptranspose()");
1596 
1597  // check valid dims
1598  if (d < 2)
1599  throw exception::DimsInvalid("qpp::ptranspose()");
1600  // END EXCEPTION CHECKS
1601 
1602  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1603  std::vector<idx> dims(n, d); // local dimensions vector
1604 
1605  return ptranspose(rA, target, dims);
1606 }
1607 
1620 template <typename Derived>
1621 dyn_mat<typename Derived::Scalar>
1622 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1623  const std::vector<idx>& dims) {
1624  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1625 
1626  // EXCEPTION CHECKS
1627 
1628  // check zero-size
1629  if (!internal::check_nonzero_size(rA))
1630  throw exception::ZeroSize("qpp::syspermute()");
1631 
1632  // check that dims is a valid dimension vector
1633  if (!internal::check_dims(dims))
1634  throw exception::DimsInvalid("qpp::syspermute()");
1635 
1636  // check that we have a valid permutation
1637  if (!internal::check_perm(perm))
1638  throw exception::PermInvalid("qpp::syspermute()");
1639 
1640  // check that permutation match dimensions
1641  if (perm.size() != dims.size())
1642  throw exception::PermMismatchDims("qpp::syspermute()");
1643 
1644  // check valid state and matching dimensions
1645  if (internal::check_cvector(rA)) {
1646  if (!internal::check_dims_match_cvect(dims, rA))
1647  throw exception::DimsMismatchCvector("qpp::syspermute()");
1648  } else if (internal::check_square_mat(rA)) {
1649  if (!internal::check_dims_match_mat(dims, rA))
1650  throw exception::DimsMismatchMatrix("qpp::syspermute()");
1651  } else
1652  throw exception::MatrixNotSquareNorCvector("qpp::syspermute()");
1653  // END EXCEPTION CHECKS
1654 
1655  idx D = static_cast<idx>(rA.rows());
1656  idx n = dims.size();
1657 
1658  dyn_mat<typename Derived::Scalar> result;
1659 
1660  //************ ket ************//
1661  if (internal::check_cvector(rA)) // we have a column vector
1662  {
1663  idx Cdims[maxn];
1664  idx Cperm[maxn];
1665 
1666  // copy dims in Cdims and perm in Cperm
1667  for (idx i = 0; i < n; ++i) {
1668  Cdims[i] = dims[i];
1669  Cperm[i] = perm[i];
1670  }
1671  result.resize(D, 1);
1672 
1673  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1674  // use static allocation for speed,
1675  // double the size for matrices reshaped as vectors
1676  idx midx[maxn];
1677  idx midxtmp[maxn];
1678  idx permdims[maxn];
1679 
1680  /* compute the multi-index */
1681  internal::n2multiidx(i, n, Cdims, midx);
1682 
1683  for (idx k = 0; k < n; ++k) {
1684  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1685  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1686  }
1687 
1688  return internal::multiidx2n(midxtmp, n, permdims);
1689  }; /* end worker */
1690 
1691 #ifdef WITH_OPENMP_
1692 #pragma omp parallel for
1693 #endif // WITH_OPENMP_
1694  for (idx i = 0; i < D; ++i)
1695  result(worker(i)) = rA(i);
1696  }
1697  //************ density matrix ************//
1698  else // we have a density operator
1699  {
1700  idx Cdims[2 * maxn];
1701  idx Cperm[2 * maxn];
1702 
1703  // copy dims in Cdims and perm in Cperm
1704  for (idx i = 0; i < n; ++i) {
1705  Cdims[i] = dims[i];
1706  Cdims[i + n] = dims[i];
1707  Cperm[i] = perm[i];
1708  Cperm[i + n] = perm[i] + n;
1709  }
1710  result.resize(D * D, 1);
1711  // map A to a column vector
1712  dyn_mat<typename Derived::Scalar> vectA =
1713  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1714  const_cast<typename Derived::Scalar*>(rA.data()), D * D, 1);
1715 
1716  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1717  // use static allocation for speed,
1718  // double the size for matrices reshaped as vectors
1719  idx midx[2 * maxn];
1720  idx midxtmp[2 * maxn];
1721  idx permdims[2 * maxn];
1722 
1723  /* compute the multi-index */
1724  internal::n2multiidx(i, 2 * n, Cdims, midx);
1725 
1726  for (idx k = 0; k < 2 * n; ++k) {
1727  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1728  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1729  }
1730 
1731  return internal::multiidx2n(midxtmp, 2 * n, permdims);
1732  }; /* end worker */
1733 
1734 #ifdef WITH_OPENMP_
1735 #pragma omp parallel for
1736 #endif // WITH_OPENMP_
1737  for (idx i = 0; i < D * D; ++i)
1738  result(worker(i)) = rA(i);
1739 
1740  result = reshape(result, D, D);
1741  }
1742 
1743  return result;
1744 }
1745 
1758 template <typename Derived>
1759 dyn_mat<typename Derived::Scalar>
1760 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1761  idx d = 2) {
1762  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1763 
1764  // EXCEPTION CHECKS
1765 
1766  // check zero size
1767  if (!internal::check_nonzero_size(rA))
1768  throw exception::ZeroSize("qpp::syspermute()");
1769 
1770  // check valid dims
1771  if (d < 2)
1772  throw exception::DimsInvalid("qpp::syspermute()");
1773  // END EXCEPTION CHECKS
1774 
1775  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1776  std::vector<idx> dims(n, d); // local dimensions vector
1777 
1778  return syspermute(rA, perm, dims);
1779 }
1780 
1781 // as in https://arxiv.org/abs/1707.08834
1792 template <typename Derived>
1793 dyn_mat<typename Derived::Scalar> applyQFT(const Eigen::MatrixBase<Derived>& A,
1794  const std::vector<idx>& target,
1795  idx d = 2, bool swap = true) {
1796  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1797 
1798  // EXCEPTION CHECKS
1799 
1800  // check zero sizes
1801  if (!internal::check_nonzero_size(rA))
1802  throw exception::ZeroSize("qpp::applyQFT()");
1803 
1804  // check valid subsystem dimension
1805  if (d < 2)
1806  throw exception::DimsInvalid("qpp::applyQFT()");
1807 
1808  // total number of qubits/qudits in the state
1809  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1810 
1811  std::vector<idx> dims(n, d); // local dimensions vector
1812 
1813  // check that target is valid w.r.t. dims
1814  if (!internal::check_subsys_match_dims(target, dims))
1815  throw exception::SubsysMismatchDims("qpp::applyQFT()");
1816 
1817  // check valid state and matching dimensions
1818  if (internal::check_cvector(rA)) {
1819  if (!internal::check_dims_match_cvect(dims, rA))
1820  throw exception::DimsMismatchCvector("qpp::applyQFT()");
1821  } else if (internal::check_square_mat(rA)) {
1822  if (!internal::check_dims_match_mat(dims, rA))
1823  throw exception::DimsMismatchMatrix("qpp::applyQFT()");
1824  } else
1825  throw exception::MatrixNotSquareNorCvector("qpp::applyQFT()");
1826  // END EXCEPTION CHECKS
1827 
1828  dyn_mat<typename Derived::Scalar> result = rA;
1829 
1830  idx n_subsys = target.size();
1831 
1832  if (d == 2) // qubits
1833  {
1834  for (idx i = 0; i < n_subsys; ++i) {
1835  // apply Hadamard on qubit i
1836  result = apply(result, Gates::get_instance().H, {target[i]});
1837  // apply controlled rotations
1838  for (idx j = 2; j <= n_subsys - i; ++j) {
1839  // construct Rj
1840  cmat Rj(2, 2);
1841  Rj << 1, 0, 0, exp(2.0 * pi * 1_i / std::pow(2, j));
1842  result =
1843  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1844  }
1845  }
1846  if (swap) {
1847  // we have the qubits in reversed order, we must swap them
1848  for (idx i = 0; i < n_subsys / 2; ++i) {
1849  result = apply(result, Gates::get_instance().SWAP,
1850  {target[i], target[n_subsys - i - 1]});
1851  }
1852  }
1853 
1854  } else { // qudits
1855  for (idx i = 0; i < n_subsys; ++i) {
1856  // apply qudit Fourier on qudit i
1857  result = apply(result, Gates::get_instance().Fd(d), {target[i]}, d);
1858  // apply controlled rotations
1859  for (idx j = 2; j <= n_subsys - i; ++j) {
1860  // construct Rj
1861  cmat Rj = cmat::Zero(d, d);
1862  for (idx m = 0; m < d; ++m) {
1863  Rj(m, m) = exp(2.0 * pi * m * 1_i / std::pow(d, j));
1864  }
1865  result =
1866  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1867  }
1868  }
1869  if (swap) {
1870  // we have the qudits in reversed order, we must swap them
1871  for (idx i = 0; i < n_subsys / 2; ++i) {
1872  result = apply(result, Gates::get_instance().SWAPd(d),
1873  {target[i], target[n_subsys - i - 1]}, d);
1874  }
1875  }
1876  }
1877 
1878  return result;
1879 }
1880 
1881 // as in https://arxiv.org/abs/1707.08834
1893 template <typename Derived>
1894 dyn_mat<typename Derived::Scalar> applyTFQ(const Eigen::MatrixBase<Derived>& A,
1895  const std::vector<idx>& target,
1896  idx d = 2, bool swap = true) {
1897  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1898 
1899  // EXCEPTION CHECKS
1900 
1901  // check zero sizes
1902  if (!internal::check_nonzero_size(rA))
1903  throw exception::ZeroSize("qpp::applyTFQ()");
1904 
1905  // check valid subsystem dimension
1906  if (d < 2)
1907  throw exception::DimsInvalid("qpp::applyTFQ()");
1908 
1909  // total number of qubits/qudits in the state
1910  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1911 
1912  std::vector<idx> dims(n, d); // local dimensions vector
1913 
1914  // check that target is valid w.r.t. dims
1915  if (!internal::check_subsys_match_dims(target, dims))
1916  throw exception::SubsysMismatchDims("qpp::applyTFQ()");
1917 
1918  // check valid state and matching dimensions
1919  if (internal::check_cvector(rA)) {
1920  if (!internal::check_dims_match_cvect(dims, rA))
1921  throw exception::DimsMismatchCvector("qpp::applyTFQ()");
1922  } else if (internal::check_square_mat(rA)) {
1923  if (!internal::check_dims_match_mat(dims, rA))
1924  throw exception::DimsMismatchMatrix("qpp::applyTFQ()");
1925  } else
1926  throw exception::MatrixNotSquareNorCvector("qpp::applyTFQ()");
1927  // END EXCEPTION CHECKS
1928 
1929  dyn_mat<typename Derived::Scalar> result = rA;
1930 
1931  idx n_subsys = target.size();
1932 
1933  if (d == 2) // qubits
1934  {
1935  if (swap) {
1936  // we have the qubits in reversed order, we must swap them
1937  for (idx i = n_subsys / 2; i-- > 0;) {
1938  result = apply(result, Gates::get_instance().SWAP,
1939  {target[i], target[n_subsys - i - 1]});
1940  }
1941  }
1942  for (idx i = n_subsys; i-- > 0;) {
1943  // apply controlled rotations
1944  for (idx j = n_subsys - i + 1; j-- > 2;) {
1945  // construct Rj
1946  cmat Rj(2, 2);
1947  Rj << 1, 0, 0, exp(-2.0 * pi * 1_i / std::pow(2, j));
1948  result =
1949  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1950  }
1951  // apply Hadamard on qubit i
1952  result = apply(result, Gates::get_instance().H, {target[i]});
1953  }
1954  } else { // qudits
1955  if (swap) {
1956  // we have the qudits in reversed order, we must swap them
1957  for (idx i = n_subsys / 2; i-- > 0;) {
1958  result = apply(result, Gates::get_instance().SWAPd(d),
1959  {target[i], target[n_subsys - i - 1]}, d);
1960  }
1961  }
1962  for (idx i = n_subsys; i-- > 0;) {
1963  // apply controlled rotations
1964  for (idx j = n_subsys - i + 1; j-- > 2;) {
1965  // construct Rj
1966  cmat Rj = cmat::Zero(d, d);
1967  for (idx m = 0; m < d; ++m) {
1968  Rj(m, m) = exp(-2.0 * pi * m * 1_i / std::pow(d, j));
1969  }
1970  result =
1971  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1972  }
1973  // apply qudit Fourier on qudit i
1974  result = apply(result, adjoint(Gates::get_instance().Fd(d)),
1975  {target[i]}, d);
1976  }
1977  }
1978 
1979  return result;
1980 }
1981 
1982 // as in https://arxiv.org/abs/1707.08834
1991 template <typename Derived>
1992 dyn_col_vect<typename Derived::Scalar> QFT(const Eigen::MatrixBase<Derived>& A,
1993  idx d = 2, bool swap = true) {
1994  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1995 
1996  // EXCEPTION CHECKS
1997 
1998  // check zero-size
1999  if (!internal::check_nonzero_size(rA))
2000  throw exception::ZeroSize("qpp::QFT()");
2001 
2002  // check valid subsystem dimension
2003  if (d < 2)
2004  throw exception::DimsInvalid("qpp::QFT()");
2005 
2006  // total number of qubits/qudits in the state
2007  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
2008 
2009  std::vector<idx> dims(n, d); // local dimensions vector
2010 
2011  // check valid state and matching dimensions
2012  if (internal::check_cvector(rA)) {
2013  if (!internal::check_dims_match_cvect(dims, rA))
2014  throw exception::DimsMismatchCvector("qpp::QFT()");
2015  } else if (internal::check_square_mat(rA)) {
2016  if (!internal::check_dims_match_mat(dims, rA))
2017  throw exception::DimsMismatchMatrix("qpp::QFT()");
2018  } else
2019  throw exception::MatrixNotSquareNorCvector("qpp::QFT()");
2020  // END EXCEPTION CHECKS
2021 
2022  std::vector<idx> subsys(n);
2023  std::iota(std::begin(subsys), std::end(subsys), 0);
2024  ket result = applyQFT(rA, subsys, d, swap);
2025 
2026  return result;
2027 }
2028 
2029 // as in https://arxiv.org/abs/1707.08834
2038 template <typename Derived>
2039 dyn_col_vect<typename Derived::Scalar> TFQ(const Eigen::MatrixBase<Derived>& A,
2040  idx d = 2, bool swap = true) {
2041  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2042 
2043  // EXCEPTION CHECKS
2044 
2045  // check zero-size
2046  if (!internal::check_nonzero_size(rA))
2047  throw exception::ZeroSize("qpp::TFQ()");
2048 
2049  // check valid subsystem dimension
2050  if (d < 2)
2051  throw exception::DimsInvalid("qpp::TFQ()");
2052 
2053  // total number of qubits/qudits in the state
2054  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
2055 
2056  std::vector<idx> dims(n, d); // local dimensions vector
2057 
2058  // check valid state and matching dimensions
2059  if (internal::check_cvector(rA)) {
2060  if (!internal::check_dims_match_cvect(dims, rA))
2061  throw exception::DimsMismatchCvector("qpp::QFT()");
2062  } else if (internal::check_square_mat(rA)) {
2063  if (!internal::check_dims_match_mat(dims, rA))
2064  throw exception::DimsMismatchMatrix("qpp::QFT()");
2065  } else
2066  throw exception::MatrixNotSquareNorCvector("qpp::QFT()");
2067  // END EXCEPTION CHECKS
2068 
2069  std::vector<idx> subsys(n);
2070  std::iota(std::begin(subsys), std::end(subsys), 0);
2071  ket result = applyTFQ(rA, subsys, d, swap);
2072 
2073  return result;
2074 }
2075 
2076 } /* namespace qpp */
2077 
2078 #endif /* OPERATIONS_H_ */
qpp
Quantum++ main namespace.
Definition: circuits.h:35
Generated by   doxygen 1.8.11