Quantum++  v1.1
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>& subsys, 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 square matrix for the gate
79  throw exception::MatrixNotSquare("qpp::applyCTRL()");
80 
81  // check that all control subsystems have the same dimension
82  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
83  for (idx i = 1; i < ctrl.size(); ++i)
84  if (dims[ctrl[i]] != d)
85  throw exception::DimsNotEqual("qpp::applyCTRL()");
86 
87  // check that dimension is valid
88  if (!internal::check_dims(dims))
89  throw exception::DimsInvalid("qpp::applyCTRL()");
90 
91  // check subsys is valid w.r.t. dims
92  if (!internal::check_subsys_match_dims(subsys, dims))
93  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
94 
95  // check that gate matches the dimensions of the subsys
96  std::vector<idx> subsys_dims(subsys.size());
97  for (idx i = 0; i < subsys.size(); ++i)
98  subsys_dims[i] = dims[subsys[i]];
99  if (!internal::check_dims_match_mat(subsys_dims, rA))
100  throw exception::MatrixMismatchSubsys("qpp::applyCTRL()");
101 
102  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
103  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
104  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
105 
106  // check that ctrl + gate subsystem is valid
107  // with respect to local dimensions
108  if (!internal::check_subsys_match_dims(ctrlgate, dims))
109  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
110  // END EXCEPTION CHECKS
111 
112  // construct the table of A^i and (A^dagger)^i
113  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
114  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
115  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i) {
116  Ai.push_back(powm(rA, i));
117  Aidagger.push_back(powm(adjoint(rA), i));
118  }
119 
120  idx D = static_cast<idx>(rstate.rows()); // total dimension
121  idx n = dims.size(); // total number of subsystems
122  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
123  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
124  idx subsyssize = subsys.size(); // number of subsystems of the target
125  // dimension of ctrl subsystem
126  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
127  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
128 
129  idx Cdims[maxn]; // local dimensions
130  idx CdimsA[maxn]; // local dimensions
131  idx CdimsCTRL[maxn]; // local dimensions
132  idx CdimsCTRLA_bar[maxn]; // local dimensions
133 
134  // compute the complementary subsystem of ctrlgate w.r.t. dims
135  std::vector<idx> ctrlgate_bar = complement(ctrlgate, n);
136  // number of subsystems that are complementary to the ctrl+gate
137  idx ctrlgate_barsize = ctrlgate_bar.size();
138 
139  idx DCTRLA_bar = 1; // dimension of the rest
140  for (idx i = 0; i < ctrlgate_barsize; ++i)
141  DCTRLA_bar *= dims[ctrlgate_bar[i]];
142 
143  for (idx k = 0; k < n; ++k)
144  Cdims[k] = dims[k];
145  for (idx k = 0; k < subsyssize; ++k)
146  CdimsA[k] = dims[subsys[k]];
147  for (idx k = 0; k < ctrlsize; ++k)
148  CdimsCTRL[k] = d;
149  for (idx k = 0; k < ctrlgate_barsize; ++k)
150  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
151 
152  // worker, computes the coefficient and the index for the ket case
153  // used in #pragma omp parallel for collapse
154  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
155  ->std::pair<typename Derived1::Scalar, idx> {
156  idx indx = 0;
157  typename Derived1::Scalar coeff = 0;
158 
159  idx Cmidx[maxn]; // the total multi-index
160  idx CmidxA[maxn]; // the gate part multi-index
161  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
162 
163  // compute the index
164 
165  // set the CTRL part
166  for (idx k = 0; k < ctrlsize; ++k) {
167  Cmidx[ctrl[k]] = i_;
168  }
169 
170  // set the rest
171  internal::n2multiidx(r_, n - ctrlgatesize, CdimsCTRLA_bar,
172  CmidxCTRLA_bar);
173  for (idx k = 0; k < n - ctrlgatesize; ++k) {
174  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
175  }
176 
177  // set the A part
178  internal::n2multiidx(m_, subsyssize, CdimsA, CmidxA);
179  for (idx k = 0; k < subsyssize; ++k) {
180  Cmidx[subsys[k]] = CmidxA[k];
181  }
182 
183  // we now got the total index
184  indx = internal::multiidx2n(Cmidx, n, Cdims);
185 
186  // compute the coefficient
187  for (idx n_ = 0; n_ < DA; ++n_) {
188  internal::n2multiidx(n_, subsyssize, CdimsA, CmidxA);
189  for (idx k = 0; k < subsyssize; ++k) {
190  Cmidx[subsys[k]] = CmidxA[k];
191  }
192  coeff +=
193  Ai[i_](m_, n_) * rstate(internal::multiidx2n(Cmidx, n, Cdims));
194  }
195 
196  return std::make_pair(coeff, indx);
197  }; /* end coeff_idx_ket */
198 
199  // worker, computes the coefficient and the index
200  // for the density matrix case
201  // used in #pragma omp parallel for collapse
202  auto coeff_idx_rho = [&](idx i1_, idx m1_, idx r1_, idx i2_, idx m2_,
203  idx r2_) noexcept
204  ->std::tuple<typename Derived1::Scalar, idx, idx> {
205  idx idxrow = 0;
206  idx idxcol = 0;
207  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
208 
209  idx Cmidxrow[maxn]; // the total row multi-index
210  idx Cmidxcol[maxn]; // the total col multi-index
211  idx CmidxArow[maxn]; // the gate part row multi-index
212  idx CmidxAcol[maxn]; // the gate part col multi-index
213  idx CmidxCTRLrow[maxn]; // the control row multi-index
214  idx CmidxCTRLcol[maxn]; // the control col multi-index
215  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
216  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
217 
218  // compute the ket/bra indexes
219 
220  // set the CTRL part
221  internal::n2multiidx(i1_, ctrlsize, CdimsCTRL, CmidxCTRLrow);
222  internal::n2multiidx(i2_, ctrlsize, CdimsCTRL, CmidxCTRLcol);
223 
224  for (idx k = 0; k < ctrlsize; ++k) {
225  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
226  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
227  }
228 
229  // set the rest
230  internal::n2multiidx(r1_, n - ctrlgatesize, CdimsCTRLA_bar,
231  CmidxCTRLA_barrow);
232  internal::n2multiidx(r2_, n - ctrlgatesize, CdimsCTRLA_bar,
233  CmidxCTRLA_barcol);
234  for (idx k = 0; k < n - ctrlgatesize; ++k) {
235  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
236  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
237  }
238 
239  // set the A part
240  internal::n2multiidx(m1_, subsyssize, CdimsA, CmidxArow);
241  internal::n2multiidx(m2_, subsyssize, CdimsA, CmidxAcol);
242  for (idx k = 0; k < subsys.size(); ++k) {
243  Cmidxrow[subsys[k]] = CmidxArow[k];
244  Cmidxcol[subsys[k]] = CmidxAcol[k];
245  }
246 
247  // we now got the total row/col indexes
248  idxrow = internal::multiidx2n(Cmidxrow, n, Cdims);
249  idxcol = internal::multiidx2n(Cmidxcol, n, Cdims);
250 
251  // check whether all CTRL row and col multi indexes are equal
252  bool all_ctrl_rows_equal = true;
253  bool all_ctrl_cols_equal = true;
254 
255  idx first_ctrl_row, first_ctrl_col;
256  if (ctrlsize > 0) {
257  first_ctrl_row = CmidxCTRLrow[0];
258  first_ctrl_col = CmidxCTRLcol[0];
259  } else {
260  first_ctrl_row = first_ctrl_col = 1;
261  }
262 
263  for (idx k = 1; k < ctrlsize; ++k) {
264  if (CmidxCTRLrow[k] != first_ctrl_row) {
265  all_ctrl_rows_equal = false;
266  break;
267  }
268  }
269  for (idx k = 1; k < ctrlsize; ++k) {
270  if (CmidxCTRLcol[k] != first_ctrl_col) {
271  all_ctrl_cols_equal = false;
272  break;
273  }
274  }
275 
276  // at least one control activated, compute the coefficient
277  for (idx n1_ = 0; n1_ < DA; ++n1_) {
278  internal::n2multiidx(n1_, subsyssize, CdimsA, CmidxArow);
279  for (idx k = 0; k < subsyssize; ++k) {
280  Cmidxrow[subsys[k]] = CmidxArow[k];
281  }
282  idx idxrowtmp = internal::multiidx2n(Cmidxrow, n, Cdims);
283 
284  if (all_ctrl_rows_equal) {
285  lhs = Ai[first_ctrl_row](m1_, n1_);
286  } else {
287  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
288  }
289 
290  for (idx n2_ = 0; n2_ < DA; ++n2_) {
291  internal::n2multiidx(n2_, subsyssize, CdimsA, CmidxAcol);
292  for (idx k = 0; k < subsyssize; ++k) {
293  Cmidxcol[subsys[k]] = CmidxAcol[k];
294  }
295 
296  if (all_ctrl_cols_equal) {
297  rhs = Aidagger[first_ctrl_col](n2_, m2_);
298  } else {
299  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
300  }
301 
302  idx idxcoltmp = internal::multiidx2n(Cmidxcol, n, Cdims);
303 
304  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
305  }
306  }
307 
308  return std::make_tuple(coeff, idxrow, idxcol);
309  }; /* end coeff_idx_rho */
310 
311  //************ ket ************//
312  if (internal::check_cvector(rstate)) // we have a ket
313  {
314  // check that dims match state vector
315  if (!internal::check_dims_match_cvect(dims, rstate))
316  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
317  if (D == 1)
318  return rstate;
319 
320  dyn_mat<typename Derived1::Scalar> result = rstate;
321 
322 #ifdef WITH_OPENMP_
323 #pragma omp parallel for collapse(2)
324 #endif // WITH_OPENMP_
325  for (idx m = 0; m < DA; ++m)
326  for (idx r = 0; r < DCTRLA_bar; ++r) {
327  if (ctrlsize == 0) // no control
328  {
329  result(coeff_idx_ket(1, m, r).second) =
330  coeff_idx_ket(1, m, r).first;
331  } else
332  for (idx i = 0; i < d; ++i) {
333  result(coeff_idx_ket(i, m, r).second) =
334  coeff_idx_ket(i, m, r).first;
335  }
336  }
337 
338  return result;
339  }
340  //************ density matrix ************//
341  else if (internal::check_square_mat(rstate)) // we have a density operator
342  {
343  // check that dims match state matrix
344  if (!internal::check_dims_match_mat(dims, rstate))
345  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
346 
347  if (D == 1)
348  return rstate;
349 
350  dyn_mat<typename Derived1::Scalar> result = rstate;
351 
352 #ifdef WITH_OPENMP_
353 #pragma omp parallel for collapse(4)
354 #endif // WITH_OPENMP_
355  for (idx m1 = 0; m1 < DA; ++m1)
356  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
357  for (idx m2 = 0; m2 < DA; ++m2)
358  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
359  if (ctrlsize == 0) // no control
360  {
361  auto coeff_idxes =
362  coeff_idx_rho(1, m1, r1, 1, m2, r2);
363  result(std::get<1>(coeff_idxes),
364  std::get<2>(coeff_idxes)) =
365  std::get<0>(coeff_idxes);
366  } else {
367  for (idx i1 = 0; i1 < Dctrl; ++i1)
368  for (idx i2 = 0; i2 < Dctrl; ++i2) {
369  auto coeff_idxes =
370  coeff_idx_rho(i1, m1, r1, i2, m2, r2);
371  result(std::get<1>(coeff_idxes),
372  std::get<2>(coeff_idxes)) =
373  std::get<0>(coeff_idxes);
374  }
375  }
376 
377  return result;
378  }
379  //************ Exception: not ket nor density matrix ************//
380  else
381  throw exception::MatrixNotSquareNorCvector("qpp::applyCTRL()");
382 }
383 
399 template <typename Derived1, typename Derived2>
401 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
402  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
403  const std::vector<idx>& subsys, idx d = 2) {
404  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
405  state.derived();
406  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
407 
408  // EXCEPTION CHECKS
409 
410  // check zero size
411  if (!internal::check_nonzero_size(rstate))
412  throw exception::ZeroSize("qpp::applyCTRL()");
413 
414  // check valid dims
415  if (d < 2)
416  throw exception::DimsInvalid("qpp::applyCTRL()");
417  // END EXCEPTION CHECKS
418 
419  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
420  std::vector<idx> dims(n, d); // local dimensions vector
421 
422  return applyCTRL(rstate, rA, ctrl, subsys, dims);
423 }
424 
438 template <typename Derived1, typename Derived2>
440 apply(const Eigen::MatrixBase<Derived1>& state,
441  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& subsys,
442  const std::vector<idx>& dims) {
443  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
444  state.derived();
445  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
446 
447  // EXCEPTION CHECKS
448 
449  // check types
450  if (!std::is_same<typename Derived1::Scalar,
451  typename Derived2::Scalar>::value)
452  throw exception::TypeMismatch("qpp::apply()");
453 
454  // check zero sizes
456  throw exception::ZeroSize("qpp::apply()");
457 
458  // check zero sizes
459  if (!internal::check_nonzero_size(rstate))
460  throw exception::ZeroSize("qpp::apply()");
461 
462  // check square matrix for the gate
464  throw exception::MatrixNotSquare("qpp::apply()");
465 
466  // check that dimension is valid
467  if (!internal::check_dims(dims))
468  throw exception::DimsInvalid("qpp::apply()");
469 
470  // check subsys is valid w.r.t. dims
471  if (!internal::check_subsys_match_dims(subsys, dims))
472  throw exception::SubsysMismatchDims("qpp::apply()");
473 
474  // check that gate matches the dimensions of the subsys
475  std::vector<idx> subsys_dims(subsys.size());
476  for (idx i = 0; i < subsys.size(); ++i)
477  subsys_dims[i] = dims[subsys[i]];
478  if (!internal::check_dims_match_mat(subsys_dims, rA))
479  throw exception::MatrixMismatchSubsys("qpp::apply()");
480  // END EXCEPTION CHECKS
481 
482  //************ ket ************//
483  if (internal::check_cvector(rstate)) // we have a ket
484  {
485  // check that dims match state vector
486  if (!internal::check_dims_match_cvect(dims, rstate))
487  throw exception::DimsMismatchCvector("qpp::apply()");
488 
489  return applyCTRL(rstate, rA, {}, subsys, dims);
490  }
491  //************ density matrix ************//
492  else if (internal::check_square_mat(rstate)) // we have a density operator
493  {
494  // check that dims match state matrix
495  if (!internal::check_dims_match_mat(dims, rstate))
496  throw exception::DimsMismatchMatrix("qpp::apply()");
497 
498  return applyCTRL(rstate, rA, {}, subsys, dims);
499  }
500  //************ Exception: not ket nor density matrix ************//
501  else
502  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
503 }
504 
518 template <typename Derived1, typename Derived2>
520 apply(const Eigen::MatrixBase<Derived1>& state,
521  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& subsys,
522  idx d = 2) {
523  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
524  state.derived();
525  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
526 
527  // EXCEPTION CHECKS
528 
529  // check zero size
530  if (!internal::check_nonzero_size(rstate))
531  throw exception::ZeroSize("qpp::apply()");
532 
533  // check valid dims
534  if (d < 2)
535  throw exception::DimsInvalid("qpp::apply()");
536  // END EXCEPTION CHECKS
537 
538  idx n = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
539  std::vector<idx> dims(n, d); // local dimensions vector
540 
541  return apply(rstate, rA, subsys, dims);
542 }
543 
552 template <typename Derived>
553 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
554  const cmat& rA = A.derived();
555 
556  // EXCEPTION CHECKS
557 
559  throw exception::ZeroSize("qpp::apply()");
561  throw exception::MatrixNotSquare("qpp::apply()");
562  if (Ks.size() == 0)
563  throw exception::ZeroSize("qpp::apply()");
564  if (!internal::check_square_mat(Ks[0]))
565  throw exception::MatrixNotSquare("qpp::apply()");
566  if (Ks[0].rows() != rA.rows())
567  throw exception::DimsMismatchMatrix("qpp::apply()");
568  for (auto&& it : Ks)
569  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
570  throw exception::DimsNotEqual("qpp::apply()");
571  // END EXCEPTION CHECKS
572 
573  cmat result = cmat::Zero(rA.rows(), rA.rows());
574 
575 #ifdef WITH_OPENMP_
576 #pragma omp parallel for
577 #endif // WITH_OPENMP_
578  for (idx i = 0; i < Ks.size(); ++i) {
579 #ifdef WITH_OPENMP_
580 #pragma omp critical
581 #endif // WITH_OPENMP_
582  { result += Ks[i] * rA * adjoint(Ks[i]); }
583  }
584 
585  return result;
586 }
587 
598 template <typename Derived>
599 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
600  const std::vector<idx>& subsys, const std::vector<idx>& dims) {
601  const cmat& rA = A.derived();
602 
603  // EXCEPTION CHECKS
604 
605  // check zero sizes
607  throw exception::ZeroSize("qpp::apply()");
608 
609  // check square matrix for the A
611  throw exception::MatrixNotSquare("qpp::apply()");
612 
613  // check that dimension is valid
614  if (!internal::check_dims(dims))
615  throw exception::DimsInvalid("qpp::apply()");
616 
617  // check that dims match A matrix
618  if (!internal::check_dims_match_mat(dims, rA))
619  throw exception::DimsMismatchMatrix("qpp::apply()");
620 
621  // check subsys is valid w.r.t. dims
622  if (!internal::check_subsys_match_dims(subsys, dims))
623  throw exception::SubsysMismatchDims("qpp::apply()");
624 
625  std::vector<idx> subsys_dims(subsys.size());
626  for (idx i = 0; i < subsys.size(); ++i)
627  subsys_dims[i] = dims[subsys[i]];
628 
629  // check the Kraus operators
630  if (Ks.size() == 0)
631  throw exception::ZeroSize("qpp::apply()");
632  if (!internal::check_square_mat(Ks[0]))
633  throw exception::MatrixNotSquare("qpp::apply()");
634  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
635  throw exception::MatrixMismatchSubsys("qpp::apply()");
636  for (auto&& it : Ks)
637  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
638  throw exception::DimsNotEqual("qpp::apply()");
639  // END EXCEPTION CHECKS
640 
641  cmat result = cmat::Zero(rA.rows(), rA.rows());
642 
643  for (idx i = 0; i < Ks.size(); ++i)
644  result += apply(rA, Ks[i], subsys, dims);
645 
646  return result;
647 }
648 
659 template <typename Derived>
660 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
661  const std::vector<idx>& subsys, idx d = 2) {
662  const cmat& rA = A.derived();
663 
664  // EXCEPTION CHECKS
665 
666  // check zero sizes
668  throw exception::ZeroSize("qpp::apply()");
669 
670  // check valid dims
671  if (d < 2)
672  throw exception::DimsInvalid("qpp::apply()");
673  // END EXCEPTION CHECKS
674 
675  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
676  std::vector<idx> dims(n, d); // local dimensions vector
677 
678  return apply(rA, Ks, subsys, dims);
679 }
680 
692 inline cmat kraus2super(const std::vector<cmat>& Ks) {
693  // EXCEPTION CHECKS
694 
695  if (Ks.size() == 0)
696  throw exception::ZeroSize("qpp::kraus2super()");
697  if (!internal::check_nonzero_size(Ks[0]))
698  throw exception::ZeroSize("qpp::kraus2super()");
699  if (!internal::check_square_mat(Ks[0]))
700  throw exception::MatrixNotSquare("qpp::kraus2super()");
701  for (auto&& it : Ks)
702  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
703  throw exception::DimsNotEqual("qpp::kraus2super()");
704  // END EXCEPTION CHECKS
705 
706  idx D = static_cast<idx>(Ks[0].rows());
707 
708  cmat result(D * D, D * D);
709  cmat MN = cmat::Zero(D, D);
710  bra A = bra::Zero(D);
711  ket B = ket::Zero(D);
712  cmat EMN = cmat::Zero(D, D);
713 
714 #ifdef WITH_OPENMP_
715 #pragma omp parallel for collapse(2)
716 #endif // WITH_OPENMP_
717  for (idx m = 0; m < D; ++m) {
718  for (idx n = 0; n < D; ++n) {
719 #ifdef WITH_OPENMP_
720 #pragma omp critical
721 #endif // WITH_OPENMP_
722  {
723  // compute E(|m><n|)
724  MN(m, n) = 1;
725  for (idx i = 0; i < Ks.size(); ++i)
726  EMN += Ks[i] * MN * adjoint(Ks[i]);
727  MN(m, n) = 0;
728 
729  for (idx a = 0; a < D; ++a) {
730  A(a) = 1;
731  for (idx b = 0; b < D; ++b) {
732  // compute result(ab,mn)=<a|E(|m><n)|b>
733  B(b) = 1;
734  result(a * D + b, m * D + n) =
735  static_cast<cmat>(A * EMN * B).value();
736  B(b) = 0;
737  }
738  A(a) = 0;
739  }
740  EMN = cmat::Zero(D, D);
741  }
742  }
743  }
744 
745  return result;
746 }
747 
763 inline cmat kraus2choi(const std::vector<cmat>& Ks) {
764  // EXCEPTION CHECKS
765 
766  if (Ks.size() == 0)
767  throw exception::ZeroSize("qpp::kraus2choi()");
768  if (!internal::check_nonzero_size(Ks[0]))
769  throw exception::ZeroSize("qpp::kraus2choi()");
770  if (!internal::check_square_mat(Ks[0]))
771  throw exception::MatrixNotSquare("qpp::kraus2choi()");
772  for (auto&& it : Ks)
773  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
774  throw exception::DimsNotEqual("qpp::kraus2choi()");
775  // END EXCEPTION CHECKS
776 
777  idx D = static_cast<idx>(Ks[0].rows());
778 
779  // construct the D x D \sum |jj> vector
780  // (un-normalized maximally entangled state)
781  cmat MES = cmat::Zero(D * D, 1);
782  for (idx a = 0; a < D; ++a)
783  MES(a * D + a) = 1;
784 
785  cmat Omega = MES * adjoint(MES);
786 
787  cmat result = cmat::Zero(D * D, D * D);
788 
789 #ifdef WITH_OPENMP_
790 #pragma omp parallel for
791 #endif // WITH_OPENMP_
792  for (idx i = 0; i < Ks.size(); ++i) {
793 #ifdef WITH_OPENMP_
794 #pragma omp critical
795 #endif // WITH_OPENMP_
796  {
797  result += kron(cmat::Identity(D, D), Ks[i]) * Omega *
798  adjoint(kron(cmat::Identity(D, D), Ks[i]));
799  }
800  }
801 
802  return result;
803 }
804 
818 inline std::vector<cmat> choi2kraus(const cmat& A) {
819  // EXCEPTION CHECKS
820 
822  throw exception::ZeroSize("qpp::choi2kraus()");
824  throw exception::MatrixNotSquare("qpp::choi2kraus()");
825  idx D = internal::get_dim_subsys(A.rows(), 2);
826  // check equal dimensions
827  if (D * D != static_cast<idx>(A.rows()))
828  throw exception::DimsInvalid("qpp::choi2kraus()");
829  // END EXCEPTION CHECKS
830 
831  dmat ev = hevals(A);
832  cmat evec = hevects(A);
833  std::vector<cmat> result;
834 
835  for (idx i = 0; i < D * D; ++i) {
836  if (std::abs(ev(i)) > eps)
837  result.push_back(std::sqrt(std::abs(ev(i))) *
838  reshape(evec.col(i), D, D));
839  }
840 
841  return result;
842 }
843 
851 inline cmat choi2super(const cmat& A) {
852  // EXCEPTION CHECKS
853 
855  throw exception::ZeroSize("qpp::choi2super()");
857  throw exception::MatrixNotSquare("qpp::choi2super()");
858  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
859  // check equal dimensions
860  if (D * D != static_cast<idx>(A.rows()))
861  throw exception::DimsInvalid("qpp::choi2super()");
862  // END EXCEPTION CHECKS
863 
864  cmat result(D * D, D * D);
865 
866 #ifdef WITH_OPENMP_
867 #pragma omp parallel for collapse(4)
868 #endif // WITH_OPENMP_
869  for (idx a = 0; a < D; ++a)
870  for (idx b = 0; b < D; ++b)
871  for (idx m = 0; m < D; ++m)
872  for (idx n = 0; n < D; ++n)
873  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
874 
875  return result;
876 }
877 
885 inline cmat super2choi(const cmat& A) {
886  // EXCEPTION CHECKS
887 
889  throw exception::ZeroSize("qpp::super2choi()");
891  throw exception::MatrixNotSquare("qpp::super2choi()");
892  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
893  // check equal dimensions
894  if (D * D != static_cast<idx>(A.rows()))
895  throw exception::DimsInvalid("qpp::super2choi()");
896  // END EXCEPTION CHECKS
897 
898  cmat result(D * D, D * D);
899 
900 #ifdef WITH_OPENMP_
901 #pragma omp parallel for collapse(4)
902 #endif // WITH_OPENMP_
903  for (idx a = 0; a < D; ++a)
904  for (idx b = 0; b < D; ++b)
905  for (idx m = 0; m < D; ++m)
906  for (idx n = 0; n < D; ++n)
907  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
908 
909  return result;
910 }
911 
925 template <typename Derived>
926 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
927  const std::vector<idx>& dims) {
928  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
929 
930  // EXCEPTION CHECKS
931 
932  // check zero-size
934  throw exception::ZeroSize("qpp::ptrace1()");
935 
936  // check that dims is a valid dimension vector
937  if (!internal::check_dims(dims))
938  throw exception::DimsInvalid("qpp::ptrace1()");
939 
940  // check dims has only 2 elements
941  if (dims.size() != 2)
942  throw exception::NotBipartite("qpp::ptrace1()");
943  // END EXCEPTION CHECKS
944 
945  idx DA = dims[0];
946  idx DB = dims[1];
947 
950 
951  //************ ket ************//
952  if (internal::check_cvector(rA)) // we have a ket
953  {
954  // check that dims match the dimension of A
955  if (!internal::check_dims_match_cvect(dims, rA))
956  throw exception::DimsMismatchCvector("qpp::ptrace1()");
957 
958  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
959  typename Derived::Scalar sum = 0;
960  for (idx m = 0; m < DA; ++m)
961  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
962 
963  return sum;
964  }; /* end worker */
965 
966 #ifdef WITH_OPENMP_
967 #pragma omp parallel for collapse(2)
968 #endif // WITH_OPENMP_
969  // column major order for speed
970  for (idx j = 0; j < DB; ++j)
971  for (idx i = 0; i < DB; ++i)
972  result(i, j) = worker(i, j);
973 
974  return result;
975  }
976  //************ density matrix ************//
977  else if (internal::check_square_mat(rA)) // we have a density operator
978  {
979  // check that dims match the dimension of A
980  if (!internal::check_dims_match_mat(dims, rA))
981  throw exception::DimsMismatchMatrix("qpp::ptrace1()");
982 
983  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
984  typename Derived::Scalar sum = 0;
985  for (idx m = 0; m < DA; ++m)
986  sum += rA(m * DB + i, m * DB + j);
987 
988  return sum;
989  }; /* end worker */
990 
991 #ifdef WITH_OPENMP_
992 #pragma omp parallel for collapse(2)
993 #endif // WITH_OPENMP_
994  // column major order for speed
995  for (idx j = 0; j < DB; ++j)
996  for (idx i = 0; i < DB; ++i)
997  result(i, j) = worker(i, j);
998 
999  return result;
1000  }
1001  //************ Exception: not ket nor density matrix ************//
1002  else
1003  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1004 }
1005 
1019 template <typename Derived>
1020 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1021  idx d = 2) {
1022  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1023 
1024  // EXCEPTION CHECKS
1025 
1026  // check zero size
1028  throw exception::ZeroSize("qpp::ptrace1()");
1029 
1030  // check valid dims
1031  if (d == 0)
1032  throw exception::DimsInvalid("qpp::ptrace1()");
1033  // END EXCEPTION CHECKS
1034 
1035  std::vector<idx> dims(2, d); // local dimensions vector
1036 
1037  return ptrace1(A, dims);
1038 }
1039 
1053 template <typename Derived>
1054 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1055  const std::vector<idx>& dims) {
1056  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1057 
1058  // EXCEPTION CHECKS
1059 
1060  // check zero-size
1062  throw exception::ZeroSize("qpp::ptrace2()");
1063 
1064  // check that dims is a valid dimension vector
1065  if (!internal::check_dims(dims))
1066  throw exception::DimsInvalid("qpp::ptrace2()");
1067 
1068  // check dims has only 2 elements
1069  if (dims.size() != 2)
1070  throw exception::NotBipartite("qpp::ptrace2()");
1071  // END EXCEPTION CHECKS
1072 
1073  idx DA = dims[0];
1074  idx DB = dims[1];
1075 
1078 
1079  //************ ket ************//
1080  if (internal::check_cvector(rA)) // we have a ket
1081  {
1082  // check that dims match the dimension of A
1083  if (!internal::check_dims_match_cvect(dims, rA))
1084  throw exception::DimsMismatchCvector("qpp::ptrace2()");
1085 
1086  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
1087  typename Derived::Scalar sum = 0;
1088  for (idx m = 0; m < DB; ++m)
1089  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1090 
1091  return sum;
1092  }; /* end worker */
1093 
1094 #ifdef WITH_OPENMP_
1095 #pragma omp parallel for collapse(2)
1096 #endif // WITH_OPENMP_
1097  // column major order for speed
1098  for (idx j = 0; j < DA; ++j)
1099  for (idx i = 0; i < DA; ++i)
1100  result(i, j) = worker(i, j);
1101 
1102  return result;
1103  }
1104  //************ density matrix ************//
1105  else if (internal::check_square_mat(rA)) // we have a density operator
1106  {
1107  // check that dims match the dimension of A
1108  if (!internal::check_dims_match_mat(dims, rA))
1109  throw exception::DimsMismatchMatrix("qpp::ptrace2()");
1110 
1111 #ifdef WITH_OPENMP_
1112 #pragma omp parallel for collapse(2)
1113 #endif // WITH_OPENMP_
1114  // column major order for speed
1115  for (idx j = 0; j < DA; ++j)
1116  for (idx i = 0; i < DA; ++i)
1117  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1118 
1119  return result;
1120  }
1121  //************ Exception: not ket nor density matrix ************//
1122  else
1123  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1124 }
1125 
1139 template <typename Derived>
1140 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1141  idx d = 2) {
1142  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1143 
1144  // EXCEPTION CHECKS
1145 
1146  // check zero size
1148  throw exception::ZeroSize("qpp::ptrace2()");
1149 
1150  // check valid dims
1151  if (d == 0)
1152  throw exception::DimsInvalid("qpp::ptrace2()");
1153  // END EXCEPTION CHECKS
1154 
1155  std::vector<idx> dims(2, d); // local dimensions vector
1156 
1157  return ptrace2(A, dims);
1158 }
1159 
1174 template <typename Derived>
1175 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1176  const std::vector<idx>& subsys,
1177  const std::vector<idx>& dims) {
1178  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1179 
1180  // EXCEPTION CHECKS
1181 
1182  // check zero-size
1184  throw exception::ZeroSize("qpp::ptrace()");
1185 
1186  // check that dims is a valid dimension vector
1187  if (!internal::check_dims(dims))
1188  throw exception::DimsInvalid("qpp::ptrace()");
1189 
1190  // check that subsys are valid
1191  if (!internal::check_subsys_match_dims(subsys, dims))
1192  throw exception::SubsysMismatchDims("qpp::ptrace()");
1193  // END EXCEPTION CHECKS
1194 
1195  idx D = static_cast<idx>(rA.rows());
1196  idx n = dims.size();
1197  idx n_subsys = subsys.size();
1198  idx n_subsys_bar = n - n_subsys;
1199  idx Dsubsys = 1;
1200  for (idx i = 0; i < n_subsys; ++i)
1201  Dsubsys *= dims[subsys[i]];
1202  idx Dsubsys_bar = D / Dsubsys;
1203 
1204  idx Cdims[maxn];
1205  idx Csubsys[maxn];
1206  idx Cdimssubsys[maxn];
1207  idx Csubsys_bar[maxn];
1208  idx Cdimssubsys_bar[maxn];
1209 
1210  idx Cmidxcolsubsys_bar[maxn];
1211 
1212  std::vector<idx> subsys_bar = complement(subsys, n);
1213  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1214  std::begin(Csubsys_bar));
1215 
1216  for (idx i = 0; i < n; ++i) {
1217  Cdims[i] = dims[i];
1218  }
1219  for (idx i = 0; i < n_subsys; ++i) {
1220  Csubsys[i] = subsys[i];
1221  Cdimssubsys[i] = dims[subsys[i]];
1222  }
1223  for (idx i = 0; i < n_subsys_bar; ++i) {
1224  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1225  }
1226 
1228  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1229 
1230  //************ ket ************//
1231  if (internal::check_cvector(rA)) // we have a ket
1232  {
1233  // check that dims match the dimension of A
1234  if (!internal::check_dims_match_cvect(dims, rA))
1235  throw exception::DimsMismatchCvector("qpp::ptrace()");
1236 
1237  if (subsys.size() == dims.size()) {
1238  result(0, 0) = (adjoint(rA) * rA).value();
1239  return result;
1240  }
1241 
1242  if (subsys.size() == 0)
1243  return rA * adjoint(rA);
1244 
1245  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1246  // use static allocation for speed!
1247 
1248  idx Cmidxrow[maxn];
1249  idx Cmidxcol[maxn];
1250  idx Cmidxrowsubsys_bar[maxn];
1251  idx Cmidxsubsys[maxn];
1252 
1253  /* get the row multi-indexes of the complement */
1254  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1255  Cmidxrowsubsys_bar);
1256  /* write them in the global row/col multi-indexes */
1257  for (idx k = 0; k < n_subsys_bar; ++k) {
1258  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1259  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1260  }
1261  typename Derived::Scalar sm = 0;
1262  for (idx a = 0; a < Dsubsys; ++a) {
1263  // get the multi-index over which we do the summation
1264  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1265  // write it into the global row/col multi-indexes
1266  for (idx k = 0; k < n_subsys; ++k)
1267  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1268  Cmidxsubsys[k];
1269 
1270  // now do the sum
1271  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims)) *
1272  std::conj(rA(internal::multiidx2n(Cmidxcol, n, Cdims)));
1273  }
1274 
1275  return sm;
1276  }; /* end worker */
1277 
1278  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1279  {
1280  // compute the column multi-indexes of the complement
1281  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1282  Cmidxcolsubsys_bar);
1283 #ifdef WITH_OPENMP_
1284 #pragma omp parallel for
1285 #endif // WITH_OPENMP_
1286  for (idx i = 0; i < Dsubsys_bar; ++i) {
1287  result(i, j) = worker(i);
1288  }
1289  }
1290 
1291  return result;
1292  }
1293  //************ density matrix ************//
1294  else if (internal::check_square_mat(rA)) // we have a density operator
1295  {
1296  // check that dims match the dimension of A
1297  if (!internal::check_dims_match_mat(dims, rA))
1298  throw exception::DimsMismatchMatrix("qpp::ptrace()");
1299 
1300  if (subsys.size() == dims.size()) {
1301  result(0, 0) = rA.trace();
1302  return result;
1303  }
1304 
1305  if (subsys.size() == 0)
1306  return rA;
1307 
1308  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1309  // use static allocation for speed!
1310 
1311  idx Cmidxrow[maxn];
1312  idx Cmidxcol[maxn];
1313  idx Cmidxrowsubsys_bar[maxn];
1314  idx Cmidxsubsys[maxn];
1315 
1316  /* get the row/col multi-indexes of the complement */
1317  internal::n2multiidx(i, n_subsys_bar, Cdimssubsys_bar,
1318  Cmidxrowsubsys_bar);
1319  /* write them in the global row/col multi-indexes */
1320  for (idx k = 0; k < n_subsys_bar; ++k) {
1321  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1322  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1323  }
1324  typename Derived::Scalar sm = 0;
1325  for (idx a = 0; a < Dsubsys; ++a) {
1326  // get the multi-index over which we do the summation
1327  internal::n2multiidx(a, n_subsys, Cdimssubsys, Cmidxsubsys);
1328  // write it into the global row/col multi-indexes
1329  for (idx k = 0; k < n_subsys; ++k)
1330  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1331  Cmidxsubsys[k];
1332 
1333  // now do the sum
1334  sm += rA(internal::multiidx2n(Cmidxrow, n, Cdims),
1335  internal::multiidx2n(Cmidxcol, n, Cdims));
1336  }
1337 
1338  return sm;
1339  }; /* end worker */
1340 
1341  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1342  {
1343  // compute the column multi-indexes of the complement
1344  internal::n2multiidx(j, n_subsys_bar, Cdimssubsys_bar,
1345  Cmidxcolsubsys_bar);
1346 #ifdef WITH_OPENMP_
1347 #pragma omp parallel for
1348 #endif // WITH_OPENMP_
1349  for (idx i = 0; i < Dsubsys_bar; ++i) {
1350  result(i, j) = worker(i);
1351  }
1352  }
1353 
1354  return result;
1355  }
1356  //************ Exception: not ket nor density matrix ************//
1357  else
1358  throw exception::MatrixNotSquareNorCvector("qpp::ptrace()");
1359 }
1360 
1375 template <typename Derived>
1376 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1377  const std::vector<idx>& subsys,
1378  idx d = 2) {
1379  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1380 
1381  // EXCEPTION CHECKS
1382 
1383  // check zero size
1385  throw exception::ZeroSize("qpp::ptrace()");
1386 
1387  // check valid dims
1388  if (d < 2)
1389  throw exception::DimsInvalid("qpp::ptrace()");
1390  // END EXCEPTION CHECKS
1391 
1392  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1393  std::vector<idx> dims(n, d); // local dimensions vector
1394 
1395  return ptrace(rA, subsys, dims);
1396 }
1397 
1411 template <typename Derived>
1413 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& subsys,
1414  const std::vector<idx>& dims) {
1415  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1416 
1417  // EXCEPTION CHECKS
1418 
1419  // check zero-size
1421  throw exception::ZeroSize("qpp::ptranspose()");
1422 
1423  // check that dims is a valid dimension vector
1424  if (!internal::check_dims(dims))
1425  throw exception::DimsInvalid("qpp::ptranspose()");
1426 
1427  // check that subsys are valid
1428  if (!internal::check_subsys_match_dims(subsys, dims))
1429  throw exception::SubsysMismatchDims("qpp::ptranspose()");
1430  // END EXCEPTION CHECKS
1431 
1432  idx D = static_cast<idx>(rA.rows());
1433  idx n = dims.size();
1434  idx n_subsys = subsys.size();
1435  idx Cdims[maxn];
1436  idx Cmidxcol[maxn];
1437  idx Csubsys[maxn];
1438 
1439  // copy dims in Cdims and subsys in Csubsys
1440  for (idx i = 0; i < n; ++i)
1441  Cdims[i] = dims[i];
1442  for (idx i = 0; i < n_subsys; ++i)
1443  Csubsys[i] = subsys[i];
1444 
1445  dyn_mat<typename Derived::Scalar> result(D, D);
1446 
1447  //************ ket ************//
1448  if (internal::check_cvector(rA)) // we have a ket
1449  {
1450  // check that dims match the dimension of A
1451  if (!internal::check_dims_match_cvect(dims, rA))
1452  throw exception::DimsMismatchCvector("qpp::ptranspose()");
1453 
1454  if (subsys.size() == dims.size())
1455  return (rA * adjoint(rA)).transpose();
1456 
1457  if (subsys.size() == 0)
1458  return rA * adjoint(rA);
1459 
1460  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1461  // use static allocation for speed!
1462  idx midxcoltmp[maxn];
1463  idx midxrow[maxn];
1464 
1465  for (idx k = 0; k < n; ++k)
1466  midxcoltmp[k] = Cmidxcol[k];
1467 
1468  /* compute the row multi-index */
1469  internal::n2multiidx(i, n, Cdims, midxrow);
1470 
1471  for (idx k = 0; k < n_subsys; ++k)
1472  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1473 
1474  /* writes the result */
1475  return rA(internal::multiidx2n(midxrow, n, Cdims)) *
1476  std::conj(rA(internal::multiidx2n(midxcoltmp, n, Cdims)));
1477  }; /* end worker */
1478 
1479  for (idx j = 0; j < D; ++j) {
1480  // compute the column multi-index
1481  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1482 
1483 #ifdef WITH_OPENMP_
1484 #pragma omp parallel for
1485 #endif // WITH_OPENMP_
1486  for (idx i = 0; i < D; ++i)
1487  result(i, j) = worker(i);
1488  }
1489 
1490  return result;
1491  }
1492  //************ density matrix ************//
1493  else if (internal::check_square_mat(rA)) // we have a density operator
1494  {
1495  // check that dims match the dimension of A
1496  if (!internal::check_dims_match_mat(dims, rA))
1497  throw exception::DimsMismatchMatrix("qpp::ptranspose()");
1498 
1499  if (subsys.size() == dims.size())
1500  return rA.transpose();
1501 
1502  if (subsys.size() == 0)
1503  return rA;
1504 
1505  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1506  // use static allocation for speed!
1507  idx midxcoltmp[maxn];
1508  idx midxrow[maxn];
1509 
1510  for (idx k = 0; k < n; ++k)
1511  midxcoltmp[k] = Cmidxcol[k];
1512 
1513  /* compute the row multi-index */
1514  internal::n2multiidx(i, n, Cdims, midxrow);
1515 
1516  for (idx k = 0; k < n_subsys; ++k)
1517  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1518 
1519  /* writes the result */
1520  return rA(internal::multiidx2n(midxrow, n, Cdims),
1521  internal::multiidx2n(midxcoltmp, n, Cdims));
1522  }; /* end worker */
1523 
1524  for (idx j = 0; j < D; ++j) {
1525  // compute the column multi-index
1526  internal::n2multiidx(j, n, Cdims, Cmidxcol);
1527 
1528 #ifdef WITH_OPENMP_
1529 #pragma omp parallel for
1530 #endif // WITH_OPENMP_
1531  for (idx i = 0; i < D; ++i)
1532  result(i, j) = worker(i);
1533  }
1534 
1535  return result;
1536  }
1537  //************ Exception: not ket nor density matrix ************//
1538  else
1539  throw exception::MatrixNotSquareNorCvector("qpp::ptranspose()");
1540 }
1541 
1555 template <typename Derived>
1557 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& subsys,
1558  idx d = 2) {
1559  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1560 
1561  // EXCEPTION CHECKS
1562 
1563  // check zero size
1565  throw exception::ZeroSize("qpp::ptranspose()");
1566 
1567  // check valid dims
1568  if (d < 2)
1569  throw exception::DimsInvalid("qpp::ptranspose()");
1570  // END EXCEPTION CHECKS
1571 
1572  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1573  std::vector<idx> dims(n, d); // local dimensions vector
1574 
1575  return ptranspose(rA, subsys, dims);
1576 }
1577 
1590 template <typename Derived>
1592 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1593  const std::vector<idx>& dims) {
1594  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1595 
1596  // EXCEPTION CHECKS
1597 
1598  // check zero-size
1600  throw exception::ZeroSize("qpp::syspermute()");
1601 
1602  // check that dims is a valid dimension vector
1603  if (!internal::check_dims(dims))
1604  throw exception::DimsInvalid("qpp::syspermute()");
1605 
1606  // check that we have a valid permutation
1607  if (!internal::check_perm(perm))
1608  throw exception::PermInvalid("qpp::syspermute()");
1609 
1610  // check that permutation match dimensions
1611  if (perm.size() != dims.size())
1612  throw exception::PermMismatchDims("qpp::syspermute()");
1613  // END EXCEPTION CHECKS
1614 
1615  idx D = static_cast<idx>(rA.rows());
1616  idx n = dims.size();
1617 
1619 
1620  //************ ket ************//
1621  if (internal::check_cvector(rA)) // we have a column vector
1622  {
1623  idx Cdims[maxn];
1624  idx Cperm[maxn];
1625 
1626  // check that dims match the dimension of rA
1627  if (!internal::check_dims_match_cvect(dims, rA))
1628  throw exception::DimsMismatchCvector("qpp::syspermute()");
1629 
1630  // copy dims in Cdims and perm in Cperm
1631  for (idx i = 0; i < n; ++i) {
1632  Cdims[i] = dims[i];
1633  Cperm[i] = perm[i];
1634  }
1635  result.resize(D, 1);
1636 
1637  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1638  // use static allocation for speed,
1639  // double the size for matrices reshaped as vectors
1640  idx midx[maxn];
1641  idx midxtmp[maxn];
1642  idx permdims[maxn];
1643 
1644  /* compute the multi-index */
1645  internal::n2multiidx(i, n, Cdims, midx);
1646 
1647  for (idx k = 0; k < n; ++k) {
1648  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1649  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1650  }
1651 
1652  return internal::multiidx2n(midxtmp, n, permdims);
1653  }; /* end worker */
1654 
1655 #ifdef WITH_OPENMP_
1656 #pragma omp parallel for
1657 #endif // WITH_OPENMP_
1658  for (idx i = 0; i < D; ++i)
1659  result(worker(i)) = rA(i);
1660 
1661  return result;
1662  }
1663  //************ density matrix ************//
1664  else if (internal::check_square_mat(rA)) // we have a density operator
1665  {
1666  idx Cdims[2 * maxn];
1667  idx Cperm[2 * maxn];
1668 
1669  // check that dims match the dimension of rA
1670  if (!internal::check_dims_match_mat(dims, rA))
1671  throw exception::DimsMismatchMatrix("qpp::syspermute()");
1672 
1673  // copy dims in Cdims and perm in Cperm
1674  for (idx i = 0; i < n; ++i) {
1675  Cdims[i] = dims[i];
1676  Cdims[i + n] = dims[i];
1677  Cperm[i] = perm[i];
1678  Cperm[i + n] = perm[i] + n;
1679  }
1680  result.resize(D * D, 1);
1681  // map A to a column vector
1683  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1684  const_cast<typename Derived::Scalar*>(rA.data()), D * D, 1);
1685 
1686  auto worker = [&Cdims, &Cperm, n ](idx i) noexcept->idx {
1687  // use static allocation for speed,
1688  // double the size for matrices reshaped as vectors
1689  idx midx[2 * maxn];
1690  idx midxtmp[2 * maxn];
1691  idx permdims[2 * maxn];
1692 
1693  /* compute the multi-index */
1694  internal::n2multiidx(i, 2 * n, Cdims, midx);
1695 
1696  for (idx k = 0; k < 2 * n; ++k) {
1697  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1698  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1699  }
1700 
1701  return internal::multiidx2n(midxtmp, 2 * n, permdims);
1702  }; /* end worker */
1703 
1704 #ifdef WITH_OPENMP_
1705 #pragma omp parallel for
1706 #endif // WITH_OPENMP_
1707  for (idx i = 0; i < D * D; ++i)
1708  result(worker(i)) = rA(i);
1709 
1710  return reshape(result, D, D);
1711  }
1712  //************ Exception: not ket nor density matrix ************//
1713  else
1714  throw exception::MatrixNotSquareNorCvector("qpp::syspermute()");
1715 }
1716 
1729 template <typename Derived>
1731 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1732  idx d = 2) {
1733  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1734 
1735  // EXCEPTION CHECKS
1736 
1737  // check zero size
1739  throw exception::ZeroSize("qpp::syspermute()");
1740 
1741  // check valid dims
1742  if (d < 2)
1743  throw exception::DimsInvalid("qpp::syspermute()");
1744  // END EXCEPTION CHECKS
1745 
1746  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1747  std::vector<idx> dims(n, d); // local dimensions vector
1748 
1749  return syspermute(rA, perm, dims);
1750 }
1751 
1752 // as in https://arxiv.org/abs/1707.08834
1763 template <typename Derived>
1764 dyn_mat<typename Derived::Scalar> applyQFT(const Eigen::MatrixBase<Derived>& A,
1765  const std::vector<idx>& subsys,
1766  idx d = 2) {
1767  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1768 
1769  // EXCEPTION CHECKS
1770 
1771  // check zero sizes
1773  throw exception::ZeroSize("qpp::applyQFT()");
1774 
1775  // check valid subsystem dimension
1776  if (d < 2)
1777  throw exception::DimsInvalid("qpp::applyQFT()");
1778 
1779  // total number of qubits/qudits in the state
1780  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1781 
1782  std::vector<idx> dims(n, d); // local dimensions vector
1783  // check subsys is valid w.r.t. dims
1784  if (!internal::check_subsys_match_dims(subsys, dims))
1785  throw exception::SubsysMismatchDims("qpp::applyQFT()");
1786 
1787  //************ ket ************//
1788  if (internal::check_cvector(rA)) // we have a ket
1789  {
1790  // check that dims match state vector
1791  if (!internal::check_dims_match_cvect(dims, rA))
1792  throw exception::DimsMismatchCvector("qpp::applyQFT()");
1793  }
1794  //************ density matrix ************//
1795  else if (internal::check_square_mat(rA)) // we have a density operator
1796  {
1797  // check that dims match state matrix
1798  if (!internal::check_dims_match_mat(dims, rA))
1799  throw exception::DimsMismatchMatrix("qpp::applyQFT()");
1800  }
1801  //************ Exception: not ket nor density matrix ************//
1802  else
1803  throw exception::MatrixNotSquareNorCvector("qpp::applyQFT()");
1804  // END EXCEPTION CHECKS
1805 
1807 
1808  idx n_subsys = subsys.size();
1809 
1810  if (d == 2) // qubits
1811  {
1812  for (idx i = 0; i < n_subsys; ++i) {
1813  // apply qudit Fourier on qubit i
1814  result = apply(result, Gates::get_instance().H, {subsys[i]});
1815  // apply controlled rotations
1816  for (idx j = 2; j <= n_subsys - i; ++j) {
1817  // construct Rj
1818  cmat Rj(2, 2);
1819  Rj << 1, 0, 0, omega(std::pow(2, j));
1820  result =
1821  applyCTRL(result, Rj, {subsys[i + j - 1]}, {subsys[i]});
1822  }
1823  }
1824  // we have the qubits in reversed order, we must swap them
1825  for (idx i = 0; i < n_subsys / 2; ++i) {
1826  result = apply(result, Gates::get_instance().SWAP,
1827  {subsys[i], subsys[n_subsys - i - 1]});
1828  }
1829 
1830  } else { // qudits
1831  for (idx i = 0; i < n_subsys; ++i) {
1832  // apply qudit Fourier on qudit i
1833  result = apply(result, Gates::get_instance().Fd(d), {subsys[i]}, d);
1834  // apply controlled rotations
1835  for (idx j = 2; j <= n_subsys - i; ++j) {
1836  // construct Rj
1837  cmat Rj = cmat::Zero(d, d);
1838  for (idx m = 0; m < d; ++m) {
1839  Rj(m, m) = exp(2.0 * pi * m * 1_i / std::pow(d, j));
1840  }
1841  result =
1842  applyCTRL(result, Rj, {subsys[i + j - 1]}, {subsys[i]}, d);
1843  }
1844  }
1845  // we have the qudits in reversed order, we must swap them
1846  for (idx i = 0; i < n_subsys / 2; ++i) {
1847  result = apply(result, Gates::get_instance().SWAPd(d),
1848  {subsys[i], subsys[n_subsys - i - 1]}, d);
1849  }
1850  }
1851 
1852  return result;
1853 }
1854 
1855 // as in https://arxiv.org/abs/1707.08834
1863 template <typename Derived>
1864 dyn_col_vect<typename Derived::Scalar> QFT(const Eigen::MatrixBase<Derived>& A,
1865  idx d = 2) {
1866  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1867 
1868  // EXCEPTION CHECKS
1869 
1870  // check zero-size
1872  throw exception::ZeroSize("qpp::QFT()");
1873 
1874  // check valid subsystem dimension
1875  if (d < 2)
1876  throw exception::DimsInvalid("qpp::QFT()");
1877 
1878  // total number of qubits/qudits in the state
1879  idx n = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1880 
1881  std::vector<idx> dims(n, d); // local dimensions vector
1882 
1883  //************ ket ************//
1884  if (internal::check_cvector(rA)) // we have a ket
1885  {
1886  // check that dims match state vector
1887  if (!internal::check_dims_match_cvect(dims, rA))
1888  throw exception::DimsMismatchCvector("qpp::QFT()");
1889  }
1890  //************ density matrix ************//
1891  else if (internal::check_square_mat(rA)) // we have a density operator
1892  {
1893  // check that dims match state matrix
1894  if (!internal::check_dims_match_mat(dims, rA))
1895  throw exception::DimsMismatchMatrix("qpp::QFT()");
1896  }
1897  //************ Exception: not ket nor density matrix ************//
1898  else
1899  throw exception::MatrixNotSquareNorCvector("qpp::QFT()");
1900  // END EXCEPTION CHECKS
1901 
1902  std::vector<idx> subsys(n);
1903  std::iota(std::begin(subsys), std::end(subsys), 0);
1904  ket result = applyQFT(rA, subsys, d);
1905 
1906  return result;
1907 }
1908 
1909 } /* namespace qpp */
1910 
1911 #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:75
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
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:818
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
cplx omega(idx D)
D-th root of unity.
Definition: constants.h:97
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:68
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
dyn_mat< typename Derived::Scalar > syspermute(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &perm, const std::vector< idx > &dims)
Subsystem permutation.
Definition: operations.h:1592
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
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
dyn_mat< typename Derived1::Scalar > apply(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Applies the gate A to the part subsys of the multi-partite state vector or density matrix state...
Definition: operations.h:440
Quantum++ main namespace.
Definition: codes.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:259
idx get_num_subsys(idx sz, idx d)
Definition: util.h:366
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:376
Invalid dimension(s) exception.
Definition: exception.h:269
dyn_col_vect< typename Derived::Scalar > QFT(const Eigen::MatrixBase< Derived > &A, idx d=2)
Qudit quantum Fourier transform.
Definition: operations.h:1864
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:851
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:885
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1733
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 Derived::Scalar > applyQFT(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, idx d=2)
Applies the qudit quantum Fourier transform to the part subsys of the multi-partite state vector or d...
Definition: operations.h:1764
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:926
Invalid permutation exception.
Definition: exception.h:380
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:83
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:868
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1054
Matrix mismatch subsystems exception.
Definition: exception.h:254
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_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:751
Type mismatch exception.
Definition: exception.h:530
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:378
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:401
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:763
dyn_mat< typename Derived::Scalar > ptrace(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1175
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
constexpr double pi
Definition: constants.h:80
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:692
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 > &subsys, const std::vector< idx > &dims)
Applies the controlled-gate A to the part subsys of the multi-partite state vector or density matrix ...
Definition: operations.h:55
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
static const Gates & get_instance() noexcept(std::is_nothrow_constructible< const Gates >::value)
Definition: singleton.h:92
dyn_mat< typename Derived::Scalar > ptranspose(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial transpose.
Definition: operations.h:1413
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
std::size_t idx
Non-negative integer index.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:149
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:201
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1095
Object has zero size exception.
Definition: exception.h:134