Quantum++  v1.0-rc2
A modern C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef OPERATIONS_H_
28 #define OPERATIONS_H_
29 
30 namespace qpp
31 {
32 
49 template<typename Derived1, typename Derived2>
51  const Eigen::MatrixBase<Derived1>& state,
52  const Eigen::MatrixBase<Derived2>& A,
53  const std::vector<idx>& ctrl,
54  const std::vector<idx>& subsys,
55  const std::vector<idx>& dims)
56 {
57  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
58  = state.derived();
59  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
60 
61  // EXCEPTION CHECKS
62 
63  // check types
64  if (!std::is_same<typename Derived1::Scalar,
65  typename Derived2::Scalar>::value)
66  throw exception::TypeMismatch("qpp::applyCTRL()");
67 
68  // check zero sizes
70  throw exception::ZeroSize("qpp::applyCTRL()");
71 
72  // check zero sizes
73  if (!internal::check_nonzero_size(rstate))
74  throw exception::ZeroSize("qpp::applyCTRL()");
75 
76  // check square matrix for the gate
78  throw exception::MatrixNotSquare("qpp::applyCTRL()");
79 
80  // check that all control subsystems have the same dimension
81  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
82  for (idx i = 1; i < ctrl.size(); ++i)
83  if (dims[ctrl[i]] != d)
84  throw exception::DimsNotEqual("qpp::applyCTRL()");
85 
86  // check that dimension is valid
87  if (!internal::check_dims(dims))
88  throw exception::DimsInvalid("qpp::applyCTRL()");
89 
90  // check subsys is valid w.r.t. dims
91  if (!internal::check_subsys_match_dims(subsys, dims))
92  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
93 
94  // check that gate matches the dimensions of the subsys
95  std::vector<idx> subsys_dims(subsys.size());
96  for (idx i = 0; i < subsys.size(); ++i)
97  subsys_dims[i] = dims[subsys[i]];
98  if (!internal::check_dims_match_mat(subsys_dims, rA))
99  throw exception::MatrixMismatchSubsys("qpp::applyCTRL()");
100 
101  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
102  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
103  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
104 
105  // check that ctrl + gate subsystem is valid
106  // with respect to local dimensions
107  if (!internal::check_subsys_match_dims(ctrlgate, dims))
108  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
109  // END EXCEPTION CHECKS
110 
111  // construct the table of A^i and (A^dagger)^i
112  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
113  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
114  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i)
115  {
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  {
157  idx indx = 0;
158  typename Derived1::Scalar coeff = 0;
159 
160  idx Cmidx[maxn]; // the total multi-index
161  idx CmidxA[maxn]; // the gate part multi-index
162  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
163 
164  // compute the index
165 
166  // set the CTRL part
167  for (idx k = 0; k < ctrlsize; ++k)
168  {
169  Cmidx[ctrl[k]] = i_;
170  }
171 
172  // set the rest
173  internal::n2multiidx(r_, N - ctrlgatesize,
174  CdimsCTRLA_bar, CmidxCTRLA_bar);
175  for (idx k = 0; k < N - ctrlgatesize; ++k)
176  {
177  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
178  }
179 
180  // set the A part
181  internal::n2multiidx(m_, subsyssize, CdimsA, CmidxA);
182  for (idx k = 0; k < subsyssize; ++k)
183  {
184  Cmidx[subsys[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  {
193  internal::n2multiidx(n_, subsyssize, CdimsA, CmidxA);
194  for (idx k = 0; k < subsyssize; ++k)
195  {
196  Cmidx[subsys[k]] = CmidxA[k];
197  }
198  coeff += Ai[i_](m_, n_) *
199  rstate(internal::multiidx2n(Cmidx, N, Cdims));
200  }
201 
202  return std::make_pair(coeff, indx);
203  }; /* end coeff_idx_ket */
204 
205  // worker, computes the coefficient and the index
206  // for the density matrix case
207  // used in #pragma omp parallel for collapse
208  auto coeff_idx_rho = [&](idx i1_, idx m1_,
209  idx r1_, idx i2_, idx m2_, idx r2_) noexcept
210  -> std::tuple<typename Derived1::Scalar, idx, idx>
211  {
212  idx idxrow = 0;
213  idx idxcol = 0;
214  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
215 
216  idx Cmidxrow[maxn]; // the total row multi-index
217  idx Cmidxcol[maxn]; // the total col multi-index
218  idx CmidxArow[maxn]; // the gate part row multi-index
219  idx CmidxAcol[maxn]; // the gate part col multi-index
220  idx CmidxCTRLrow[maxn]; // the control row multi-index
221  idx CmidxCTRLcol[maxn]; // the control col multi-index
222  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
223  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
224 
225  // compute the ket/bra indexes
226 
227  // set the CTRL part
228  internal::n2multiidx(i1_, ctrlsize,
229  CdimsCTRL, CmidxCTRLrow);
230  internal::n2multiidx(i2_, ctrlsize,
231  CdimsCTRL, CmidxCTRLcol);
232 
233  for (idx k = 0; k < ctrlsize; ++k)
234  {
235  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
236  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
237  }
238 
239  // set the rest
240  internal::n2multiidx(r1_, N - ctrlgatesize,
241  CdimsCTRLA_bar, CmidxCTRLA_barrow);
242  internal::n2multiidx(r2_, N - ctrlgatesize,
243  CdimsCTRLA_bar, CmidxCTRLA_barcol);
244  for (idx k = 0; k < N - ctrlgatesize; ++k)
245  {
246  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
247  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
248  }
249 
250  // set the A part
251  internal::n2multiidx(m1_, subsyssize, CdimsA, CmidxArow);
252  internal::n2multiidx(m2_, subsyssize, CdimsA, CmidxAcol);
253  for (idx k = 0; k < subsys.size(); ++k)
254  {
255  Cmidxrow[subsys[k]] = CmidxArow[k];
256  Cmidxcol[subsys[k]] = CmidxAcol[k];
257  }
258 
259  // we now got the total row/col indexes
260  idxrow = internal::multiidx2n(Cmidxrow, N, Cdims);
261  idxcol = internal::multiidx2n(Cmidxcol, N, Cdims);
262 
263  // check whether all CTRL row and col multi indexes are equal
264  bool all_ctrl_rows_equal = true;
265  bool all_ctrl_cols_equal = true;
266 
267  idx first_ctrl_row, first_ctrl_col;
268  if (ctrlsize > 0)
269  {
270  first_ctrl_row = CmidxCTRLrow[0];
271  first_ctrl_col = CmidxCTRLcol[0];
272  } else
273  {
274  first_ctrl_row = first_ctrl_col = 1;
275  }
276 
277  for (idx k = 1; k < ctrlsize; ++k)
278  {
279  if (CmidxCTRLrow[k] != first_ctrl_row)
280  {
281  all_ctrl_rows_equal = false;
282  break;
283  }
284  }
285  for (idx k = 1; k < ctrlsize; ++k)
286  {
287  if (CmidxCTRLcol[k] != first_ctrl_col)
288  {
289  all_ctrl_cols_equal = false;
290  break;
291  }
292  }
293 
294  // at least one control activated, compute the coefficient
295  for (idx n1_ = 0; n1_ < DA; ++n1_)
296  {
297  internal::n2multiidx(n1_, subsyssize, CdimsA, CmidxArow);
298  for (idx k = 0; k < subsyssize; ++k)
299  {
300  Cmidxrow[subsys[k]] = CmidxArow[k];
301  }
302  idx idxrowtmp = internal::multiidx2n(Cmidxrow, N, Cdims);
303 
304  if (all_ctrl_rows_equal)
305  {
306  lhs = Ai[first_ctrl_row](m1_, n1_);
307  } else
308  {
309  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
310  }
311 
312  for (idx n2_ = 0; n2_ < DA; ++n2_)
313  {
314  internal::n2multiidx(n2_, subsyssize, CdimsA, CmidxAcol);
315  for (idx k = 0; k < subsyssize; ++k)
316  {
317  Cmidxcol[subsys[k]] = CmidxAcol[k];
318  }
319 
320  if (all_ctrl_cols_equal)
321  {
322  rhs = Aidagger[first_ctrl_col](n2_, m2_);
323  } else
324  {
325  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
326  }
327 
328  idx idxcoltmp = internal::multiidx2n(Cmidxcol, N, Cdims);
329 
330  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
331  }
332  }
333 
334  return std::make_tuple(coeff, idxrow, idxcol);
335  }; /* end coeff_idx_rho */
336 
337  //************ ket ************//
338  if (internal::check_cvector(rstate)) // we have a ket
339  {
340  // check that dims match state vector
341  if (!internal::check_dims_match_cvect(dims, rstate))
342  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
343  if (D == 1)
344  return rstate;
345 
346  dyn_mat<typename Derived1::Scalar> result = rstate;
347 
348 #ifdef WITH_OPENMP_
349 #pragma omp parallel for collapse(2)
350 #endif // WITH_OPENMP_
351  for (idx m = 0; m < DA; ++m)
352  for (idx r = 0; r < DCTRLA_bar; ++r)
353  {
354  if (ctrlsize == 0) // no control
355  {
356  result(coeff_idx_ket(1, m, r).second) =
357  coeff_idx_ket(1, m, r).first;
358  } else
359  for (idx i = 0; i < d; ++i)
360  {
361  result(coeff_idx_ket(i, m, r).second) =
362  coeff_idx_ket(i, m, r).first;
363  }
364  }
365 
366  return result;
367  }
368  //************ density matrix ************//
369  else if (internal::check_square_mat(rstate)) // we have a density operator
370  {
371  // check that dims match state matrix
372  if (!internal::check_dims_match_mat(dims, rstate))
373  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
374 
375  if (D == 1)
376  return rstate;
377 
378  dyn_mat<typename Derived1::Scalar> result = rstate;
379 
380 #ifdef WITH_OPENMP_
381 #pragma omp parallel for collapse(4)
382 #endif // WITH_OPENMP_
383  for (idx m1 = 0; m1 < DA; ++m1)
384  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
385  for (idx m2 = 0; m2 < DA; ++m2)
386  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
387  if (ctrlsize == 0) // no control
388  {
389  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
390  1, m2, r2);
391  result(std::get<1>(coeff_idxes),
392  std::get<2>(coeff_idxes)) =
393  std::get<0>(coeff_idxes);
394  } else
395  {
396  for (idx i1 = 0; i1 < Dctrl; ++i1)
397  for (idx i2 = 0; i2 < Dctrl; ++i2)
398  {
399  auto coeff_idxes = coeff_idx_rho(
400  i1, m1, r1,
401  i2, m2, r2);
402  result(std::get<1>(coeff_idxes),
403  std::get<2>(coeff_idxes)) =
404  std::get<0>(coeff_idxes);
405  }
406  }
407 
408  return result;
409  }
410  //************ Exception: not ket nor density matrix ************//
411  else
412  throw exception::MatrixNotSquareNorCvector("qpp::applyCTRL()");
413 }
414 
430 template<typename Derived1, typename Derived2>
432  const Eigen::MatrixBase<Derived1>& state,
433  const Eigen::MatrixBase<Derived2>& A,
434  const std::vector<idx>& ctrl,
435  const std::vector<idx>& subsys,
436  idx d = 2)
437 {
438  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
439  = state.derived();
440  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
441 
442  // EXCEPTION CHECKS
443 
444  // check zero size
445  if (!internal::check_nonzero_size(rstate))
446  throw exception::ZeroSize("qpp::applyCTRL()");
447 
448  // check valid dims
449  if (d < 2)
450  throw exception::DimsInvalid("qpp::applyCTRL()");
451  // END EXCEPTION CHECKS
452 
453  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
454  std::vector<idx> dims(N, d); // local dimensions vector
455 
456  return applyCTRL(rstate, rA, ctrl, subsys, dims);
457 }
458 
472 template<typename Derived1, typename Derived2>
474  const Eigen::MatrixBase<Derived1>& state,
475  const Eigen::MatrixBase<Derived2>& A,
476  const std::vector<idx>& subsys,
477  const std::vector<idx>& dims)
478 {
479  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
480  = state.derived();
481  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
482 
483  // EXCEPTION CHECKS
484 
485  // check types
486  if (!std::is_same<typename Derived1::Scalar,
487  typename Derived2::Scalar>::value)
488  throw exception::TypeMismatch("qpp::apply()");
489 
490  // check zero sizes
492  throw exception::ZeroSize("qpp::apply()");
493 
494  // check zero sizes
495  if (!internal::check_nonzero_size(rstate))
496  throw exception::ZeroSize("qpp::apply()");
497 
498  // check square matrix for the gate
500  throw exception::MatrixNotSquare("qpp::apply()");
501 
502  // check that dimension is valid
503  if (!internal::check_dims(dims))
504  throw exception::DimsInvalid("qpp::apply()");
505 
506  // check subsys is valid w.r.t. dims
507  if (!internal::check_subsys_match_dims(subsys, dims))
508  throw exception::SubsysMismatchDims("qpp::apply()");
509 
510  // check that gate matches the dimensions of the subsys
511  std::vector<idx> subsys_dims(subsys.size());
512  for (idx i = 0; i < subsys.size(); ++i)
513  subsys_dims[i] = dims[subsys[i]];
514  if (!internal::check_dims_match_mat(subsys_dims, rA))
515  throw exception::MatrixMismatchSubsys("qpp::apply()");
516  // END EXCEPTION CHECKS
517 
518  //************ ket ************//
519  if (internal::check_cvector(rstate)) // we have a ket
520  {
521  // check that dims match state vector
522  if (!internal::check_dims_match_cvect(dims, rstate))
523  throw exception::DimsMismatchCvector("qpp::apply()");
524 
525  return applyCTRL(rstate, rA, {}, subsys, dims);
526  }
527  //************ density matrix ************//
528  else if (internal::check_square_mat(rstate)) // we have a density operator
529  {
530 
531  // check that dims match state matrix
532  if (!internal::check_dims_match_mat(dims, rstate))
533  throw exception::DimsMismatchMatrix("qpp::apply()");
534 
535  return applyCTRL(rstate, rA, {}, subsys, dims);
536  }
537  //************ Exception: not ket nor density matrix ************//
538  else
539  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
540 }
541 
555 template<typename Derived1, typename Derived2>
557  const Eigen::MatrixBase<Derived1>& state,
558  const Eigen::MatrixBase<Derived2>& A,
559  const std::vector<idx>& subsys,
560  idx d = 2)
561 {
562  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
563  = state.derived();
564  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
565 
566  // EXCEPTION CHECKS
567 
568  // check zero size
569  if (!internal::check_nonzero_size(rstate))
570  throw exception::ZeroSize("qpp::apply()");
571 
572  // check valid dims
573  if (d < 2)
574  throw exception::DimsInvalid("qpp::apply()");
575  // END EXCEPTION CHECKS
576 
577  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
578  std::vector<idx> dims(N, d); // local dimensions vector
579 
580  return apply(rstate, rA, subsys, dims);
581 }
582 
591 template<typename Derived>
592 cmat apply(const Eigen::MatrixBase<Derived>& A,
593  const std::vector<cmat>& Ks)
594 {
595  const cmat& rA = A.derived();
596 
597  // EXCEPTION CHECKS
598 
600  throw exception::ZeroSize("qpp::apply()");
602  throw exception::MatrixNotSquare("qpp::apply()");
603  if (Ks.size() == 0)
604  throw exception::ZeroSize("qpp::apply()");
605  if (!internal::check_square_mat(Ks[0]))
606  throw exception::MatrixNotSquare("qpp::apply()");
607  if (Ks[0].rows() != rA.rows())
608  throw exception::DimsMismatchMatrix("qpp::apply()");
609  for (auto&& it : Ks)
610  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
611  throw exception::DimsNotEqual("qpp::apply()");
612  // END EXCEPTION CHECKS
613 
614  cmat result = cmat::Zero(rA.rows(), rA.rows());
615 
616 #ifdef WITH_OPENMP_
617 #pragma omp parallel for
618 #endif // WITH_OPENMP_
619  for (idx i = 0; i < Ks.size(); ++i)
620  {
621 #ifdef WITH_OPENMP_
622 #pragma omp critical
623 #endif // WITH_OPENMP_
624  {
625  result += Ks[i] * rA * adjoint(Ks[i]);
626  }
627  }
628 
629  return result;
630 }
631 
642 template<typename Derived>
643 cmat apply(const Eigen::MatrixBase<Derived>& A,
644  const std::vector<cmat>& Ks,
645  const std::vector<idx>& subsys,
646  const std::vector<idx>& dims)
647 {
648  const cmat& rA = A.derived();
649 
650  // EXCEPTION CHECKS
651 
652  // check zero sizes
654  throw exception::ZeroSize("qpp::apply()");
655 
656  // check square matrix for the A
658  throw exception::MatrixNotSquare("qpp::apply()");
659 
660  // check that dimension is valid
661  if (!internal::check_dims(dims))
662  throw exception::DimsInvalid("qpp::apply()");
663 
664  // check that dims match A matrix
665  if (!internal::check_dims_match_mat(dims, rA))
666  throw exception::DimsMismatchMatrix("qpp::apply()");
667 
668  // check subsys is valid w.r.t. dims
669  if (!internal::check_subsys_match_dims(subsys, dims))
670  throw exception::SubsysMismatchDims("qpp::apply()");
671 
672  std::vector<idx> subsys_dims(subsys.size());
673  for (idx i = 0; i < subsys.size(); ++i)
674  subsys_dims[i] = dims[subsys[i]];
675 
676  // check the Kraus operators
677  if (Ks.size() == 0)
678  throw exception::ZeroSize("qpp::apply()");
679  if (!internal::check_square_mat(Ks[0]))
680  throw exception::MatrixNotSquare("qpp::apply()");
681  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
682  throw exception::MatrixMismatchSubsys("qpp::apply()");
683  for (auto&& it : Ks)
684  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
685  throw exception::DimsNotEqual("qpp::apply()");
686  // END EXCEPTION CHECKS
687 
688  cmat result = cmat::Zero(rA.rows(), rA.rows());
689 
690  for (idx i = 0; i < Ks.size(); ++i)
691  result += apply(rA, Ks[i], subsys, dims);
692 
693  return result;
694 }
695 
706 template<typename Derived>
707 cmat apply(const Eigen::MatrixBase<Derived>& A,
708  const std::vector<cmat>& Ks,
709  const std::vector<idx>& subsys,
710  idx d = 2)
711 {
712  const cmat& rA = A.derived();
713 
714  // EXCEPTION CHECKS
715 
716  // check zero sizes
718  throw exception::ZeroSize("qpp::apply()");
719 
720  // check valid dims
721  if (d < 2)
722  throw exception::DimsInvalid("qpp::apply()");
723  // END EXCEPTION CHECKS
724 
725  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
726  std::vector<idx> dims(N, d); // local dimensions vector
727 
728  return apply(rA, Ks, subsys, dims);
729 }
730 
742 inline cmat kraus2super(const std::vector<cmat>& Ks)
743 {
744  // EXCEPTION CHECKS
745 
746  if (Ks.size() == 0)
747  throw exception::ZeroSize("qpp::kraus2super()");
748  if (!internal::check_nonzero_size(Ks[0]))
749  throw exception::ZeroSize("qpp::kraus2super()");
750  if (!internal::check_square_mat(Ks[0]))
751  throw exception::MatrixNotSquare("qpp::kraus2super()");
752  for (auto&& it : Ks)
753  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
754  throw exception::DimsNotEqual("qpp::kraus2super()");
755  // END EXCEPTION CHECKS
756 
757  idx D = static_cast<idx>(Ks[0].rows());
758 
759  cmat result(D * D, D * D);
760  cmat MN = cmat::Zero(D, D);
761  bra A = bra::Zero(D);
762  ket B = ket::Zero(D);
763  cmat EMN = cmat::Zero(D, D);
764 
765 #ifdef WITH_OPENMP_
766 #pragma omp parallel for collapse(2)
767 #endif // WITH_OPENMP_
768  for (idx m = 0; m < D; ++m)
769  {
770  for (idx n = 0; n < D; ++n)
771  {
772 #ifdef WITH_OPENMP_
773 #pragma omp critical
774 #endif // WITH_OPENMP_
775  {
776  // compute E(|m><n|)
777  MN(m, n) = 1;
778  for (idx i = 0; i < Ks.size(); ++i)
779  EMN += Ks[i] * MN * adjoint(Ks[i]);
780  MN(m, n) = 0;
781 
782  for (idx a = 0; a < D; ++a)
783  {
784  A(a) = 1;
785  for (idx b = 0; b < D; ++b)
786  {
787  // compute result(ab,mn)=<a|E(|m><n)|b>
788  B(b) = 1;
789  result(a * D + b, m * D + n) =
790  static_cast<cmat>(A * EMN * B).value();
791  B(b) = 0;
792  }
793  A(a) = 0;
794  }
795  EMN = cmat::Zero(D, D);
796  }
797  }
798  }
799 
800  return result;
801 }
802 
818 inline cmat kraus2choi(const std::vector<cmat>& Ks)
819 {
820  // EXCEPTION CHECKS
821 
822  if (Ks.size() == 0)
823  throw exception::ZeroSize("qpp::kraus2choi()");
824  if (!internal::check_nonzero_size(Ks[0]))
825  throw exception::ZeroSize("qpp::kraus2choi()");
826  if (!internal::check_square_mat(Ks[0]))
827  throw exception::MatrixNotSquare("qpp::kraus2choi()");
828  for (auto&& it : Ks)
829  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
830  throw exception::DimsNotEqual("qpp::kraus2choi()");
831  // END EXCEPTION CHECKS
832 
833  idx D = static_cast<idx>(Ks[0].rows());
834 
835  // construct the D x D \sum |jj> vector
836  // (un-normalized maximally entangled state)
837  cmat MES = cmat::Zero(D * D, 1);
838  for (idx a = 0; a < D; ++a)
839  MES(a * D + a) = 1;
840 
841  cmat Omega = MES * adjoint(MES);
842 
843  cmat result = cmat::Zero(D * D, D * D);
844 
845 #ifdef WITH_OPENMP_
846 #pragma omp parallel for
847 #endif // WITH_OPENMP_
848  for (idx i = 0; i < Ks.size(); ++i)
849  {
850 #ifdef WITH_OPENMP_
851 #pragma omp critical
852 #endif // WITH_OPENMP_
853  {
854  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
855  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
856  }
857  }
858 
859  return result;
860 }
861 
875 inline std::vector<cmat> choi2kraus(const cmat& A)
876 {
877  // EXCEPTION CHECKS
878 
880  throw exception::ZeroSize("qpp::choi2kraus()");
882  throw exception::MatrixNotSquare("qpp::choi2kraus()");
883  idx D = internal::get_dim_subsys(A.rows(), 2);
884  // check equal dimensions
885  if (D * D != static_cast<idx>(A.rows()))
886  throw exception::DimsInvalid("qpp::choi2kraus()");
887  // END EXCEPTION CHECKS
888 
889  dmat ev = hevals(A);
890  cmat evec = hevects(A);
891  std::vector<cmat> result;
892 
893  for (idx i = 0; i < D * D; ++i)
894  {
895  if (std::abs(ev(i)) > eps)
896  result.push_back(
897  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
898  }
899 
900  return result;
901 }
902 
910 inline cmat choi2super(const cmat& A)
911 {
912  // EXCEPTION CHECKS
913 
915  throw exception::ZeroSize("qpp::choi2super()");
917  throw exception::MatrixNotSquare("qpp::choi2super()");
918  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
919  // check equal dimensions
920  if (D * D != static_cast<idx>(A.rows()))
921  throw exception::DimsInvalid("qpp::choi2super()");
922  // END EXCEPTION CHECKS
923 
924  cmat result(D * D, D * D);
925 
926 #ifdef WITH_OPENMP_
927 #pragma omp parallel for collapse(4)
928 #endif // WITH_OPENMP_
929  for (idx a = 0; a < D; ++a)
930  for (idx b = 0; b < D; ++b)
931  for (idx m = 0; m < D; ++m)
932  for (idx n = 0; n < D; ++n)
933  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
934 
935  return result;
936 }
937 
945 inline cmat super2choi(const cmat& A)
946 {
947  // EXCEPTION CHECKS
948 
950  throw exception::ZeroSize("qpp::super2choi()");
952  throw exception::MatrixNotSquare("qpp::super2choi()");
953  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
954  // check equal dimensions
955  if (D * D != static_cast<idx>(A.rows()))
956  throw exception::DimsInvalid("qpp::super2choi()");
957  // END EXCEPTION CHECKS
958 
959  cmat result(D * D, D * D);
960 
961 #ifdef WITH_OPENMP_
962 #pragma omp parallel for collapse(4)
963 #endif // WITH_OPENMP_
964  for (idx a = 0; a < D; ++a)
965  for (idx b = 0; b < D; ++b)
966  for (idx m = 0; m < D; ++m)
967  for (idx n = 0; n < D; ++n)
968  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
969 
970  return result;
971 }
972 
986 template<typename Derived>
987 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
988  const std::vector<idx>& dims)
989 {
990  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
991  = A.derived();
992 
993  // EXCEPTION CHECKS
994 
995  // check zero-size
997  throw exception::ZeroSize("qpp::ptrace1()");
998 
999  // check that dims is a valid dimension vector
1000  if (!internal::check_dims(dims))
1001  throw exception::DimsInvalid("qpp::ptrace1()");
1002 
1003  // check dims has only 2 elements
1004  if (dims.size() != 2)
1005  throw exception::NotBipartite("qpp::ptrace1()");
1006  // END EXCEPTION CHECKS
1007 
1008  idx DA = dims[0];
1009  idx DB = dims[1];
1010 
1013 
1014  //************ ket ************//
1015  if (internal::check_cvector(rA)) // we have a ket
1016  {
1017  // check that dims match the dimension of A
1018  if (!internal::check_dims_match_cvect(dims, rA))
1019  throw exception::DimsMismatchCvector("qpp::ptrace1()");
1020 
1021  auto worker = [&](idx i, idx j) noexcept -> typename Derived::Scalar
1022  {
1023  typename Derived::Scalar sum = 0;
1024  for (idx m = 0; m < DA; ++m)
1025  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
1026 
1027  return sum;
1028  }; /* end worker */
1029 
1030 #ifdef WITH_OPENMP_
1031 #pragma omp parallel for collapse(2)
1032 #endif // WITH_OPENMP_
1033  for (idx j = 0; j < DB; ++j) // column major order for speed
1034  for (idx i = 0; i < DB; ++i)
1035  result(i, j) = worker(i, j);
1036 
1037  return result;
1038  }
1039  //************ density matrix ************//
1040  else if (internal::check_square_mat(rA)) // we have a density operator
1041  {
1042  // check that dims match the dimension of A
1043  if (!internal::check_dims_match_mat(dims, rA))
1044  throw exception::DimsMismatchMatrix("qpp::ptrace1()");
1045 
1046  auto worker = [&](idx i, idx j) noexcept -> typename Derived::Scalar
1047  {
1048  typename Derived::Scalar sum = 0;
1049  for (idx m = 0; m < DA; ++m)
1050  sum += rA(m * DB + i, m * DB + j);
1051 
1052  return sum;
1053  }; /* end worker */
1054 
1055 #ifdef WITH_OPENMP_
1056 #pragma omp parallel for collapse(2)
1057 #endif // WITH_OPENMP_
1058  for (idx j = 0; j < DB; ++j) // column major order for speed
1059  for (idx i = 0; i < DB; ++i)
1060  result(i, j) = worker(i, j);
1061 
1062  return result;
1063  }
1064  //************ Exception: not ket nor density matrix ************//
1065  else
1066  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1067 }
1068 
1082 template<typename Derived>
1083 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1084  idx d = 2)
1085 {
1086  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1087 
1088  // EXCEPTION CHECKS
1089 
1090  // check zero size
1092  throw exception::ZeroSize("qpp::ptrace1()");
1093 
1094  // check valid dims
1095  if (d == 0)
1096  throw exception::DimsInvalid("qpp::ptrace1()");
1097  // END EXCEPTION CHECKS
1098 
1099  std::vector<idx> dims(2, d); // local dimensions vector
1100 
1101  return ptrace1(A, dims);
1102 }
1103 
1117 template<typename Derived>
1118 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1119  const std::vector<idx>& dims)
1120 {
1121  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1122  = A.derived();
1123 
1124  // EXCEPTION CHECKS
1125 
1126  // check zero-size
1128  throw exception::ZeroSize("qpp::ptrace2()");
1129 
1130  // check that dims is a valid dimension vector
1131  if (!internal::check_dims(dims))
1132  throw exception::DimsInvalid("qpp::ptrace2()");
1133 
1134  // check dims has only 2 elements
1135  if (dims.size() != 2)
1136  throw exception::NotBipartite("qpp::ptrace2()");
1137  // END EXCEPTION CHECKS
1138 
1139  idx DA = dims[0];
1140  idx DB = dims[1];
1141 
1144 
1145  //************ ket ************//
1146  if (internal::check_cvector(rA)) // we have a ket
1147  {
1148  // check that dims match the dimension of A
1149  if (!internal::check_dims_match_cvect(dims, rA))
1150  throw exception::DimsMismatchCvector("qpp::ptrace2()");
1151 
1152  auto worker = [&](idx i, idx j) noexcept -> typename Derived::Scalar
1153  {
1154  typename Derived::Scalar sum = 0;
1155  for (idx m = 0; m < DB; ++m)
1156  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1157 
1158  return sum;
1159  }; /* end worker */
1160 
1161 #ifdef WITH_OPENMP_
1162 #pragma omp parallel for collapse(2)
1163 #endif // WITH_OPENMP_
1164  for (idx j = 0; j < DA; ++j) // column major order for speed
1165  for (idx i = 0; i < DA; ++i)
1166  result(i, j) = worker(i, j);
1167 
1168  return result;
1169  }
1170  //************ density matrix ************//
1171  else if (internal::check_square_mat(rA)) // we have a density operator
1172  {
1173  // check that dims match the dimension of A
1174  if (!internal::check_dims_match_mat(dims, rA))
1175  throw exception::DimsMismatchMatrix("qpp::ptrace2()");
1176 
1177 #ifdef WITH_OPENMP_
1178 #pragma omp parallel for collapse(2)
1179 #endif // WITH_OPENMP_
1180  for (idx j = 0; j < DA; ++j) // column major order for speed
1181  for (idx i = 0; i < DA; ++i)
1182  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1183 
1184  return result;
1185  }
1186  //************ Exception: not ket nor density matrix ************//
1187  else
1188  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1189 }
1190 
1204 template<typename Derived>
1205 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1206  idx d = 2)
1207 {
1208  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1209 
1210  // EXCEPTION CHECKS
1211 
1212  // check zero size
1214  throw exception::ZeroSize("qpp::ptrace2()");
1215 
1216  // check valid dims
1217  if (d == 0)
1218  throw exception::DimsInvalid("qpp::ptrace2()");
1219  // END EXCEPTION CHECKS
1220 
1221  std::vector<idx> dims(2, d); // local dimensions vector
1222 
1223  return ptrace2(A, dims);
1224 }
1225 
1240 template<typename Derived>
1241 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1242  const std::vector<idx>& subsys,
1243  const std::vector<idx>& dims)
1244 {
1245  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1246  = A.derived();
1247 
1248  // EXCEPTION CHECKS
1249 
1250  // check zero-size
1252  throw exception::ZeroSize("qpp::ptrace()");
1253 
1254  // check that dims is a valid dimension vector
1255  if (!internal::check_dims(dims))
1256  throw exception::DimsInvalid("qpp::ptrace()");
1257 
1258  // check that subsys are valid
1259  if (!internal::check_subsys_match_dims(subsys, dims))
1260  throw exception::SubsysMismatchDims("qpp::ptrace()");
1261  // END EXCEPTION CHECKS
1262 
1263  idx D = static_cast<idx>(rA.rows());
1264  idx N = dims.size();
1265  idx Nsubsys = subsys.size();
1266  idx Nsubsys_bar = N - Nsubsys;
1267  idx Dsubsys = 1;
1268  for (idx i = 0; i < Nsubsys; ++i)
1269  Dsubsys *= dims[subsys[i]];
1270  idx Dsubsys_bar = D / Dsubsys;
1271 
1272  idx Cdims[maxn];
1273  idx Csubsys[maxn];
1274  idx Cdimssubsys[maxn];
1275  idx Csubsys_bar[maxn];
1276  idx Cdimssubsys_bar[maxn];
1277 
1278  idx Cmidxcolsubsys_bar[maxn];
1279 
1280  std::vector<idx> subsys_bar = complement(subsys, N);
1281  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1282  std::begin(Csubsys_bar));
1283 
1284  for (idx i = 0; i < N; ++i)
1285  {
1286  Cdims[i] = dims[i];
1287  }
1288  for (idx i = 0; i < Nsubsys; ++i)
1289  {
1290  Csubsys[i] = subsys[i];
1291  Cdimssubsys[i] = dims[subsys[i]];
1292  }
1293  for (idx i = 0; i < Nsubsys_bar; ++i)
1294  {
1295  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1296  }
1297 
1299  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1300 
1301  //************ ket ************//
1302  if (internal::check_cvector(rA)) // we have a ket
1303  {
1304  // check that dims match the dimension of A
1305  if (!internal::check_dims_match_cvect(dims, rA))
1306  throw exception::DimsMismatchCvector("qpp::ptrace()");
1307 
1308  if (subsys.size() == dims.size())
1309  {
1310  result(0, 0) = (adjoint(rA) * rA).value();
1311  return result;
1312  }
1313 
1314  if (subsys.size() == 0)
1315  return rA * adjoint(rA);
1316 
1317  auto worker = [&](idx i) noexcept -> typename Derived::Scalar
1318  {
1319  // use static allocation for speed!
1320 
1321  idx Cmidxrow[maxn];
1322  idx Cmidxcol[maxn];
1323  idx Cmidxrowsubsys_bar[maxn];
1324  idx Cmidxsubsys[maxn];
1325 
1326  /* get the row multi-indexes of the complement */
1327  internal::n2multiidx(i, Nsubsys_bar,
1328  Cdimssubsys_bar, Cmidxrowsubsys_bar);
1329  /* write them in the global row/col multi-indexes */
1330  for (idx k = 0; k < Nsubsys_bar; ++k)
1331  {
1332  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1333  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1334  }
1335  typename Derived::Scalar sm = 0;
1336  for (idx a = 0; a < Dsubsys; ++a)
1337  {
1338  // get the multi-index over which we do the summation
1339  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1340  // write it into the global row/col multi-indexes
1341  for (idx k = 0; k < Nsubsys; ++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  std::conj(rA(internal::multiidx2n(Cmidxcol, N,
1348  Cdims)));
1349  }
1350 
1351  return sm;
1352  }; /* end worker */
1353 
1354  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1355  {
1356  // compute the column multi-indexes of the complement
1357  internal::n2multiidx(j, Nsubsys_bar,
1358  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1359 #ifdef WITH_OPENMP_
1360 #pragma omp parallel for
1361 #endif // WITH_OPENMP_
1362  for (idx i = 0; i < Dsubsys_bar; ++i)
1363  {
1364  result(i, j) = worker(i);
1365  }
1366  }
1367 
1368  return result;
1369  }
1370  //************ density matrix ************//
1371  else if (internal::check_square_mat(rA)) // we have a density operator
1372  {
1373  // check that dims match the dimension of A
1374  if (!internal::check_dims_match_mat(dims, rA))
1375  throw exception::DimsMismatchMatrix("qpp::ptrace()");
1376 
1377  if (subsys.size() == dims.size())
1378  {
1379  result(0, 0) = rA.trace();
1380  return result;
1381  }
1382 
1383  if (subsys.size() == 0)
1384  return rA;
1385 
1386  auto worker = [&](idx i) noexcept -> typename Derived::Scalar
1387  {
1388  // use static allocation for speed!
1389 
1390  idx Cmidxrow[maxn];
1391  idx Cmidxcol[maxn];
1392  idx Cmidxrowsubsys_bar[maxn];
1393  idx Cmidxsubsys[maxn];
1394 
1395  /* get the row/col multi-indexes of the complement */
1396  internal::n2multiidx(i, Nsubsys_bar,
1397  Cdimssubsys_bar, Cmidxrowsubsys_bar);
1398  /* write them in the global row/col multi-indexes */
1399  for (idx k = 0; k < Nsubsys_bar; ++k)
1400  {
1401  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1402  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1403  }
1404  typename Derived::Scalar sm = 0;
1405  for (idx a = 0; a < Dsubsys; ++a)
1406  {
1407  // get the multi-index over which we do the summation
1408  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1409  // write it into the global row/col multi-indexes
1410  for (idx k = 0; k < Nsubsys; ++k)
1411  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1412  = Cmidxsubsys[k];
1413 
1414  // now do the sum
1415  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims),
1416  internal::multiidx2n(Cmidxcol, N, Cdims));
1417  }
1418 
1419  return sm;
1420  }; /* end worker */
1421 
1422  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1423  {
1424  // compute the column multi-indexes of the complement
1425  internal::n2multiidx(j, Nsubsys_bar,
1426  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1427 #ifdef WITH_OPENMP_
1428 #pragma omp parallel for
1429 #endif // WITH_OPENMP_
1430  for (idx i = 0; i < Dsubsys_bar; ++i)
1431  {
1432  result(i, j) = worker(i);
1433  }
1434  }
1435 
1436  return result;
1437  }
1438  //************ Exception: not ket nor density matrix ************//
1439  else
1440  throw exception::MatrixNotSquareNorCvector("qpp::ptrace()");
1441 }
1442 
1457 template<typename Derived>
1458 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1459  const std::vector<idx>& subsys,
1460  idx d = 2)
1461 {
1462  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1463 
1464  // EXCEPTION CHECKS
1465 
1466  // check zero size
1468  throw exception::ZeroSize("qpp::ptrace()");
1469 
1470  // check valid dims
1471  if (d < 2)
1472  throw exception::DimsInvalid("qpp::ptrace()");
1473  // END EXCEPTION CHECKS
1474 
1475  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1476  std::vector<idx> dims(N, d); // local dimensions vector
1477 
1478  return ptrace(rA, subsys, dims);
1479 }
1480 
1494 template<typename Derived>
1496  const Eigen::MatrixBase<Derived>& A,
1497  const std::vector<idx>& subsys,
1498  const std::vector<idx>& dims)
1499 {
1500  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1501 
1502  // EXCEPTION CHECKS
1503 
1504  // check zero-size
1506  throw exception::ZeroSize("qpp::ptranspose()");
1507 
1508  // check that dims is a valid dimension vector
1509  if (!internal::check_dims(dims))
1510  throw exception::DimsInvalid("qpp::ptranspose()");
1511 
1512  // check that subsys are valid
1513  if (!internal::check_subsys_match_dims(subsys, dims))
1514  throw exception::SubsysMismatchDims("qpp::ptranspose()");
1515  // END EXCEPTION CHECKS
1516 
1517  idx D = static_cast<idx>(rA.rows());
1518  idx N = dims.size();
1519  idx Nsubsys = subsys.size();
1520  idx Cdims[maxn];
1521  idx Cmidxcol[maxn];
1522  idx Csubsys[maxn];
1523 
1524  // copy dims in Cdims and subsys in Csubsys
1525  for (idx i = 0; i < N; ++i)
1526  Cdims[i] = dims[i];
1527  for (idx i = 0; i < Nsubsys; ++i)
1528  Csubsys[i] = subsys[i];
1529 
1530  dyn_mat<typename Derived::Scalar> result(D, D);
1531 
1532  //************ ket ************//
1533  if (internal::check_cvector(rA)) // we have a ket
1534  {
1535  // check that dims match the dimension of A
1536  if (!internal::check_dims_match_cvect(dims, rA))
1537  throw exception::DimsMismatchCvector("qpp::ptranspose()");
1538 
1539  if (subsys.size() == dims.size())
1540  return (rA * adjoint(rA)).transpose();
1541 
1542  if (subsys.size() == 0)
1543  return rA * adjoint(rA);
1544 
1545  auto worker = [&](idx i) noexcept -> typename Derived::Scalar
1546  {
1547  // use static allocation for speed!
1548  idx midxcoltmp[maxn];
1549  idx midxrow[maxn];
1550 
1551  for (idx k = 0; k < N; ++k)
1552  midxcoltmp[k] = Cmidxcol[k];
1553 
1554  /* compute the row multi-index */
1555  internal::n2multiidx(i, N, Cdims, midxrow);
1556 
1557  for (idx k = 0; k < Nsubsys; ++k)
1558  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1559 
1560  /* writes the result */
1561  return rA(internal::multiidx2n(midxrow, N, Cdims)) *
1562  std::conj(rA(internal::multiidx2n(midxcoltmp, N,
1563  Cdims)));
1564  }; /* end worker */
1565 
1566  for (idx j = 0; j < D; ++j)
1567  {
1568  // compute the column multi-index
1569  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1570 
1571 #ifdef WITH_OPENMP_
1572 #pragma omp parallel for
1573 #endif // WITH_OPENMP_
1574  for (idx i = 0; i < D; ++i)
1575  result(i, j) = worker(i);
1576  }
1577 
1578  return result;
1579  }
1580  //************ density matrix ************//
1581  else if (internal::check_square_mat(rA)) // we have a density operator
1582  {
1583  // check that dims match the dimension of A
1584  if (!internal::check_dims_match_mat(dims, rA))
1585  throw exception::DimsMismatchMatrix("qpp::ptranspose()");
1586 
1587  if (subsys.size() == dims.size())
1588  return rA.transpose();
1589 
1590  if (subsys.size() == 0)
1591  return rA;
1592 
1593  auto worker = [&](idx i) noexcept -> typename Derived::Scalar
1594  {
1595  // use static allocation for speed!
1596  idx midxcoltmp[maxn];
1597  idx midxrow[maxn];
1598 
1599  for (idx k = 0; k < N; ++k)
1600  midxcoltmp[k] = Cmidxcol[k];
1601 
1602  /* compute the row multi-index */
1603  internal::n2multiidx(i, N, Cdims, midxrow);
1604 
1605  for (idx k = 0; k < Nsubsys; ++k)
1606  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1607 
1608  /* writes the result */
1609  return rA(internal::multiidx2n(midxrow, N, Cdims),
1610  internal::multiidx2n(midxcoltmp, N, Cdims));
1611  }; /* end worker */
1612 
1613  for (idx j = 0; j < D; ++j)
1614  {
1615  // compute the column multi-index
1616  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1617 
1618 #ifdef WITH_OPENMP_
1619 #pragma omp parallel for
1620 #endif // WITH_OPENMP_
1621  for (idx i = 0; i < D; ++i)
1622  result(i, j) = worker(i);
1623  }
1624 
1625  return result;
1626  }
1627  //************ Exception: not ket nor density matrix ************//
1628  else
1629  throw exception::MatrixNotSquareNorCvector("qpp::ptranspose()");
1630 }
1631 
1645 template<typename Derived>
1647  const Eigen::MatrixBase<Derived>& A,
1648  const std::vector<idx>& subsys,
1649  idx d = 2)
1650 {
1651  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1652 
1653  // EXCEPTION CHECKS
1654 
1655  // check zero size
1657  throw exception::ZeroSize("qpp::ptranspose()");
1658 
1659  // check valid dims
1660  if (d < 2)
1661  throw exception::DimsInvalid("qpp::ptranspose()");
1662  // END EXCEPTION CHECKS
1663 
1664  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1665  std::vector<idx> dims(N, d); // local dimensions vector
1666 
1667  return ptranspose(rA, subsys, dims);
1668 }
1669 
1682 template<typename Derived>
1684  const Eigen::MatrixBase<Derived>& A,
1685  const std::vector<idx>& perm,
1686  const std::vector<idx>& dims)
1687 {
1688  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1689 
1690  // EXCEPTION CHECKS
1691 
1692  // check zero-size
1694  throw exception::ZeroSize("qpp::syspermute()");
1695 
1696  // check that dims is a valid dimension vector
1697  if (!internal::check_dims(dims))
1698  throw exception::DimsInvalid("qpp::syspermute()");
1699 
1700  // check that we have a valid permutation
1701  if (!internal::check_perm(perm))
1702  throw exception::PermInvalid("qpp::syspermute()");
1703 
1704  // check that permutation match dimensions
1705  if (perm.size() != dims.size())
1706  throw exception::PermMismatchDims("qpp::syspermute()");
1707  // END EXCEPTION CHECKS
1708 
1709  idx D = static_cast<idx>(rA.rows());
1710  idx N = dims.size();
1711 
1713 
1714  //************ ket ************//
1715  if (internal::check_cvector(rA)) // we have a column vector
1716  {
1717  idx Cdims[maxn];
1718  idx Cperm[maxn];
1719 
1720  // check that dims match the dimension of rA
1721  if (!internal::check_dims_match_cvect(dims, rA))
1722  throw exception::DimsMismatchCvector("qpp::syspermute()");
1723 
1724  // copy dims in Cdims and perm in Cperm
1725  for (idx i = 0; i < N; ++i)
1726  {
1727  Cdims[i] = dims[i];
1728  Cperm[i] = perm[i];
1729  }
1730  result.resize(D, 1);
1731 
1732  auto worker = [&Cdims, &Cperm, N](idx i) noexcept -> idx
1733  {
1734  // use static allocation for speed,
1735  // double the size for matrices reshaped as vectors
1736  idx midx[maxn];
1737  idx midxtmp[maxn];
1738  idx permdims[maxn];
1739 
1740  /* compute the multi-index */
1741  internal::n2multiidx(i, N, Cdims, midx);
1742 
1743  for (idx k = 0; k < N; ++k)
1744  {
1745  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1746  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1747  }
1748 
1749  return internal::multiidx2n(midxtmp, N, permdims);
1750  }; /* end worker */
1751 
1752 #ifdef WITH_OPENMP_
1753 #pragma omp parallel for
1754 #endif // WITH_OPENMP_
1755  for (idx i = 0; i < D; ++i)
1756  result(worker(i)) = rA(i);
1757 
1758  return result;
1759  }
1760  //************ density matrix ************//
1761  else if (internal::check_square_mat(rA)) // we have a density operator
1762  {
1763  idx Cdims[2 * maxn];
1764  idx Cperm[2 * maxn];
1765 
1766  // check that dims match the dimension of rA
1767  if (!internal::check_dims_match_mat(dims, rA))
1768  throw exception::DimsMismatchMatrix("qpp::syspermute()");
1769 
1770  // copy dims in Cdims and perm in Cperm
1771  for (idx i = 0; i < N; ++i)
1772  {
1773  Cdims[i] = dims[i];
1774  Cdims[i + N] = dims[i];
1775  Cperm[i] = perm[i];
1776  Cperm[i + N] = perm[i] + N;
1777  }
1778  result.resize(D * D, 1);
1779  // map A to a column vector
1781  Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1782  const_cast<typename Derived::Scalar*>(rA.data()),
1783  D * D, 1);
1784 
1785  auto worker = [&Cdims, &Cperm, N](idx i) noexcept -> idx
1786  {
1787  // use static allocation for speed,
1788  // double the size for matrices reshaped as vectors
1789  idx midx[2 * maxn];
1790  idx midxtmp[2 * maxn];
1791  idx permdims[2 * maxn];
1792 
1793  /* compute the multi-index */
1794  internal::n2multiidx(i, 2 * N, Cdims, midx);
1795 
1796  for (idx k = 0; k < 2 * N; ++k)
1797  {
1798  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1799  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1800  }
1801 
1802  return internal::multiidx2n(midxtmp, 2 * N, permdims);
1803  }; /* end worker */
1804 
1805 #ifdef WITH_OPENMP_
1806 #pragma omp parallel for
1807 #endif // WITH_OPENMP_
1808  for (idx i = 0; i < D * D; ++i)
1809  result(worker(i)) = rA(i);
1810 
1811  return reshape(result, D, D);
1812  }
1813  //************ Exception: not ket nor density matrix ************//
1814  else
1815  throw exception::MatrixNotSquareNorCvector("qpp::syspermute()");
1816 }
1817 
1830 template<typename Derived>
1832  const Eigen::MatrixBase<Derived>& A,
1833  const std::vector<idx>& perm,
1834  idx d = 2)
1835 {
1836  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1837 
1838  // EXCEPTION CHECKS
1839 
1840  // check zero size
1842  throw exception::ZeroSize("qpp::syspermute()");
1843 
1844  // check valid dims
1845  if (d < 2)
1846  throw exception::DimsInvalid("qpp::syspermute()");
1847  // END EXCEPTION CHECKS
1848 
1849  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1850  std::vector<idx> dims(N, d); // local dimensions vector
1851 
1852  return syspermute(rA, perm, dims);
1853 }
1854 
1855 } /* namespace qpp */
1856 
1857 #endif /* OPERATIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
Dimensions not equal exception.
Definition: exception.h:304
Dimension(s) mismatch matrix size exception.
Definition: exception.h:322
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:71
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:221
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:122
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:875
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:65
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:64
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:55
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:1683
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:77
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:158
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:70
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:40
Subsystems mismatch dimensions exception.
Definition: exception.h:395
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:50
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:473
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:278
idx get_num_subsys(idx sz, idx d)
Definition: util.h:391
Matrix is not square nor column vector exception.
Definition: exception.h:219
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:387
Invalid dimension(s) exception.
Definition: exception.h:287
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:910
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:945
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1811
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:174
Not bi-partite exception.
Definition: exception.h:532
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:987
Invalid permutation exception.
Definition: exception.h:412
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:81
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:899
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1118
Matrix mismatch subsystems exception.
Definition: exception.h:270
Dimension(s) mismatch column vector size exception.
Definition: exception.h:340
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:41
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:777
Type mismatch exception.
Definition: exception.h:585
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:404
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:413
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:818
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:1241
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:742
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:50
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:1495
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:142
Permutation mismatch dimensions exception.
Definition: exception.h:430
std::size_t idx
Non-negative integer index.
Definition: types.h:35
Matrix is not square exception.
Definition: exception.h:151
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:205
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1139
Object has zero size exception.
Definition: exception.h:134