Quantum++  v1.2
A modern C++11 quantum computing library
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 - 2019 Vlad Gheorghiu (vgheorgh@gmail.com)
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 
53 template <typename Derived1, typename Derived2>
54 dyn_mat<typename Derived1::Scalar>
55 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
56  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
57  const std::vector<idx>& target, const std::vector<idx>& dims) {
58  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
59  state.derived();
60  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
61 
62  // EXCEPTION CHECKS
63 
64  // check types
65  if (!std::is_same<typename Derived1::Scalar,
66  typename Derived2::Scalar>::value)
67  throw exception::TypeMismatch("qpp::applyCTRL()");
68 
69  // check zero sizes
71  throw exception::ZeroSize("qpp::applyCTRL()");
72 
73  // check zero sizes
74  if (!internal::check_nonzero_size(rstate))
75  throw exception::ZeroSize("qpp::applyCTRL()");
76 
77  // check zero sizes
78  if (!internal::check_nonzero_size(target))
79  throw exception::ZeroSize("qpp::applyCTRL()");
80 
81  // check square matrix for the gate
83  throw exception::MatrixNotSquare("qpp::applyCTRL()");
84 
85  // check that all control subsystems have the same dimension
86  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
87  for (idx i = 1; i < ctrl.size(); ++i)
88  if (dims[ctrl[i]] != d)
89  throw exception::DimsNotEqual("qpp::applyCTRL()");
90 
91  // check that dimension is valid
92  if (!internal::check_dims(dims))
93  throw exception::DimsInvalid("qpp::applyCTRL()");
94 
95  // check that target is valid w.r.t. dims
96  if (!internal::check_subsys_match_dims(target, dims))
97  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
98 
99  // check that gate matches the dimensions of the target
100  std::vector<idx> target_dims(target.size());
101  for (idx i = 0; i < target.size(); ++i)
102  target_dims[i] = dims[target[i]];
103  if (!internal::check_dims_match_mat(target_dims, rA))
104  throw exception::MatrixMismatchSubsys("qpp::applyCTRL()");
105 
106  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
107  ctrlgate.insert(std::end(ctrlgate), std::begin(target), std::end(target));
108  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
109 
110  // check that ctrl + gate subsystem is valid
111  // with respect to local dimensions
112  if (!internal::check_subsys_match_dims(ctrlgate, dims))
113  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
114  // END EXCEPTION CHECKS
115 
116  // construct the table of A^i and (A^dagger)^i
117  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
118  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
119  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i) {
120  Ai.push_back(powm(rA, i));
121  Aidagger.push_back(powm(adjoint(rA), i));
122  }
123 
124  idx D = static_cast<idx>(rstate.rows()); // total dimension
125  idx n = dims.size(); // total number of subsystems
126  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
127  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
128  idx targetsize = target.size(); // number of subsystems of the target
129  // dimension of ctrl subsystem
130  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
131  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
132 
133  idx Cdims[maxn]; // local dimensions
134  idx CdimsA[maxn]; // local dimensions
135  idx CdimsCTRL[maxn]; // local dimensions
136  idx CdimsCTRLA_bar[maxn]; // local dimensions
137 
138  // compute the complementary subsystem of ctrlgate w.r.t. dims
139  std::vector<idx> ctrlgate_bar = complement(ctrlgate, n);
140  // number of subsystems that are complementary to the ctrl+gate
141  idx ctrlgate_barsize = ctrlgate_bar.size();
142 
143  idx DCTRLA_bar = 1; // dimension of the rest
144  for (idx i = 0; i < ctrlgate_barsize; ++i)
145  DCTRLA_bar *= dims[ctrlgate_bar[i]];
146 
147  for (idx k = 0; k < n; ++k)
148  Cdims[k] = dims[k];
149  for (idx k = 0; k < targetsize; ++k)
150  CdimsA[k] = dims[target[k]];
151  for (idx k = 0; k < ctrlsize; ++k)
152  CdimsCTRL[k] = d;
153  for (idx k = 0; k < ctrlgate_barsize; ++k)
154  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
155 
156  // worker, computes the coefficient and the index for the ket case
157  // used in #pragma omp parallel for collapse
158  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
159  ->std::pair<typename Derived1::Scalar, idx> {
160  idx indx = 0;
161  typename Derived1::Scalar coeff = 0;
162 
163  idx Cmidx[maxn]; // the total multi-index
164  idx CmidxA[maxn]; // the gate part multi-index
165  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
166 
167  // compute the index
168 
169  // set the CTRL part
170  for (idx k = 0; k < ctrlsize; ++k) {
171  Cmidx[ctrl[k]] = i_;
172  }
173 
174  // set the rest
175  internal::n2multiidx(r_, n - ctrlgatesize, CdimsCTRLA_bar,
176  CmidxCTRLA_bar);
177  for (idx k = 0; k < n - ctrlgatesize; ++k) {
178  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
179  }
180 
181  // set the A part
182  internal::n2multiidx(m_, targetsize, CdimsA, CmidxA);
183  for (idx k = 0; k < targetsize; ++k) {
184  Cmidx[target[k]] = CmidxA[k];
185  }
186 
187  // we now got the total index
188  indx = internal::multiidx2n(Cmidx, n, Cdims);
189 
190  // compute the coefficient
191  for (idx n_ = 0; n_ < DA; ++n_) {
192  internal::n2multiidx(n_, targetsize, CdimsA, CmidxA);
193  for (idx k = 0; k < targetsize; ++k) {
194  Cmidx[target[k]] = CmidxA[k];
195  }
196  coeff +=
197  Ai[i_](m_, n_) * rstate(internal::multiidx2n(Cmidx, n, Cdims));
198  }
199 
200  return std::make_pair(coeff, indx);
201  }; /* end coeff_idx_ket */
202 
203  // worker, computes the coefficient and the index
204  // for the density matrix case
205  // used in #pragma omp parallel for collapse
206  auto coeff_idx_rho = [&](idx i1_, idx m1_, idx r1_, idx i2_, idx m2_,
207  idx r2_) noexcept
208  ->std::tuple<typename Derived1::Scalar, idx, idx> {
209  idx idxrow = 0;
210  idx idxcol = 0;
211  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
212 
213  idx Cmidxrow[maxn]; // the total row multi-index
214  idx Cmidxcol[maxn]; // the total col multi-index
215  idx CmidxArow[maxn]; // the gate part row multi-index
216  idx CmidxAcol[maxn]; // the gate part col multi-index
217  idx CmidxCTRLrow[maxn]; // the control row multi-index
218  idx CmidxCTRLcol[maxn]; // the control col multi-index
219  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
220  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
221 
222  // compute the ket/bra indexes
223 
224  // set the CTRL part
225  internal::n2multiidx(i1_, ctrlsize, CdimsCTRL, CmidxCTRLrow);
226  internal::n2multiidx(i2_, ctrlsize, CdimsCTRL, CmidxCTRLcol);
227 
228  for (idx k = 0; k < ctrlsize; ++k) {
229  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
230  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
231  }
232 
233  // set the rest
234  internal::n2multiidx(r1_, n - ctrlgatesize, CdimsCTRLA_bar,
235  CmidxCTRLA_barrow);
236  internal::n2multiidx(r2_, n - ctrlgatesize, CdimsCTRLA_bar,
237  CmidxCTRLA_barcol);
238  for (idx k = 0; k < n - ctrlgatesize; ++k) {
239  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
240  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
241  }
242 
243  // set the A part
244  internal::n2multiidx(m1_, targetsize, CdimsA, CmidxArow);
245  internal::n2multiidx(m2_, targetsize, CdimsA, CmidxAcol);
246  for (idx k = 0; k < target.size(); ++k) {
247  Cmidxrow[target[k]] = CmidxArow[k];
248  Cmidxcol[target[k]] = CmidxAcol[k];
249  }
250 
251  // we now got the total row/col indexes
252  idxrow = internal::multiidx2n(Cmidxrow, n, Cdims);
253  idxcol = internal::multiidx2n(Cmidxcol, n, Cdims);
254 
255  // check whether all CTRL row and col multi indexes are equal
256  bool all_ctrl_rows_equal = true;
257  bool all_ctrl_cols_equal = true;
258 
259  idx first_ctrl_row, first_ctrl_col;
260  if (ctrlsize > 0) {
261  first_ctrl_row = CmidxCTRLrow[0];
262  first_ctrl_col = CmidxCTRLcol[0];
263  } else {
264  first_ctrl_row = first_ctrl_col = 1;
265  }
266 
267  for (idx k = 1; k < ctrlsize; ++k) {
268  if (CmidxCTRLrow[k] != first_ctrl_row) {
269  all_ctrl_rows_equal = false;
270  break;
271  }
272  }
273  for (idx k = 1; k < ctrlsize; ++k) {
274  if (CmidxCTRLcol[k] != first_ctrl_col) {
275  all_ctrl_cols_equal = false;
276  break;
277  }
278  }
279 
280  // at least one control activated, compute the coefficient
281  for (idx n1_ = 0; n1_ < DA; ++n1_) {
282  internal::n2multiidx(n1_, targetsize, CdimsA, CmidxArow);
283  for (idx k = 0; k < targetsize; ++k) {
284  Cmidxrow[target[k]] = CmidxArow[k];
285  }
286  idx idxrowtmp = internal::multiidx2n(Cmidxrow, n, Cdims);
287 
288  if (all_ctrl_rows_equal) {
289  lhs = Ai[first_ctrl_row](m1_, n1_);
290  } else {
291  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
292  }
293 
294  for (idx n2_ = 0; n2_ < DA; ++n2_) {
295  internal::n2multiidx(n2_, targetsize, CdimsA, CmidxAcol);
296  for (idx k = 0; k < targetsize; ++k) {
297  Cmidxcol[target[k]] = CmidxAcol[k];
298  }
299 
300  if (all_ctrl_cols_equal) {
301  rhs = Aidagger[first_ctrl_col](n2_, m2_);
302  } else {
303  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
304  }
305 
306  idx idxcoltmp = internal::multiidx2n(Cmidxcol, n, Cdims);
307 
308  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
309  }
310  }
311 
312  return std::make_tuple(coeff, idxrow, idxcol);
313  }; /* end coeff_idx_rho */
314 
315  //************ ket ************//
316  if (internal::check_cvector(rstate)) // we have a ket
317  {
318  // check that dims match state vector
319  if (!internal::check_dims_match_cvect(dims, rstate))
320  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
321  if (D == 1)
322  return rstate;
323 
324  dyn_mat<typename Derived1::Scalar> result = rstate;
325 
326 #ifdef WITH_OPENMP_
327 #pragma omp parallel for collapse(2)
328 #endif // WITH_OPENMP_
329  for (idx m = 0; m < DA; ++m)
330  for (idx r = 0; r < DCTRLA_bar; ++r) {
331  if (ctrlsize == 0) // no control
332  {
333  result(coeff_idx_ket(1, m, r).second) =
334  coeff_idx_ket(1, m, r).first;
335  } else
336  for (idx i = 0; i < d; ++i) {
337  result(coeff_idx_ket(i, m, r).second) =
338  coeff_idx_ket(i, m, r).first;
339  }
340  }
341 
342  return result;
343  }
344  //************ density matrix ************//
345  else if (internal::check_square_mat(rstate)) // we have a density operator
346  {
347  // check that dims match state matrix
348  if (!internal::check_dims_match_mat(dims, rstate))
349  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
350 
351  if (D == 1)
352  return rstate;
353 
354  dyn_mat<typename Derived1::Scalar> result = rstate;
355 
356 #ifdef WITH_OPENMP_
357 #pragma omp parallel for collapse(4)
358 #endif // WITH_OPENMP_
359  for (idx m1 = 0; m1 < DA; ++m1)
360  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
361  for (idx m2 = 0; m2 < DA; ++m2)
362  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
363  if (ctrlsize == 0) // no control
364  {
365  auto coeff_idxes =
366  coeff_idx_rho(1, m1, r1, 1, m2, r2);
367  result(std::get<1>(coeff_idxes),
368  std::get<2>(coeff_idxes)) =
369  std::get<0>(coeff_idxes);
370  } else {
371  for (idx i1 = 0; i1 < Dctrl; ++i1)
372  for (idx i2 = 0; i2 < Dctrl; ++i2) {
373  auto coeff_idxes =
374  coeff_idx_rho(i1, m1, r1, i2, m2, r2);
375  result(std::get<1>(coeff_idxes),
376  std::get<2>(coeff_idxes)) =
377  std::get<0>(coeff_idxes);
378  }
379  }
380 
381  return result;
382  }
383  //************ Exception: not ket nor density matrix ************//
384  else
385  throw exception::MatrixNotSquareNorCvector("qpp::applyCTRL()");
386 }
387 
403 template <typename Derived1, typename Derived2>
404 dyn_mat<typename Derived1::Scalar>
405 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
406  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
407  const std::vector<idx>& target, idx d = 2) {
408  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
409  state.derived();
410  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
411 
412  // EXCEPTION CHECKS
413 
414  // check zero size
415  if (!internal::check_nonzero_size(rstate))
416  throw exception::ZeroSize("qpp::applyCTRL()");
417 
418  // check valid dims
419  if (d < 2)
420  throw exception::DimsInvalid("qpp::applyCTRL()");
421  // END EXCEPTION CHECKS
422 
423  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
424  std::vector<idx> dims(n, d); // local dimensions vector
425 
426  return applyCTRL(rstate, rA, ctrl, target, dims);
427 }
428 
442 template <typename Derived1, typename Derived2>
443 dyn_mat<typename Derived1::Scalar>
444 apply(const Eigen::MatrixBase<Derived1>& state,
445  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& target,
446  const std::vector<idx>& dims) {
447  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
448  state.derived();
449  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
450 
451  // EXCEPTION CHECKS
452 
453  // check types
454  if (!std::is_same<typename Derived1::Scalar,
455  typename Derived2::Scalar>::value)
456  throw exception::TypeMismatch("qpp::apply()");
457 
458  // check zero sizes
460  throw exception::ZeroSize("qpp::apply()");
461 
462  // check zero sizes
463  if (!internal::check_nonzero_size(rstate))
464  throw exception::ZeroSize("qpp::apply()");
465 
466  // check zero sizes
467  if (!internal::check_nonzero_size(target))
468  throw exception::ZeroSize("qpp::apply()");
469 
470  // check square matrix for the gate
472  throw exception::MatrixNotSquare("qpp::apply()");
473 
474  // check that dimension is valid
475  if (!internal::check_dims(dims))
476  throw exception::DimsInvalid("qpp::apply()");
477 
478  // check that target is valid w.r.t. dims
479  if (!internal::check_subsys_match_dims(target, dims))
480  throw exception::SubsysMismatchDims("qpp::apply()");
481 
482  // check that gate matches the dimensions of the target
483  std::vector<idx> subsys_dims(target.size());
484  for (idx i = 0; i < target.size(); ++i)
485  subsys_dims[i] = dims[target[i]];
486  if (!internal::check_dims_match_mat(subsys_dims, rA))
487  throw exception::MatrixMismatchSubsys("qpp::apply()");
488  // END EXCEPTION CHECKS
489 
490  //************ ket ************//
491  if (internal::check_cvector(rstate)) // we have a ket
492  {
493  // check that dims match state vector
494  if (!internal::check_dims_match_cvect(dims, rstate))
495  throw exception::DimsMismatchCvector("qpp::apply()");
496 
497  return applyCTRL(rstate, rA, {}, target, dims);
498  }
499  //************ density matrix ************//
500  else if (internal::check_square_mat(rstate)) // we have a density operator
501  {
502  // check that dims match state matrix
503  if (!internal::check_dims_match_mat(dims, rstate))
504  throw exception::DimsMismatchMatrix("qpp::apply()");
505 
506  return applyCTRL(rstate, rA, {}, target, dims);
507  }
508  //************ Exception: not ket nor density matrix ************//
509  else
510  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
511 }
512 
526 template <typename Derived1, typename Derived2>
527 dyn_mat<typename Derived1::Scalar>
528 apply(const Eigen::MatrixBase<Derived1>& state,
529  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& target,
530  idx d = 2) {
531  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
532  state.derived();
533  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
534 
535  // EXCEPTION CHECKS
536 
537  // check zero size
538  if (!internal::check_nonzero_size(rstate))
539  throw exception::ZeroSize("qpp::apply()");
540 
541  // check valid dims
542  if (d < 2)
543  throw exception::DimsInvalid("qpp::apply()");
544  // END EXCEPTION CHECKS
545 
546  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
547  std::vector<idx> dims(n, d); // local dimensions vector
548 
549  return apply(rstate, rA, target, dims);
550 }
551 
560 template <typename Derived>
561 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
562  const cmat& rA = A.derived();
563 
564  // EXCEPTION CHECKS
565 
567  throw exception::ZeroSize("qpp::apply()");
569  throw exception::MatrixNotSquare("qpp::apply()");
570  if (Ks.size() == 0)
571  throw exception::ZeroSize("qpp::apply()");
572  if (!internal::check_square_mat(Ks[0]))
573  throw exception::MatrixNotSquare("qpp::apply()");
574  if (Ks[0].rows() != rA.rows())
575  throw exception::DimsMismatchMatrix("qpp::apply()");
576  for (auto&& elem : Ks)
577  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
578  throw exception::DimsNotEqual("qpp::apply()");
579  // END EXCEPTION CHECKS
580 
581  cmat result = cmat::Zero(rA.rows(), rA.rows());
582 
583 #ifdef WITH_OPENMP_
584 #pragma omp parallel for
585 #endif // WITH_OPENMP_
586  for (idx i = 0; i < Ks.size(); ++i) {
587 #ifdef WITH_OPENMP_
588 #pragma omp critical
589 #endif // WITH_OPENMP_
590  { result += Ks[i] * rA * adjoint(Ks[i]); }
591  }
592 
593  return result;
594 }
595 
606 template <typename Derived>
607 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
608  const std::vector<idx>& target, const std::vector<idx>& dims) {
609  const cmat& rA = A.derived();
610 
611  // EXCEPTION CHECKS
612 
613  // check zero sizes
615  throw exception::ZeroSize("qpp::apply()");
616 
617  // check zero sizes
618  if (!internal::check_nonzero_size(target))
619  throw exception::ZeroSize("qpp::apply()");
620 
621  // check square matrix for the A
623  throw exception::MatrixNotSquare("qpp::apply()");
624 
625  // check that dimension is valid
626  if (!internal::check_dims(dims))
627  throw exception::DimsInvalid("qpp::apply()");
628 
629  // check that dims match A matrix
630  if (!internal::check_dims_match_mat(dims, rA))
631  throw exception::DimsMismatchMatrix("qpp::apply()");
632 
633  // check that target is valid w.r.t. dims
634  if (!internal::check_subsys_match_dims(target, dims))
635  throw exception::SubsysMismatchDims("qpp::apply()");
636 
637  std::vector<idx> subsys_dims(target.size());
638  for (idx i = 0; i < target.size(); ++i)
639  subsys_dims[i] = dims[target[i]];
640 
641  // check the Kraus operators
642  if (Ks.size() == 0)
643  throw exception::ZeroSize("qpp::apply()");
644  if (!internal::check_square_mat(Ks[0]))
645  throw exception::MatrixNotSquare("qpp::apply()");
646  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
647  throw exception::MatrixMismatchSubsys("qpp::apply()");
648  for (auto&& elem : Ks)
649  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
650  throw exception::DimsNotEqual("qpp::apply()");
651  // END EXCEPTION CHECKS
652 
653  cmat result = cmat::Zero(rA.rows(), rA.rows());
654 
655  for (idx i = 0; i < Ks.size(); ++i)
656  result += apply(rA, Ks[i], target, dims);
657 
658  return result;
659 }
660 
671 template <typename Derived>
672 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
673  const std::vector<idx>& target, idx d = 2) {
674  const cmat& rA = A.derived();
675 
676  // EXCEPTION CHECKS
677 
678  // check zero sizes
680  throw exception::ZeroSize("qpp::apply()");
681 
682  // check valid dims
683  if (d < 2)
684  throw exception::DimsInvalid("qpp::apply()");
685  // END EXCEPTION CHECKS
686 
687  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
688  std::vector<idx> dims(n, d); // local dimensions vector
689 
690  return apply(rA, Ks, target, dims);
691 }
692 
704 inline cmat kraus2super(const std::vector<cmat>& Ks) {
705  // EXCEPTION CHECKS
706 
707  if (Ks.size() == 0)
708  throw exception::ZeroSize("qpp::kraus2super()");
709  if (!internal::check_nonzero_size(Ks[0]))
710  throw exception::ZeroSize("qpp::kraus2super()");
711  if (!internal::check_square_mat(Ks[0]))
712  throw exception::MatrixNotSquare("qpp::kraus2super()");
713  for (auto&& elem : Ks)
714  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
715  throw exception::DimsNotEqual("qpp::kraus2super()");
716  // END EXCEPTION CHECKS
717 
718  idx D = static_cast<idx>(Ks[0].rows());
719 
720  cmat result(D * D, D * D);
721  cmat MN = cmat::Zero(D, D);
722  bra A = bra::Zero(D);
723  ket B = ket::Zero(D);
724  cmat EMN = cmat::Zero(D, D);
725 
726 #ifdef WITH_OPENMP_
727 #pragma omp parallel for collapse(2)
728 #endif // WITH_OPENMP_
729  for (idx m = 0; m < D; ++m) {
730  for (idx n = 0; n < D; ++n) {
731 #ifdef WITH_OPENMP_
732 #pragma omp critical
733 #endif // WITH_OPENMP_
734  {
735  // compute E(|m><n|)
736  MN(m, n) = 1;
737  for (idx i = 0; i < Ks.size(); ++i)
738  EMN += Ks[i] * MN * adjoint(Ks[i]);
739  MN(m, n) = 0;
740 
741  for (idx a = 0; a < D; ++a) {
742  A(a) = 1;
743  for (idx b = 0; b < D; ++b) {
744  // compute result(ab,mn)=<a|E(|m><n)|b>
745  B(b) = 1;
746  result(a * D + b, m * D + n) =
747  static_cast<cmat>(A * EMN * B).value();
748  B(b) = 0;
749  }
750  A(a) = 0;
751  }
752  EMN = cmat::Zero(D, D);
753  }
754  }
755  }
756 
757  return result;
758 }
759 
775 inline cmat kraus2choi(const std::vector<cmat>& Ks) {
776  // EXCEPTION CHECKS
777 
778  if (Ks.size() == 0)
779  throw exception::ZeroSize("qpp::kraus2choi()");
780  if (!internal::check_nonzero_size(Ks[0]))
781  throw exception::ZeroSize("qpp::kraus2choi()");
782  if (!internal::check_square_mat(Ks[0]))
783  throw exception::MatrixNotSquare("qpp::kraus2choi()");
784  for (auto&& elem : Ks)
785  if (elem.rows() != Ks[0].rows() || elem.cols() != Ks[0].rows())
786  throw exception::DimsNotEqual("qpp::kraus2choi()");
787  // END EXCEPTION CHECKS
788 
789  idx D = static_cast<idx>(Ks[0].rows());
790 
791  // construct the D x D \sum |jj> vector
792  // (un-normalized maximally entangled state)
793  cmat MES = cmat::Zero(D * D, 1);
794  for (idx a = 0; a < D; ++a)
795  MES(a * D + a) = 1;
796 
797  cmat Omega = MES * adjoint(MES);
798 
799  cmat result = cmat::Zero(D * D, D * D);
800 
801 #ifdef WITH_OPENMP_
802 #pragma omp parallel for
803 #endif // WITH_OPENMP_
804  for (idx i = 0; i < Ks.size(); ++i) {
805 #ifdef WITH_OPENMP_
806 #pragma omp critical
807 #endif // WITH_OPENMP_
808  {
809  result += kron(cmat::Identity(D, D), Ks[i]) * Omega *
810  adjoint(kron(cmat::Identity(D, D), Ks[i]));
811  }
812  }
813 
814  return result;
815 }
816 
830 inline std::vector<cmat> choi2kraus(const cmat& A) {
831  // EXCEPTION CHECKS
832 
834  throw exception::ZeroSize("qpp::choi2kraus()");
836  throw exception::MatrixNotSquare("qpp::choi2kraus()");
837  idx D = internal::get_dim_subsys(A.rows(), 2);
838  // check equal dimensions
839  if (D * D != static_cast<idx>(A.rows()))
840  throw exception::DimsInvalid("qpp::choi2kraus()");
841  // END EXCEPTION CHECKS
842 
843  dmat ev = hevals(A);
844  cmat evec = hevects(A);
845  std::vector<cmat> result;
846 
847  for (idx i = 0; i < D * D; ++i) {
848  if (std::abs(ev(i)) > 0)
849  result.push_back(std::sqrt(std::abs(ev(i))) *
850  reshape(evec.col(i), D, D));
851  }
852 
853  return result;
854 }
855 
863 inline cmat choi2super(const cmat& A) {
864  // EXCEPTION CHECKS
865 
867  throw exception::ZeroSize("qpp::choi2super()");
869  throw exception::MatrixNotSquare("qpp::choi2super()");
870  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
871  // check equal dimensions
872  if (D * D != static_cast<idx>(A.rows()))
873  throw exception::DimsInvalid("qpp::choi2super()");
874  // END EXCEPTION CHECKS
875 
876  cmat result(D * D, D * D);
877 
878 #ifdef WITH_OPENMP_
879 #pragma omp parallel for collapse(4)
880 #endif // WITH_OPENMP_
881  for (idx a = 0; a < D; ++a)
882  for (idx b = 0; b < D; ++b)
883  for (idx m = 0; m < D; ++m)
884  for (idx n = 0; n < D; ++n)
885  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
886 
887  return result;
888 }
889 
897 inline cmat super2choi(const cmat& A) {
898  // EXCEPTION CHECKS
899 
901  throw exception::ZeroSize("qpp::super2choi()");
903  throw exception::MatrixNotSquare("qpp::super2choi()");
904  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
905  // check equal dimensions
906  if (D * D != static_cast<idx>(A.rows()))
907  throw exception::DimsInvalid("qpp::super2choi()");
908  // END EXCEPTION CHECKS
909 
910  cmat result(D * D, D * D);
911 
912 #ifdef WITH_OPENMP_
913 #pragma omp parallel for collapse(4)
914 #endif // WITH_OPENMP_
915  for (idx a = 0; a < D; ++a)
916  for (idx b = 0; b < D; ++b)
917  for (idx m = 0; m < D; ++m)
918  for (idx n = 0; n < D; ++n)
919  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
920 
921  return result;
922 }
923 
937 template <typename Derived>
938 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
939  const std::vector<idx>& dims) {
940  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
941 
942  // EXCEPTION CHECKS
943 
944  // check zero-size
946  throw exception::ZeroSize("qpp::ptrace1()");
947 
948  // check that dims is a valid dimension vector
949  if (!internal::check_dims(dims))
950  throw exception::DimsInvalid("qpp::ptrace1()");
951 
952  // check dims has only 2 elements
953  if (dims.size() != 2)
954  throw exception::NotBipartite("qpp::ptrace1()");
955  // END EXCEPTION CHECKS
956 
957  idx DA = dims[0];
958  idx DB = dims[1];
959 
962 
963  //************ ket ************//
964  if (internal::check_cvector(rA)) // we have a ket
965  {
966  // check that dims match the dimension of A
967  if (!internal::check_dims_match_cvect(dims, rA))
968  throw exception::DimsMismatchCvector("qpp::ptrace1()");
969 
970  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
971  typename Derived::Scalar sum = 0;
972  for (idx m = 0; m < DA; ++m)
973  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
974 
975  return sum;
976  }; /* end worker */
977 
978 #ifdef WITH_OPENMP_
979 #pragma omp parallel for collapse(2)
980 #endif // WITH_OPENMP_
981  // column major order for speed
982  for (idx j = 0; j < DB; ++j)
983  for (idx i = 0; i < DB; ++i)
984  result(i, j) = worker(i, j);
985 
986  return result;
987  }
988  //************ density matrix ************//
989  else if (internal::check_square_mat(rA)) // we have a density operator
990  {
991  // check that dims match the dimension of A
992  if (!internal::check_dims_match_mat(dims, rA))
993  throw exception::DimsMismatchMatrix("qpp::ptrace1()");
994 
995  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
996  typename Derived::Scalar sum = 0;
997  for (idx m = 0; m < DA; ++m)
998  sum += rA(m * DB + i, m * DB + j);
999 
1000  return sum;
1001  }; /* end worker */
1002 
1003 #ifdef WITH_OPENMP_
1004 #pragma omp parallel for collapse(2)
1005 #endif // WITH_OPENMP_
1006  // column major order for speed
1007  for (idx j = 0; j < DB; ++j)
1008  for (idx i = 0; i < DB; ++i)
1009  result(i, j) = worker(i, j);
1010 
1011  return result;
1012  }
1013  //************ Exception: not ket nor density matrix ************//
1014  else
1015  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1016 }
1017 
1031 template <typename Derived>
1032 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1033  idx d = 2) {
1034  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1035 
1036  // EXCEPTION CHECKS
1037 
1038  // check zero size
1040  throw exception::ZeroSize("qpp::ptrace1()");
1041 
1042  // check valid dims
1043  if (d == 0)
1044  throw exception::DimsInvalid("qpp::ptrace1()");
1045  // END EXCEPTION CHECKS
1046 
1047  std::vector<idx> dims(2, d); // local dimensions vector
1048 
1049  return ptrace1(A, dims);
1050 }
1051 
1065 template <typename Derived>
1066 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1067  const std::vector<idx>& dims) {
1068  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1069 
1070  // EXCEPTION CHECKS
1071 
1072  // check zero-size
1074  throw exception::ZeroSize("qpp::ptrace2()");
1075 
1076  // check that dims is a valid dimension vector
1077  if (!internal::check_dims(dims))
1078  throw exception::DimsInvalid("qpp::ptrace2()");
1079 
1080  // check dims has only 2 elements
1081  if (dims.size() != 2)
1082  throw exception::NotBipartite("qpp::ptrace2()");
1083  // END EXCEPTION CHECKS
1084 
1085  idx DA = dims[0];
1086  idx DB = dims[1];
1087 
1090 
1091  //************ ket ************//
1092  if (internal::check_cvector(rA)) // we have a ket
1093  {
1094  // check that dims match the dimension of A
1095  if (!internal::check_dims_match_cvect(dims, rA))
1096  throw exception::DimsMismatchCvector("qpp::ptrace2()");
1097 
1098  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
1099  typename Derived::Scalar sum = 0;
1100  for (idx m = 0; m < DB; ++m)
1101  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1102 
1103  return sum;
1104  }; /* end worker */
1105 
1106 #ifdef WITH_OPENMP_
1107 #pragma omp parallel for collapse(2)
1108 #endif // WITH_OPENMP_
1109  // column major order for speed
1110  for (idx j = 0; j < DA; ++j)
1111  for (idx i = 0; i < DA; ++i)
1112  result(i, j) = worker(i, j);
1113 
1114  return result;
1115  }
1116  //************ density matrix ************//
1117  else if (internal::check_square_mat(rA)) // we have a density operator
1118  {
1119  // check that dims match the dimension of A
1120  if (!internal::check_dims_match_mat(dims, rA))
1121  throw exception::DimsMismatchMatrix("qpp::ptrace2()");
1122 
1123 #ifdef WITH_OPENMP_
1124 #pragma omp parallel for collapse(2)
1125 #endif // WITH_OPENMP_
1126  // column major order for speed
1127  for (idx j = 0; j < DA; ++j)
1128  for (idx i = 0; i < DA; ++i)
1129  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1130 
1131  return result;
1132  }
1133  //************ Exception: not ket nor density matrix ************//
1134  else
1135  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1136 }
1137 
1151 template <typename Derived>
1152 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1153  idx d = 2) {
1154  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1155 
1156  // EXCEPTION CHECKS
1157 
1158  // check zero size
1160  throw exception::ZeroSize("qpp::ptrace2()");
1161 
1162  // check valid dims
1163  if (d == 0)
1164  throw exception::DimsInvalid("qpp::ptrace2()");
1165  // END EXCEPTION CHECKS
1166 
1167  std::vector<idx> dims(2, d); // local dimensions vector
1168 
1169  return ptrace2(A, dims);
1170 }
1171 
1186 template <typename Derived>
1187 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1188  const std::vector<idx>& target,
1189  const std::vector<idx>& dims) {
1190  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1191 
1192  // EXCEPTION CHECKS
1193 
1194  // check zero-size
1196  throw exception::ZeroSize("qpp::ptrace()");
1197 
1198  // check that dims is a valid dimension vector
1199  if (!internal::check_dims(dims))
1200  throw exception::DimsInvalid("qpp::ptrace()");
1201 
1202  // check that target are valid
1203  if (!internal::check_subsys_match_dims(target, dims))
1204  throw exception::SubsysMismatchDims("qpp::ptrace()");
1205  // END EXCEPTION CHECKS
1206 
1207  idx D = static_cast<idx>(rA.rows());
1208  idx n = dims.size();
1209  idx n_subsys = target.size();
1210  idx n_subsys_bar = n - n_subsys;
1211  idx Dsubsys = 1;
1212  for (idx i = 0; i < n_subsys; ++i)
1213  Dsubsys *= dims[target[i]];
1214  idx Dsubsys_bar = D / Dsubsys;
1215 
1216  idx Cdims[maxn];
1217  idx Csubsys[maxn];
1218  idx Cdimssubsys[maxn];
1219  idx Csubsys_bar[maxn];
1220  idx Cdimssubsys_bar[maxn];
1221 
1222  idx Cmidxcolsubsys_bar[maxn];
1223 
1224  std::vector<idx> subsys_bar = complement(target, n);
1225  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1226  std::begin(Csubsys_bar));
1227 
1228  for (idx i = 0; i < n; ++i) {
1229  Cdims[i] = dims[i];
1230  }
1231  for (idx i = 0; i < n_subsys; ++i) {
1232  Csubsys[i] = target[i];
1233  Cdimssubsys[i] = dims[target[i]];
1234  }
1235  for (idx i = 0; i < n_subsys_bar; ++i) {
1236  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1237  }
1238 
1240  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1241 
1242  //************ ket ************//
1243  if (internal::check_cvector(rA)) // we have a ket
1244  {
1245  // check that dims match the dimension of A
1246  if (!internal::check_dims_match_cvect(dims, rA))
1247  throw exception::DimsMismatchCvector("qpp::ptrace()");
1248 
1249  if (target.size() == dims.size()) {
1250  result(0, 0) = (adjoint(rA) * rA).value();
1251  return result;
1252  }
1253 
1254  if (target.size() == 0)
1255  return rA * adjoint(rA);
1256 
1257  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1258  // use static allocation for speed!
1259 
1260  idx Cmidxrow[maxn];
1261  idx Cmidxcol[maxn];
1262  idx Cmidxrowsubsys_bar[maxn];
1263  idx Cmidxsubsys[maxn];
1264 
1265  /* get the row multi-indexes of the complement */
1266  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1267  Cmidxrowsubsys_bar);
1268  /* write them in the global row/col multi-indexes */
1269  for (idx k = 0; k < n_subsys_bar; ++k) {
1270  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1271  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1272  }
1273  typename Derived::Scalar sm = 0;
1274  for (idx a = 0; a < Dsubsys; ++a) {
1275  // get the multi-index over which we do the summation
1276  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1277  // write it into the global row/col multi-indexes
1278  for (idx k = 0; k < n_subsys; ++k)
1279  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1280  Cmidxsubsys[k];
1281 
1282  // now do the sum
1283  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims)) *
1284  std::conj(rA(internal::multiidx2n(Cmidxcol, n, Cdims)));
1285  }
1286 
1287  return sm;
1288  }; /* end worker */
1289 
1290  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1291  {
1292  // compute the column multi-indexes of the complement
1293  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1294  Cmidxcolsubsys_bar);
1295 #ifdef WITH_OPENMP_
1296 #pragma omp parallel for
1297 #endif // WITH_OPENMP_
1298  for (idx i = 0; i < Dsubsys_bar; ++i) {
1299  result(i, j) = worker(i);
1300  }
1301  }
1302 
1303  return result;
1304  }
1305  //************ density matrix ************//
1306  else if (internal::check_square_mat(rA)) // we have a density operator
1307  {
1308  // check that dims match the dimension of A
1309  if (!internal::check_dims_match_mat(dims, rA))
1310  throw exception::DimsMismatchMatrix("qpp::ptrace()");
1311 
1312  if (target.size() == dims.size()) {
1313  result(0, 0) = rA.trace();
1314  return result;
1315  }
1316 
1317  if (target.size() == 0)
1318  return rA;
1319 
1320  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1321  // use static allocation for speed!
1322 
1323  idx Cmidxrow[maxn];
1324  idx Cmidxcol[maxn];
1325  idx Cmidxrowsubsys_bar[maxn];
1326  idx Cmidxsubsys[maxn];
1327 
1328  /* get the row/col multi-indexes of the complement */
1329  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1330  Cmidxrowsubsys_bar);
1331  /* write them in the global row/col multi-indexes */
1332  for (idx k = 0; k < n_subsys_bar; ++k) {
1333  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1334  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1335  }
1336  typename Derived::Scalar sm = 0;
1337  for (idx a = 0; a < Dsubsys; ++a) {
1338  // get the multi-index over which we do the summation
1339  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1340  // write it into the global row/col multi-indexes
1341  for (idx k = 0; k < n_subsys; ++k)
1342  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1343  Cmidxsubsys[k];
1344 
1345  // now do the sum
1346  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims),
1347  internal::multiidx2n(Cmidxcol, n, Cdims));
1348  }
1349 
1350  return sm;
1351  }; /* end worker */
1352 
1353  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1354  {
1355  // compute the column multi-indexes of the complement
1356  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1357  Cmidxcolsubsys_bar);
1358 #ifdef WITH_OPENMP_
1359 #pragma omp parallel for
1360 #endif // WITH_OPENMP_
1361  for (idx i = 0; i < Dsubsys_bar; ++i) {
1362  result(i, j) = worker(i);
1363  }
1364  }
1365 
1366  return result;
1367  }
1368  //************ Exception: not ket nor density matrix ************//
1369  else
1370  throw exception::MatrixNotSquareNorCvector("qpp::ptrace()");
1371 }
1372 
1387 template <typename Derived>
1388 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1389  const std::vector<idx>& target,
1390  idx d = 2) {
1391  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1392 
1393  // EXCEPTION CHECKS
1394 
1395  // check zero size
1397  throw exception::ZeroSize("qpp::ptrace()");
1398 
1399  // check valid dims
1400  if (d < 2)
1401  throw exception::DimsInvalid("qpp::ptrace()");
1402  // END EXCEPTION CHECKS
1403 
1404  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1405  std::vector<idx> dims(n, d); // local dimensions vector
1406 
1407  return ptrace(rA, target, dims);
1408 }
1409 
1423 template <typename Derived>
1424 dyn_mat<typename Derived::Scalar>
1425 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& target,
1426  const std::vector<idx>& dims) {
1427  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1428 
1429  // EXCEPTION CHECKS
1430 
1431  // check zero-size
1433  throw exception::ZeroSize("qpp::ptranspose()");
1434 
1435  // check that dims is a valid dimension vector
1436  if (!internal::check_dims(dims))
1437  throw exception::DimsInvalid("qpp::ptranspose()");
1438 
1439  // check that target are valid
1440  if (!internal::check_subsys_match_dims(target, dims))
1441  throw exception::SubsysMismatchDims("qpp::ptranspose()");
1442  // END EXCEPTION CHECKS
1443 
1444  idx D = static_cast<idx>(rA.rows());
1445  idx n = dims.size();
1446  idx n_subsys = target.size();
1447  idx Cdims[maxn];
1448  idx Cmidxcol[maxn];
1449  idx Csubsys[maxn];
1450 
1451  // copy dims in Cdims and target in Csubsys
1452  for (idx i = 0; i < n; ++i)
1453  Cdims[i] = dims[i];
1454  for (idx i = 0; i < n_subsys; ++i)
1455  Csubsys[i] = target[i];
1456 
1457  dyn_mat<typename Derived::Scalar> result(D, D);
1458 
1459  //************ ket ************//
1460  if (internal::check_cvector(rA)) // we have a ket
1461  {
1462  // check that dims match the dimension of A
1463  if (!internal::check_dims_match_cvect(dims, rA))
1464  throw exception::DimsMismatchCvector("qpp::ptranspose()");
1465 
1466  if (target.size() == dims.size())
1467  return (rA * adjoint(rA)).transpose();
1468 
1469  if (target.size() == 0)
1470  return rA * adjoint(rA);
1471 
1472  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1473  // use static allocation for speed!
1474  idx midxcoltmp[maxn];
1475  idx midxrow[maxn];
1476 
1477  for (idx k = 0; k < n; ++k)
1478  midxcoltmp[k] = Cmidxcol[k];
1479 
1480  /* compute the row multi-index */
1481  internal::n2multiidx(i, n, Cdims, midxrow);
1482 
1483  for (idx k = 0; k < n_subsys; ++k)
1484  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1485 
1486  /* writes the result */
1487  return rA(internal::multiidx2n(midxrow, n, Cdims)) *
1488  std::conj(rA(internal::multiidx2n(midxcoltmp, n, Cdims)));
1489  }; /* end worker */
1490 
1491  for (idx j = 0; j < D; ++j) {
1492  // compute the column multi-index
1493  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1494 
1495 #ifdef WITH_OPENMP_
1496 #pragma omp parallel for
1497 #endif // WITH_OPENMP_
1498  for (idx i = 0; i < D; ++i)
1499  result(i, j) = worker(i);
1500  }
1501 
1502  return result;
1503  }
1504  //************ density matrix ************//
1505  else if (internal::check_square_mat(rA)) // we have a density operator
1506  {
1507  // check that dims match the dimension of A
1508  if (!internal::check_dims_match_mat(dims, rA))
1509  throw exception::DimsMismatchMatrix("qpp::ptranspose()");
1510 
1511  if (target.size() == dims.size())
1512  return rA.transpose();
1513 
1514  if (target.size() == 0)
1515  return rA;
1516 
1517  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1518  // use static allocation for speed!
1519  idx midxcoltmp[maxn];
1520  idx midxrow[maxn];
1521 
1522  for (idx k = 0; k < n; ++k)
1523  midxcoltmp[k] = Cmidxcol[k];
1524 
1525  /* compute the row multi-index */
1526  internal::n2multiidx(i, n, Cdims, midxrow);
1527 
1528  for (idx k = 0; k < n_subsys; ++k)
1529  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1530 
1531  /* writes the result */
1532  return rA(internal::multiidx2n(midxrow, n, Cdims),
1533  internal::multiidx2n(midxcoltmp, n, Cdims));
1534  }; /* end worker */
1535 
1536  for (idx j = 0; j < D; ++j) {
1537  // compute the column multi-index
1538  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1539 
1540 #ifdef WITH_OPENMP_
1541 #pragma omp parallel for
1542 #endif // WITH_OPENMP_
1543  for (idx i = 0; i < D; ++i)
1544  result(i, j) = worker(i);
1545  }
1546 
1547  return result;
1548  }
1549  //************ Exception: not ket nor density matrix ************//
1550  else
1551  throw exception::MatrixNotSquareNorCvector("qpp::ptranspose()");
1552 }
1553 
1567 template <typename Derived>
1568 dyn_mat<typename Derived::Scalar>
1569 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& target,
1570  idx d = 2) {
1571  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1572 
1573  // EXCEPTION CHECKS
1574 
1575  // check zero size
1577  throw exception::ZeroSize("qpp::ptranspose()");
1578 
1579  // check valid dims
1580  if (d < 2)
1581  throw exception::DimsInvalid("qpp::ptranspose()");
1582  // END EXCEPTION CHECKS
1583 
1584  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1585  std::vector<idx> dims(n, d); // local dimensions vector
1586 
1587  return ptranspose(rA, target, dims);
1588 }
1589 
1602 template <typename Derived>
1603 dyn_mat<typename Derived::Scalar>
1604 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1605  const std::vector<idx>& dims) {
1606  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1607 
1608  // EXCEPTION CHECKS
1609 
1610  // check zero-size
1612  throw exception::ZeroSize("qpp::syspermute()");
1613 
1614  // check that dims is a valid dimension vector
1615  if (!internal::check_dims(dims))
1616  throw exception::DimsInvalid("qpp::syspermute()");
1617 
1618  // check that we have a valid permutation
1619  if (!internal::check_perm(perm))
1620  throw exception::PermInvalid("qpp::syspermute()");
1621 
1622  // check that permutation match dimensions
1623  if (perm.size() != dims.size())
1624  throw exception::PermMismatchDims("qpp::syspermute()");
1625  // END EXCEPTION CHECKS
1626 
1627  idx D = static_cast<idx>(rA.rows());
1628  idx n = dims.size();
1629 
1631 
1632  //************ ket ************//
1633  if (internal::check_cvector(rA)) // we have a column vector
1634  {
1635  idx Cdims[maxn];
1636  idx Cperm[maxn];
1637 
1638  // check that dims match the dimension of rA
1639  if (!internal::check_dims_match_cvect(dims, rA))
1640  throw exception::DimsMismatchCvector("qpp::syspermute()");
1641 
1642  // copy dims in Cdims and perm in Cperm
1643  for (idx i = 0; i < n; ++i) {
1644  Cdims[i] = dims[i];
1645  Cperm[i] = perm[i];
1646  }
1647  result.resize(D, 1);
1648 
1649  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1650  // use static allocation for speed,
1651  // double the size for matrices reshaped as vectors
1652  idx midx[maxn];
1653  idx midxtmp[maxn];
1654  idx permdims[maxn];
1655 
1656  /* compute the multi-index */
1657  internal::n2multiidx(i, n, Cdims, midx);
1658 
1659  for (idx k = 0; k < n; ++k) {
1660  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1661  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1662  }
1663 
1664  return internal::multiidx2n(midxtmp, n, permdims);
1665  }; /* end worker */
1666 
1667 #ifdef WITH_OPENMP_
1668 #pragma omp parallel for
1669 #endif // WITH_OPENMP_
1670  for (idx i = 0; i < D; ++i)
1671  result(worker(i)) = rA(i);
1672 
1673  return result;
1674  }
1675  //************ density matrix ************//
1676  else if (internal::check_square_mat(rA)) // we have a density operator
1677  {
1678  idx Cdims[2 * maxn];
1679  idx Cperm[2 * maxn];
1680 
1681  // check that dims match the dimension of rA
1682  if (!internal::check_dims_match_mat(dims, rA))
1683  throw exception::DimsMismatchMatrix("qpp::syspermute()");
1684 
1685  // copy dims in Cdims and perm in Cperm
1686  for (idx i = 0; i < n; ++i) {
1687  Cdims[i] = dims[i];
1688  Cdims[i + n] = dims[i];
1689  Cperm[i] = perm[i];
1690  Cperm[i + n] = perm[i] + n;
1691  }
1692  result.resize(D * D, 1);
1693  // map A to a column vector
1694  dyn_mat<typename Derived::Scalar> vectA =
1695  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1696  const_cast<typename Derived::Scalar*>(rA.data()), D * D, 1);
1697 
1698  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1699  // use static allocation for speed,
1700  // double the size for matrices reshaped as vectors
1701  idx midx[2 * maxn];
1702  idx midxtmp[2 * maxn];
1703  idx permdims[2 * maxn];
1704 
1705  /* compute the multi-index */
1706  internal::n2multiidx(i, 2 * n, Cdims, midx);
1707 
1708  for (idx k = 0; k < 2 * n; ++k) {
1709  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1710  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1711  }
1712 
1713  return internal::multiidx2n(midxtmp, 2 * n, permdims);
1714  }; /* end worker */
1715 
1716 #ifdef WITH_OPENMP_
1717 #pragma omp parallel for
1718 #endif // WITH_OPENMP_
1719  for (idx i = 0; i < D * D; ++i)
1720  result(worker(i)) = rA(i);
1721 
1722  return reshape(result, D, D);
1723  }
1724  //************ Exception: not ket nor density matrix ************//
1725  else
1726  throw exception::MatrixNotSquareNorCvector("qpp::syspermute()");
1727 }
1728 
1741 template <typename Derived>
1742 dyn_mat<typename Derived::Scalar>
1743 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1744  idx d = 2) {
1745  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1746 
1747  // EXCEPTION CHECKS
1748 
1749  // check zero size
1751  throw exception::ZeroSize("qpp::syspermute()");
1752 
1753  // check valid dims
1754  if (d < 2)
1755  throw exception::DimsInvalid("qpp::syspermute()");
1756  // END EXCEPTION CHECKS
1757 
1758  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1759  std::vector<idx> dims(n, d); // local dimensions vector
1760 
1761  return syspermute(rA, perm, dims);
1762 }
1763 
1764 // as in https://arxiv.org/abs/1707.08834
1776 template <typename Derived>
1777 dyn_mat<typename Derived::Scalar> applyQFT(const Eigen::MatrixBase<Derived>& A,
1778  const std::vector<idx>& target,
1779  idx d = 2, bool swap = true) {
1780  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1781 
1782  // EXCEPTION CHECKS
1783 
1784  // check zero sizes
1786  throw exception::ZeroSize("qpp::applyQFT()");
1787 
1788  // check valid subsystem dimension
1789  if (d < 2)
1790  throw exception::DimsInvalid("qpp::applyQFT()");
1791 
1792  // total number of qubits/qudits in the state
1793  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1794 
1795  std::vector<idx> dims(n, d); // local dimensions vector
1796  // check that target is valid w.r.t. dims
1797  if (!internal::check_subsys_match_dims(target, dims))
1798  throw exception::SubsysMismatchDims("qpp::applyQFT()");
1799 
1800  //************ ket ************//
1801  if (internal::check_cvector(rA)) // we have a ket
1802  {
1803  // check that dims match state vector
1804  if (!internal::check_dims_match_cvect(dims, rA))
1805  throw exception::DimsMismatchCvector("qpp::applyQFT()");
1806  }
1807  //************ density matrix ************//
1808  else if (internal::check_square_mat(rA)) // we have a density operator
1809  {
1810  // check that dims match state matrix
1811  if (!internal::check_dims_match_mat(dims, rA))
1812  throw exception::DimsMismatchMatrix("qpp::applyQFT()");
1813  }
1814  //************ Exception: not ket nor density matrix ************//
1815  else
1816  throw exception::MatrixNotSquareNorCvector("qpp::applyQFT()");
1817  // END EXCEPTION CHECKS
1818 
1820 
1821  idx n_subsys = target.size();
1822 
1823  if (d == 2) // qubits
1824  {
1825  for (idx i = 0; i < n_subsys; ++i) {
1826  // apply Hadamard on qubit i
1827  result = apply(result, Gates::get_instance().H, {target[i]});
1828  // apply controlled rotations
1829  for (idx j = 2; j <= n_subsys - i; ++j) {
1830  // construct Rj
1831  cmat Rj(2, 2);
1832  Rj << 1, 0, 0, exp(2.0 * pi * 1_i / std::pow(2, j));
1833  result =
1834  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1835  }
1836  }
1837  if (swap) {
1838  // we have the qubits in reversed order, we must swap them
1839  for (idx i = 0; i < n_subsys / 2; ++i) {
1840  result = apply(result, Gates::get_instance().SWAP,
1841  {target[i], target[n_subsys - i - 1]});
1842  }
1843  }
1844 
1845  } else { // qudits
1846  for (idx i = 0; i < n_subsys; ++i) {
1847  // apply qudit Fourier on qudit i
1848  result = apply(result, Gates::get_instance().Fd(d), {target[i]}, d);
1849  // apply controlled rotations
1850  for (idx j = 2; j <= n_subsys - i; ++j) {
1851  // construct Rj
1852  cmat Rj = cmat::Zero(d, d);
1853  for (idx m = 0; m < d; ++m) {
1854  Rj(m, m) = exp(2.0 * pi * m * 1_i / std::pow(d, j));
1855  }
1856  result =
1857  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1858  }
1859  }
1860  if (swap) {
1861  // we have the qudits in reversed order, we must swap them
1862  for (idx i = 0; i < n_subsys / 2; ++i) {
1863  result = apply(result, Gates::get_instance().SWAPd(d),
1864  {target[i], target[n_subsys - i - 1]}, d);
1865  }
1866  }
1867  }
1868 
1869  return result;
1870 }
1871 
1872 // as in https://arxiv.org/abs/1707.08834
1884 template <typename Derived>
1885 dyn_mat<typename Derived::Scalar> applyTFQ(const Eigen::MatrixBase<Derived>& A,
1886  const std::vector<idx>& target,
1887  idx d = 2, bool swap = true) {
1888  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1889 
1890  // EXCEPTION CHECKS
1891 
1892  // check zero sizes
1894  throw exception::ZeroSize("qpp::applyTFQ()");
1895 
1896  // check valid subsystem dimension
1897  if (d < 2)
1898  throw exception::DimsInvalid("qpp::applyTFQ()");
1899 
1900  // total number of qubits/qudits in the state
1901  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1902 
1903  std::vector<idx> dims(n, d); // local dimensions vector
1904  // check that target is valid w.r.t. dims
1905  if (!internal::check_subsys_match_dims(target, dims))
1906  throw exception::SubsysMismatchDims("qpp::applyTFQ()");
1907 
1908  //************ ket ************//
1909  if (internal::check_cvector(rA)) // we have a ket
1910  {
1911  // check that dims match state vector
1912  if (!internal::check_dims_match_cvect(dims, rA))
1913  throw exception::DimsMismatchCvector("qpp::applyTFQ()");
1914  }
1915  //************ density matrix ************//
1916  else if (internal::check_square_mat(rA)) // we have a density operator
1917  {
1918  // check that dims match state matrix
1919  if (!internal::check_dims_match_mat(dims, rA))
1920  throw exception::DimsMismatchMatrix("qpp::applyTFQ()");
1921  }
1922  //************ Exception: not ket nor density matrix ************//
1923  else
1924  throw exception::MatrixNotSquareNorCvector("qpp::applyTFQ()");
1925  // END EXCEPTION CHECKS
1926 
1928 
1929  idx n_subsys = target.size();
1930 
1931  if (d == 2) // qubits
1932  {
1933  if (swap) {
1934  // we have the qubits in reversed order, we must swap them
1935  for (idx i = 0; i < n_subsys / 2; ++i) {
1936  result = apply(result, Gates::get_instance().SWAP,
1937  {target[i], target[n_subsys - i - 1]});
1938  }
1939  }
1940  for (idx i = n_subsys; i-- > 0;) {
1941  // apply controlled rotations
1942  for (idx j = n_subsys - i + 1; j-- > 2;) {
1943  // construct Rj
1944  cmat Rj(2, 2);
1945  Rj << 1, 0, 0, exp(-2.0 * pi * 1_i / std::pow(2, j));
1946  result =
1947  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]});
1948  }
1949  // apply Hadamard on qubit i
1950  result = apply(result, Gates::get_instance().H, {target[i]});
1951  }
1952  } else { // qudits
1953  if (swap) {
1954  // we have the qudits in reversed order, we must swap them
1955  for (idx i = 0; i < n_subsys / 2; ++i) {
1956  result = apply(result, Gates::get_instance().SWAPd(d),
1957  {target[i], target[n_subsys - i - 1]}, d);
1958  }
1959  }
1960  for (idx i = n_subsys; i-- > 0;) {
1961  // apply controlled rotations
1962  for (idx j = n_subsys - i + 1; j-- > 2;) {
1963  // construct Rj
1964  cmat Rj = cmat::Zero(d, d);
1965  for (idx m = 0; m < d; ++m) {
1966  Rj(m, m) = exp(-2.0 * pi * m * 1_i / std::pow(d, j));
1967  }
1968  result =
1969  applyCTRL(result, Rj, {target[i + j - 1]}, {target[i]}, d);
1970  }
1971  // apply qudit Fourier on qudit i
1972  result = apply(result, adjoint(Gates::get_instance().Fd(d)),
1973  {target[i]}, d);
1974  }
1975  }
1976 
1977  return result;
1978 }
1979 
1980 // as in https://arxiv.org/abs/1707.08834
1989 template <typename Derived>
1990 dyn_col_vect<typename Derived::Scalar> TFQ(const Eigen::MatrixBase<Derived>& A,
1991  idx d = 2, bool swap = true) {
1992  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1993 
1994  // EXCEPTION CHECKS
1995 
1996  // check zero-size
1998  throw exception::ZeroSize("qpp::TFQ()");
1999 
2000  // check valid subsystem dimension
2001  if (d < 2)
2002  throw exception::DimsInvalid("qpp::TFQ()");
2003 
2004  // total number of qubits/qudits in the state
2005  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
2006 
2007  std::vector<idx> dims(n, d); // local dimensions vector
2008 
2009  //************ ket ************//
2010  if (internal::check_cvector(rA)) // we have a ket
2011  {
2012  // check that dims match state vector
2013  if (!internal::check_dims_match_cvect(dims, rA))
2014  throw exception::DimsMismatchCvector("qpp::TFQ()");
2015  }
2016  //************ density matrix ************//
2017  else if (internal::check_square_mat(rA)) // we have a density operator
2018  {
2019  // check that dims match state matrix
2020  if (!internal::check_dims_match_mat(dims, rA))
2021  throw exception::DimsMismatchMatrix("qpp::TFQ()");
2022  }
2023  //************ Exception: not ket nor density matrix ************//
2024  else
2025  throw exception::MatrixNotSquareNorCvector("qpp::TFQ()");
2026  // END EXCEPTION CHECKS
2027 
2028  std::vector<idx> subsys(n);
2029  std::iota(std::begin(subsys), std::end(subsys), 0);
2030  ket result = applyTFQ(rA, subsys, d, swap);
2031 
2032  return result;
2033 }
2034 
2035 // as in https://arxiv.org/abs/1707.08834
2044 template <typename Derived>
2045 dyn_col_vect<typename Derived::Scalar> QFT(const Eigen::MatrixBase<Derived>& A,
2046  idx d = 2, bool swap = true) {
2047  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
2048 
2049  // EXCEPTION CHECKS
2050 
2051  // check zero-size
2053  throw exception::ZeroSize("qpp::QFT()");
2054 
2055  // check valid subsystem dimension
2056  if (d < 2)
2057  throw exception::DimsInvalid("qpp::QFT()");
2058 
2059  // total number of qubits/qudits in the state
2060  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
2061 
2062  std::vector<idx> dims(n, d); // local dimensions vector
2063 
2064  //************ ket ************//
2065  if (internal::check_cvector(rA)) // we have a ket
2066  {
2067  // check that dims match state vector
2068  if (!internal::check_dims_match_cvect(dims, rA))
2069  throw exception::DimsMismatchCvector("qpp::QFT()");
2070  }
2071  //************ density matrix ************//
2072  else if (internal::check_square_mat(rA)) // we have a density operator
2073  {
2074  // check that dims match state matrix
2075  if (!internal::check_dims_match_mat(dims, rA))
2076  throw exception::DimsMismatchMatrix("qpp::QFT()");
2077  }
2078  //************ Exception: not ket nor density matrix ************//
2079  else
2080  throw exception::MatrixNotSquareNorCvector("qpp::QFT()");
2081  // END EXCEPTION CHECKS
2082 
2083  std::vector<idx> subsys(n);
2084  std::iota(std::begin(subsys), std::end(subsys), 0);
2085  ket result = applyQFT(rA, subsys, d, swap);
2086 
2087  return result;
2088 }
2089 
2090 } /* namespace qpp */
2091 
2092 #endif /* OPERATIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimensions not equal exception.
Definition: exception.h:284
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:67
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:219
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:830
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
dyn_mat< typename Derived::Scalar > applyTFQ(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2, bool swap=true)
Applies the inverse (adjoint) qudit quantum Fourier transform to the part target of the multi-partite...
Definition: operations.h:1885
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
dyn_mat< typename Derived::Scalar > syspermute(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &perm, idx d=2)
Subsystem permutation.
Definition: operations.h:1743
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: circuits.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:263
Matrix is not square nor column vector exception.
Definition: exception.h:209
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:425
Invalid dimension(s) exception.
Definition: exception.h:269
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:863
idx get_num_subsys(idx D, idx d)
Definition: util.h:370
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:897
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
Not bi-partite exception.
Definition: exception.h:486
dyn_mat< typename Derived1::Scalar > apply(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &target, const std::vector< idx > &dims)
Applies the gate A to the part target of the multi-partite state vector or density matrix state...
Definition: operations.h:444
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:938
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1396
Invalid permutation exception.
Definition: exception.h:380
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:89
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &target, const std::vector< idx > &dims)
Applies the controlled-gate A to the part target of the multi-partite state vector or density matrix ...
Definition: operations.h:55
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1424
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:917
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &target, idx d=2)
Applies the controlled-gate A to the part target of the multi-partite state vector or density matrix ...
Definition: operations.h:405
Matrix mismatch subsystems exception.
Definition: exception.h:254
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, idx d=2)
Partial trace.
Definition: operations.h:1152
Dimension(s) mismatch column vector size exception.
Definition: exception.h:316
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:46
dyn_col_vect< typename Derived::Scalar > TFQ(const Eigen::MatrixBase< Derived > &A, idx d=2, bool swap=true)
Inverse (adjoint) qudit quantum Fourier transform.
Definition: operations.h:1990
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Fast matrix power based on the SQUARE-AND-MULTIPLY algorithm.
Definition: functions.h:800
Type mismatch exception.
Definition: exception.h:530
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:382
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors of Hermitian matrix.
Definition: functions.h:450
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:775
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2)
Partial trace.
Definition: operations.h:1388
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:130
constexpr double pi
Definition: constants.h:72
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:704
dyn_mat< typename Derived::Scalar > applyQFT(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2, bool swap=true)
Applies the qudit quantum Fourier transform to the part target of the multi-partite state vector or d...
Definition: operations.h:1777
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
cmat apply(const Eigen::MatrixBase< Derived > &A, const std::vector< cmat > &Ks, const std::vector< idx > &target, idx d=2)
Applies the channel specified by the set of Kraus operators Ks to the part target of the multi-partit...
Definition: operations.h:672
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
Permutation mismatch dimensions exception.
Definition: exception.h:396
dyn_mat< typename Derived::Scalar > ptranspose(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &target, idx d=2)
Partial transpose.
Definition: operations.h:1569
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, idx d=2)
Partial trace.
Definition: operations.h:1032
std::size_t idx
Non-negative integer index, make sure you use an unsigned type.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:149
std::vector< idx > complement(std::vector< idx > subsys, idx n)
Constructs the complement of a subsystem vector.
Definition: functions.h:1780
dyn_col_vect< typename Derived::Scalar > QFT(const Eigen::MatrixBase< Derived > &A, idx d=2, bool swap=true)
Qudit quantum Fourier transform.
Definition: operations.h:2045
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:207
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1144
Object has zero size exception.
Definition: exception.h:134