Quantum++  v1.0.0-beta1
C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 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 // silence g++4.8 bogus warning -Wunused-but-set-variable in lambda functions
31 #if (__GNUC__ && !__clang__)
32 #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
33 #endif
34 
35 namespace qpp
36 {
53 template<typename Derived1, typename Derived2>
55  const Eigen::MatrixBase<Derived1>& state,
56  const Eigen::MatrixBase<Derived2>& A,
57  const std::vector<idx>& ctrl,
58  const std::vector<idx>& subsys,
59  const std::vector<idx>& dims)
60 {
61  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
62  = state.derived();
63  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
64 
65  // EXCEPTION CHECKS
66 
67  // check types
68  if (!std::is_same<typename Derived1::Scalar,
69  typename Derived2::Scalar>::value)
70  throw Exception("qpp::applyCTRL()", Exception::Type::TYPE_MISMATCH);
71 
72  // check zero sizes
74  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
75 
76  // check zero sizes
77  if (!internal::check_nonzero_size(rstate))
78  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
79 
80  // check square matrix for the gate
82  throw Exception("qpp::applyCTRL()",
84 
85  // check that all control subsystems have the same dimension
86  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
87  for (idx i = 1; i < ctrl.size(); ++i)
88  if (dims[ctrl[i]] != d)
89  throw Exception("qpp::applyCTRL()",
91 
92  // check that dimension is valid
93  if (!internal::check_dims(dims))
94  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
95 
96  // check subsys is valid w.r.t. dims
97  if (!internal::check_subsys_match_dims(subsys, dims))
98  throw Exception("qpp::applyCTRL()",
100 
101  // check that gate matches the dimensions of the subsys
102  std::vector<idx> subsys_dims(subsys.size());
103  for (idx i = 0; i < subsys.size(); ++i)
104  subsys_dims[i] = dims[subsys[i]];
105  if (!internal::check_dims_match_mat(subsys_dims, rA))
106  throw Exception("qpp::applyCTRL()",
108 
109  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
110  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
111  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
112 
113  // check that ctrl + gate subsystem is valid
114  // with respect to local dimensions
115  if (!internal::check_subsys_match_dims(ctrlgate, dims))
116  throw Exception("qpp::applyCTRL()",
118  // END EXCEPTION CHECKS
119 
120  // construct the table of A^i and (A^dagger)^i
121  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
122  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
123  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i)
124  {
125  Ai.push_back(powm(rA, i));
126  Aidagger.push_back(powm(adjoint(rA), i));
127  }
128 
129  idx D = static_cast<idx>(rstate.rows()); // total dimension
130  idx N = dims.size(); // total number of subsystems
131  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
132  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
133  idx subsyssize = subsys.size(); // number of subsystems of the target
134  // dimension of ctrl subsystem
135  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
136  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
137 
138  idx Cdims[maxn]; // local dimensions
139  idx CdimsA[maxn]; // local dimensions
140  idx CdimsCTRL[maxn]; // local dimensions
141  idx CdimsCTRLA_bar[maxn]; // local dimensions
142 
143  // compute the complementary subsystem of ctrlgate w.r.t. dims
144  std::vector<idx> ctrlgate_bar = complement(ctrlgate, N);
145  // number of subsystems that are complementary to the ctrl+gate
146  idx ctrlgate_barsize = ctrlgate_bar.size();
147 
148  idx DCTRLA_bar = 1; // dimension of the rest
149  for (idx i = 0; i < ctrlgate_barsize; ++i)
150  DCTRLA_bar *= dims[ctrlgate_bar[i]];
151 
152  for (idx k = 0; k < N; ++k)
153  Cdims[k] = dims[k];
154  for (idx k = 0; k < subsyssize; ++k)
155  CdimsA[k] = dims[subsys[k]];
156  for (idx k = 0; k < ctrlsize; ++k)
157  CdimsCTRL[k] = d;
158  for (idx k = 0; k < ctrlgate_barsize; ++k)
159  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
160 
161  // worker, computes the coefficient and the index for the ket case
162  // used in #pragma omp parallel for collapse
163  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
164  -> std::pair<typename Derived1::Scalar, idx>
165  {
166  idx indx = 0;
167  typename Derived1::Scalar coeff = 0;
168 
169  idx Cmidx[maxn]; // the total multi-index
170  idx CmidxA[maxn]; // the gate part multi-index
171  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
172 
173  // compute the index
174 
175  // set the CTRL part
176  for (idx k = 0; k < ctrlsize; ++k)
177  {
178  Cmidx[ctrl[k]] = i_;
179  }
180 
181  // set the rest
182  internal::n2multiidx(r_, N - ctrlgatesize,
183  CdimsCTRLA_bar, CmidxCTRLA_bar);
184  for (idx k = 0; k < N - ctrlgatesize; ++k)
185  {
186  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
187  }
188 
189  // set the A part
190  internal::n2multiidx(m_, subsyssize, CdimsA, CmidxA);
191  for (idx k = 0; k < subsyssize; ++k)
192  {
193  Cmidx[subsys[k]] = CmidxA[k];
194  }
195 
196  // we now got the total index
197  indx = internal::multiidx2n(Cmidx, N, Cdims);
198 
199  // compute the coefficient
200  for (idx n_ = 0; n_ < DA; ++n_)
201  {
202  internal::n2multiidx(n_, subsyssize, CdimsA, CmidxA);
203  for (idx k = 0; k < subsyssize; ++k)
204  {
205  Cmidx[subsys[k]] = CmidxA[k];
206  }
207  coeff += Ai[i_](m_, n_) *
208  rstate(internal::multiidx2n(Cmidx, N, Cdims));
209  }
210 
211  return std::make_pair(coeff, indx);
212  }; /* end coeff_idx_ket */
213 
214  // worker, computes the coefficient and the index
215  // for the density matrix case
216  // used in #pragma omp parallel for collapse
217  auto coeff_idx_rho = [&](idx i1_, idx m1_,
218  idx r1_, idx i2_, idx m2_, idx r2_) noexcept
219  -> std::tuple<typename Derived1::Scalar, idx, idx>
220  {
221  idx idxrow = 0;
222  idx idxcol = 0;
223  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
224 
225  idx Cmidxrow[maxn]; // the total row multi-index
226  idx Cmidxcol[maxn]; // the total col multi-index
227  idx CmidxArow[maxn]; // the gate part row multi-index
228  idx CmidxAcol[maxn]; // the gate part col multi-index
229  idx CmidxCTRLrow[maxn]; // the control row multi-index
230  idx CmidxCTRLcol[maxn]; // the control col multi-index
231  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
232  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
233 
234  // compute the ket/bra indexes
235 
236  // set the CTRL part
237  internal::n2multiidx(i1_, ctrlsize,
238  CdimsCTRL, CmidxCTRLrow);
239  internal::n2multiidx(i2_, ctrlsize,
240  CdimsCTRL, CmidxCTRLcol);
241 
242  for (idx k = 0; k < ctrlsize; ++k)
243  {
244  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
245  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
246  }
247 
248  // set the rest
249  internal::n2multiidx(r1_, N - ctrlgatesize,
250  CdimsCTRLA_bar, CmidxCTRLA_barrow);
251  internal::n2multiidx(r2_, N - ctrlgatesize,
252  CdimsCTRLA_bar, CmidxCTRLA_barcol);
253  for (idx k = 0; k < N - ctrlgatesize; ++k)
254  {
255  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
256  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
257  }
258 
259  // set the A part
260  internal::n2multiidx(m1_, subsyssize, CdimsA, CmidxArow);
261  internal::n2multiidx(m2_, subsyssize, CdimsA, CmidxAcol);
262  for (idx k = 0; k < subsys.size(); ++k)
263  {
264  Cmidxrow[subsys[k]] = CmidxArow[k];
265  Cmidxcol[subsys[k]] = CmidxAcol[k];
266  }
267 
268  // we now got the total row/col indexes
269  idxrow = internal::multiidx2n(Cmidxrow, N, Cdims);
270  idxcol = internal::multiidx2n(Cmidxcol, N, Cdims);
271 
272  // check whether all CTRL row and col multi indexes are equal
273  bool all_ctrl_rows_equal = true;
274  bool all_ctrl_cols_equal = true;
275 
276  idx first_ctrl_row, first_ctrl_col;
277  if (ctrlsize > 0)
278  {
279  first_ctrl_row = CmidxCTRLrow[0];
280  first_ctrl_col = CmidxCTRLcol[0];
281  } else
282  {
283  first_ctrl_row = first_ctrl_col = 1;
284  }
285 
286  for (idx k = 1; k < ctrlsize; ++k)
287  {
288  if (CmidxCTRLrow[k] != first_ctrl_row)
289  {
290  all_ctrl_rows_equal = false;
291  break;
292  }
293  }
294  for (idx k = 1; k < ctrlsize; ++k)
295  {
296  if (CmidxCTRLcol[k] != first_ctrl_col)
297  {
298  all_ctrl_cols_equal = false;
299  break;
300  }
301  }
302 
303  // at least one control activated, compute the coefficient
304  for (idx n1_ = 0; n1_ < DA; ++n1_)
305  {
306  internal::n2multiidx(n1_, subsyssize, CdimsA, CmidxArow);
307  for (idx k = 0; k < subsyssize; ++k)
308  {
309  Cmidxrow[subsys[k]] = CmidxArow[k];
310  }
311  idx idxrowtmp = internal::multiidx2n(Cmidxrow, N, Cdims);
312 
313  if (all_ctrl_rows_equal)
314  {
315  lhs = Ai[first_ctrl_row](m1_, n1_);
316  } else
317  {
318  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
319  }
320 
321  for (idx n2_ = 0; n2_ < DA; ++n2_)
322  {
323  internal::n2multiidx(n2_, subsyssize, CdimsA, CmidxAcol);
324  for (idx k = 0; k < subsyssize; ++k)
325  {
326  Cmidxcol[subsys[k]] = CmidxAcol[k];
327  }
328 
329  if (all_ctrl_cols_equal)
330  {
331  rhs = Aidagger[first_ctrl_col](n2_, m2_);
332  } else
333  {
334  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
335  }
336 
337  idx idxcoltmp = internal::multiidx2n(Cmidxcol, N, Cdims);
338 
339  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
340  }
341  }
342 
343  return std::make_tuple(coeff, idxrow, idxcol);
344  }; /* end coeff_idx_rho */
345 
346  //************ ket ************//
347  if (internal::check_cvector(rstate)) // we have a ket
348  {
349  // check that dims match state vector
350  if (!internal::check_dims_match_cvect(dims, rstate))
351  throw Exception("qpp::applyCTRL()",
353  if (D == 1)
354  return rstate;
355 
356  dyn_mat<typename Derived1::Scalar> result = rstate;
357 
358 #ifdef WITH_OPENMP_
359 #pragma omp parallel for collapse(2)
360 #endif // WITH_OPENMP_
361  for (idx m = 0; m < DA; ++m)
362  for (idx r = 0; r < DCTRLA_bar; ++r)
363  {
364  if (ctrlsize == 0) // no control
365  {
366  result(coeff_idx_ket(1, m, r).second) =
367  coeff_idx_ket(1, m, r).first;
368  } else
369  for (idx i = 0; i < d; ++i)
370  {
371  result(coeff_idx_ket(i, m, r).second) =
372  coeff_idx_ket(i, m, r).first;
373  }
374  }
375 
376  return result;
377  }
378  //************ density matrix ************//
379  else if (internal::check_square_mat(rstate)) // we have a density operator
380  {
381  // check that dims match state matrix
382  if (!internal::check_dims_match_mat(dims, rstate))
383  throw Exception("qpp::applyCTRL()",
385 
386  if (D == 1)
387  return rstate;
388 
389  dyn_mat<typename Derived1::Scalar> result = rstate;
390 
391 #ifdef WITH_OPENMP_
392 #pragma omp parallel for collapse(4)
393 #endif // WITH_OPENMP_
394  for (idx m1 = 0; m1 < DA; ++m1)
395  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
396  for (idx m2 = 0; m2 < DA; ++m2)
397  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
398  if (ctrlsize == 0) // no control
399  {
400  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
401  1, m2, r2);
402  result(std::get<1>(coeff_idxes),
403  std::get<2>(coeff_idxes)) =
404  std::get<0>(coeff_idxes);
405  } else
406  {
407  for (idx i1 = 0; i1 < Dctrl; ++i1)
408  for (idx i2 = 0; i2 < Dctrl; ++i2)
409  {
410  auto coeff_idxes = coeff_idx_rho(
411  i1, m1, r1,
412  i2, m2, r2);
413  result(std::get<1>(coeff_idxes),
414  std::get<2>(coeff_idxes)) =
415  std::get<0>(coeff_idxes);
416  }
417  }
418 
419  return result;
420  }
421  //************ Exception: not ket nor density matrix ************//
422  else
423  throw Exception("qpp::applyCTRL()",
425 }
426 
442 template<typename Derived1, typename Derived2>
444  const Eigen::MatrixBase<Derived1>& state,
445  const Eigen::MatrixBase<Derived2>& A,
446  const std::vector<idx>& ctrl,
447  const std::vector<idx>& subsys,
448  idx d = 2)
449 {
450  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
451  = state.derived();
452  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
453 
454  // EXCEPTION CHECKS
455 
456  // check zero size
457  if (!internal::check_nonzero_size(rstate))
458  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
459 
460  // check valid dims
461  if (d == 0)
462  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
463  // END EXCEPTION CHECKS
464 
465  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
466  std::vector<idx> dims(N, d); // local dimensions vector
467 
468  return applyCTRL(rstate, rA, ctrl, subsys, dims);
469 }
470 
484 template<typename Derived1, typename Derived2>
486  const Eigen::MatrixBase<Derived1>& state,
487  const Eigen::MatrixBase<Derived2>& A,
488  const std::vector<idx>& subsys,
489  const std::vector<idx>& dims)
490 {
491  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
492  = state.derived();
493  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
494 
495  // EXCEPTION CHECKS
496 
497  // check types
498  if (!std::is_same<typename Derived1::Scalar,
499  typename Derived2::Scalar>::value)
500  throw Exception("qpp::apply()", Exception::Type::TYPE_MISMATCH);
501 
502  // check zero sizes
504  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
505 
506  // check zero sizes
507  if (!internal::check_nonzero_size(rstate))
508  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
509 
510  // check square matrix for the gate
512  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
513 
514  // check that dimension is valid
515  if (!internal::check_dims(dims))
516  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
517 
518  // check subsys is valid w.r.t. dims
519  if (!internal::check_subsys_match_dims(subsys, dims))
520  throw Exception("qpp::apply()", Exception::Type::SUBSYS_MISMATCH_DIMS);
521 
522  // check that gate matches the dimensions of the subsys
523  std::vector<idx> subsys_dims(subsys.size());
524  for (idx i = 0; i < subsys.size(); ++i)
525  subsys_dims[i] = dims[subsys[i]];
526  if (!internal::check_dims_match_mat(subsys_dims, rA))
527  throw Exception("qpp::apply()",
529  // END EXCEPTION CHECKS
530 
531  //************ ket ************//
532  if (internal::check_cvector(rstate)) // we have a ket
533  {
534  // check that dims match state vector
535  if (!internal::check_dims_match_cvect(dims, rstate))
536  throw Exception("qpp::apply()",
538 
539  return applyCTRL(rstate, rA, {}, subsys, dims);
540  }
541  //************ density matrix ************//
542  else if (internal::check_square_mat(rstate)) // we have a density operator
543  {
544 
545  // check that dims match state matrix
546  if (!internal::check_dims_match_mat(dims, rstate))
547  throw Exception("qpp::apply()",
549 
550  return applyCTRL(rstate, rA, {}, subsys, dims);
551  }
552  //************ Exception: not ket nor density matrix ************//
553  else
554  throw Exception("qpp::apply()",
556 }
557 
571 template<typename Derived1, typename Derived2>
573  const Eigen::MatrixBase<Derived1>& state,
574  const Eigen::MatrixBase<Derived2>& A,
575  const std::vector<idx>& subsys,
576  idx d = 2)
577 {
578  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
579  = state.derived();
580  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
581 
582  // EXCEPTION CHECKS
583 
584  // check zero size
585  if (!internal::check_nonzero_size(rstate))
586  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
587 
588  // check valid dims
589  if (d == 0)
590  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
591  // END EXCEPTION CHECKS
592 
593  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
594  std::vector<idx> dims(N, d); // local dimensions vector
595 
596  return apply(rstate, rA, subsys, dims);
597 }
598 
607 template<typename Derived>
608 cmat apply(const Eigen::MatrixBase<Derived>& rho,
609  const std::vector<cmat>& Ks)
610 {
611  const cmat& rrho = rho.derived();
612 
613  // EXCEPTION CHECKS
614 
615  if (!internal::check_nonzero_size(rrho))
616  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
617  if (!internal::check_square_mat(rrho))
618  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
619  if (Ks.size() == 0)
620  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
621  if (!internal::check_square_mat(Ks[0]))
622  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
623  if (Ks[0].rows() != rrho.rows())
624  throw Exception("qpp::apply()",
626  for (auto&& it : Ks)
627  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
628  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
629  // END EXCEPTION CHECKS
630 
631  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
632 
633 #ifdef WITH_OPENMP_
634 #pragma omp parallel for
635 #endif // WITH_OPENMP_
636  for (idx i = 0; i < Ks.size(); ++i)
637  {
638 #ifdef WITH_OPENMP_
639 #pragma omp critical
640 #endif // WITH_OPENMP_
641  {
642  result += Ks[i] * rrho * adjoint(Ks[i]);
643  }
644  }
645 
646  return result;
647 }
648 
659 template<typename Derived>
660 cmat apply(const Eigen::MatrixBase<Derived>& rho,
661  const std::vector<cmat>& Ks,
662  const std::vector<idx>& subsys,
663  const std::vector<idx>& dims)
664 {
665  const cmat& rrho = rho.derived();
666 
667  // EXCEPTION CHECKS
668 
669  // check zero sizes
670  if (!internal::check_nonzero_size(rrho))
671  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
672 
673  // check square matrix for the rho
674  if (!internal::check_square_mat(rrho))
675  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
676 
677  // check that dimension is valid
678  if (!internal::check_dims(dims))
679  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
680 
681  // check that dims match rho matrix
682  if (!internal::check_dims_match_mat(dims, rrho))
683  throw Exception("qpp::apply()",
685 
686  // check subsys is valid w.r.t. dims
687  if (!internal::check_subsys_match_dims(subsys, dims))
688  throw Exception("qpp::apply()",
690 
691  std::vector<idx> subsys_dims(subsys.size());
692  for (idx i = 0; i < subsys.size(); ++i)
693  subsys_dims[i] = dims[subsys[i]];
694 
695  // check the Kraus operators
696  if (Ks.size() == 0)
697  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
698  if (!internal::check_square_mat(Ks[0]))
699  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
700  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
701  throw Exception("qpp::apply()",
703  for (auto&& it : Ks)
704  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
705  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
706  // END EXCEPTION CHECKS
707 
708  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
709 
710  for (idx i = 0; i < Ks.size(); ++i)
711  result += apply(rrho, Ks[i], subsys, dims);
712 
713  return result;
714 }
715 
726 template<typename Derived>
727 cmat apply(const Eigen::MatrixBase<Derived>& rho,
728  const std::vector<cmat>& Ks,
729  const std::vector<idx>& subsys,
730  idx d = 2)
731 {
732  const cmat& rrho = rho.derived();
733 
734  // EXCEPTION CHECKS
735 
736  // check zero sizes
737  if (!internal::check_nonzero_size(rrho))
738  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
739 
740  // check valid dims
741  if (d == 0)
742  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
743  // END EXCEPTION CHECKS
744 
745  idx N = internal::get_num_subsys(static_cast<idx>(rrho.rows()), d);
746  std::vector<idx> dims(N, d); // local dimensions vector
747 
748  return apply(rrho, Ks, subsys, dims);
749 }
750 
762 inline cmat kraus2super(const std::vector<cmat>& Ks)
763 {
764  // EXCEPTION CHECKS
765 
766  if (Ks.size() == 0)
767  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
768  if (!internal::check_nonzero_size(Ks[0]))
769  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
770  if (!internal::check_square_mat(Ks[0]))
771  throw Exception("qpp::kraus2super()",
773  for (auto&& it : Ks)
774  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
775  throw Exception("qpp::kraus2super()",
777  // END EXCEPTION CHECKS
778 
779  idx D = static_cast<idx>(Ks[0].rows());
780 
781  cmat result(D * D, D * D);
782  cmat MN = cmat::Zero(D, D);
783  bra A = bra::Zero(D);
784  ket B = ket::Zero(D);
785  cmat EMN = cmat::Zero(D, D);
786 
787 #ifdef WITH_OPENMP_
788 #pragma omp parallel for collapse(2)
789 #endif // WITH_OPENMP_
790  for (idx m = 0; m < D; ++m)
791  {
792  for (idx n = 0; n < D; ++n)
793  {
794 #ifdef WITH_OPENMP_
795 #pragma omp critical
796 #endif // WITH_OPENMP_
797  {
798  // compute E(|m><n|)
799  MN(m, n) = 1;
800  for (idx i = 0; i < Ks.size(); ++i)
801  EMN += Ks[i] * MN * adjoint(Ks[i]);
802  MN(m, n) = 0;
803 
804  for (idx a = 0; a < D; ++a)
805  {
806  A(a) = 1;
807  for (idx b = 0; b < D; ++b)
808  {
809  // compute result(ab,mn)=<a|E(|m><n)|b>
810  B(b) = 1;
811  result(a * D + b, m * D + n) =
812  static_cast<cmat>(A * EMN * B).value();
813  B(b) = 0;
814  }
815  A(a) = 0;
816  }
817  EMN = cmat::Zero(D, D);
818  }
819  }
820  }
821 
822  return result;
823 }
824 
840 inline cmat kraus2choi(const std::vector<cmat>& Ks)
841 {
842  // EXCEPTION CHECKS
843 
844  if (Ks.size() == 0)
845  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
846  if (!internal::check_nonzero_size(Ks[0]))
847  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
848  if (!internal::check_square_mat(Ks[0]))
849  throw Exception("qpp::kraus2choi()",
851  for (auto&& it : Ks)
852  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
853  throw Exception("qpp::kraus2choi()",
855  // END EXCEPTION CHECKS
856 
857  idx D = static_cast<idx>(Ks[0].rows());
858 
859  // construct the D x D \sum |jj> vector
860  // (un-normalized maximally entangled state)
861  cmat MES = cmat::Zero(D * D, 1);
862  for (idx a = 0; a < D; ++a)
863  MES(a * D + a) = 1;
864 
865  cmat Omega = MES * adjoint(MES);
866 
867  cmat result = cmat::Zero(D * D, D * D);
868 
869 #ifdef WITH_OPENMP_
870 #pragma omp parallel for
871 #endif // WITH_OPENMP_
872  for (idx i = 0; i < Ks.size(); ++i)
873  {
874 #ifdef WITH_OPENMP_
875 #pragma omp critical
876 #endif // WITH_OPENMP_
877  {
878  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
879  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
880  }
881  }
882 
883  return result;
884 }
885 
899 inline std::vector<cmat> choi2kraus(const cmat& A)
900 {
901  // EXCEPTION CHECKS
902 
904  throw Exception("qpp::choi2kraus()", Exception::Type::ZERO_SIZE);
906  throw Exception("qpp::choi2kraus()",
908  idx D = internal::get_dim_subsys(A.rows(), 2);
909  // check equal dimensions
910  if (D * D != static_cast<idx>(A.rows()))
911  throw Exception("qpp::choi2kraus()", Exception::Type::DIMS_INVALID);
912  // END EXCEPTION CHECKS
913 
914  dmat ev = hevals(A);
915  cmat evec = hevects(A);
916  std::vector<cmat> result;
917 
918  for (idx i = 0; i < D * D; ++i)
919  {
920  if (std::abs(ev(i)) > eps)
921  result.push_back(
922  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
923  }
924 
925  return result;
926 }
927 
935 inline cmat choi2super(const cmat& A)
936 {
937  // EXCEPTION CHECKS
938 
940  throw Exception("qpp::choi2super()", Exception::Type::ZERO_SIZE);
942  throw Exception("qpp::choi2super()",
944  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
945  // check equal dimensions
946  if (D * D != static_cast<idx>(A.rows()))
947  throw Exception("qpp::choi2super()", Exception::Type::DIMS_INVALID);
948  // END EXCEPTION CHECKS
949 
950  cmat result(D * D, D * D);
951 
952 #ifdef WITH_OPENMP_
953 #pragma omp parallel for collapse(4)
954 #endif // WITH_OPENMP_
955  for (idx a = 0; a < D; ++a)
956  for (idx b = 0; b < D; ++b)
957  for (idx m = 0; m < D; ++m)
958  for (idx n = 0; n < D; ++n)
959  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
960 
961  return result;
962 }
963 
971 inline cmat super2choi(const cmat& A)
972 {
973  // EXCEPTION CHECKS
974 
976  throw Exception("qpp::super2choi()", Exception::Type::ZERO_SIZE);
978  throw Exception("qpp::super2choi()",
980  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()),2);
981  // check equal dimensions
982  if (D * D != static_cast<idx>(A.rows()))
983  throw Exception("qpp::super2choi()", Exception::Type::DIMS_INVALID);
984  // END EXCEPTION CHECKS
985 
986  cmat result(D * D, D * D);
987 
988 #ifdef WITH_OPENMP_
989 #pragma omp parallel for collapse(4)
990 #endif // WITH_OPENMP_
991  for (idx a = 0; a < D; ++a)
992  for (idx b = 0; b < D; ++b)
993  for (idx m = 0; m < D; ++m)
994  for (idx n = 0; n < D; ++n)
995  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
996 
997  return result;
998 }
999 
1014 template<typename Derived>
1015 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1016  const std::vector<idx>& dims)
1017 {
1018  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1019  = A.derived();
1020 
1021  // EXCEPTION CHECKS
1022 
1023  // check zero-size
1025  throw Exception("qpp::ptrace1()", Exception::Type::ZERO_SIZE);
1026 
1027  // check that dims is a valid dimension vector
1028  if (!internal::check_dims(dims))
1029  throw Exception("qpp::ptrace1()", Exception::Type::DIMS_INVALID);
1030 
1031  // check dims has only 2 elements
1032  if (dims.size() != 2)
1033  throw Exception("qpp::ptrace1()", Exception::Type::NOT_BIPARTITE);
1034  // END EXCEPTION CHECKS
1035 
1036  idx DA = dims[0];
1037  idx DB = dims[1];
1038 
1041 
1042  //************ ket ************//
1043  if (internal::check_cvector(rA)) // we have a ket
1044  {
1045  // check that dims match the dimension of A
1046  if (!internal::check_dims_match_cvect(dims, rA))
1047  throw Exception("qpp::ptrace1()",
1049 
1050  auto worker = [=](idx i, idx j) noexcept
1051  -> typename Derived::Scalar
1052  {
1053  typename Derived::Scalar sum = 0;
1054  for (idx m = 0; m < DA; ++m)
1055  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
1056 
1057  return sum;
1058  }; /* end worker */
1059 
1060 #ifdef WITH_OPENMP_
1061 #pragma omp parallel for collapse(2)
1062 #endif // WITH_OPENMP_
1063  for (idx j = 0; j < DB; ++j) // column major order for speed
1064  for (idx i = 0; i < DB; ++i)
1065  result(i, j) = worker(i, j);
1066 
1067  return result;
1068  }
1069  //************ density matrix ************//
1070  else if (internal::check_square_mat(rA)) // we have a density operator
1071  {
1072  // check that dims match the dimension of A
1073  if (!internal::check_dims_match_mat(dims, rA))
1074  throw Exception("qpp::ptrace1()",
1076 
1077  auto worker = [=](idx i, idx j) noexcept
1078  -> typename Derived::Scalar
1079  {
1080  typename Derived::Scalar sum = 0;
1081  for (idx m = 0; m < DA; ++m)
1082  sum += rA(m * DB + i, m * DB + j);
1083 
1084  return sum;
1085  }; /* end worker */
1086 
1087 #ifdef WITH_OPENMP_
1088 #pragma omp parallel for collapse(2)
1089 #endif // WITH_OPENMP_
1090  for (idx j = 0; j < DB; ++j) // column major order for speed
1091  for (idx i = 0; i < DB; ++i)
1092  result(i, j) = worker(i, j);
1093 
1094  return result;
1095  }
1096  //************ Exception: not ket nor density matrix ************//
1097  else
1098  throw Exception("qpp::ptrace1()",
1100 }
1101 
1116 template<typename Derived>
1117 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1118  const std::vector<idx>& dims)
1119 {
1120  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1121  = A.derived();
1122 
1123  // EXCEPTION CHECKS
1124 
1125  // check zero-size
1127  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1128 
1129  // check that dims is a valid dimension vector
1130  if (!internal::check_dims(dims))
1131  throw Exception("qpp::ptrace2()", Exception::Type::DIMS_INVALID);
1132 
1133  // check dims has only 2 elements
1134  if (dims.size() != 2)
1135  throw Exception("qpp::ptrace2()", Exception::Type::NOT_BIPARTITE);
1136  // END EXCEPTION CHECKS
1137 
1138  idx DA = dims[0];
1139  idx DB = dims[1];
1140 
1143 
1144  //************ ket ************//
1145  if (internal::check_cvector(rA)) // we have a ket
1146  {
1147  // check that dims match the dimension of A
1148  if (!internal::check_dims_match_cvect(dims, rA))
1149  throw Exception("qpp::ptrace2()",
1151 
1152  auto worker = [=](idx i, idx j) noexcept
1153  -> typename Derived::Scalar
1154  {
1155  typename Derived::Scalar sum = 0;
1156  for (idx m = 0; m < DB; ++m)
1157  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1158 
1159  return sum;
1160  }; /* end worker */
1161 
1162 #ifdef WITH_OPENMP_
1163 #pragma omp parallel for collapse(2)
1164 #endif // WITH_OPENMP_
1165  for (idx j = 0; j < DA; ++j) // column major order for speed
1166  for (idx i = 0; i < DA; ++i)
1167  result(i, j) = worker(i, j);
1168 
1169  return result;
1170  }
1171  //************ density matrix ************//
1172  else if (internal::check_square_mat(rA)) // we have a density operator
1173  {
1174  // check that dims match the dimension of A
1175  if (!internal::check_dims_match_mat(dims, rA))
1176  throw Exception("qpp::ptrace2()",
1178 
1179 #ifdef WITH_OPENMP_
1180 #pragma omp parallel for collapse(2)
1181 #endif // WITH_OPENMP_
1182  for (idx j = 0; j < DA; ++j) // column major order for speed
1183  for (idx i = 0; i < DA; ++i)
1184  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1185 
1186  return result;
1187  }
1188  //************ Exception: not ket nor density matrix ************//
1189  else
1190  throw Exception("qpp::ptrace1()",
1192 }
1193 
1208 template<typename Derived>
1209 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1210  const std::vector<idx>& subsys,
1211  const std::vector<idx>& dims)
1212 {
1213  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1214  = A.derived();
1215 
1216  // EXCEPTION CHECKS
1217 
1218  // check zero-size
1220  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1221 
1222  // check that dims is a valid dimension vector
1223  if (!internal::check_dims(dims))
1224  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1225 
1226  // check that subsys are valid
1227  if (!internal::check_subsys_match_dims(subsys, dims))
1228  throw Exception("qpp::ptrace()",
1230  // END EXCEPTION CHECKS
1231 
1232  idx D = static_cast<idx>(rA.rows());
1233  idx N = dims.size();
1234  idx Nsubsys = subsys.size();
1235  idx Nsubsys_bar = N - Nsubsys;
1236  idx Dsubsys = 1;
1237  for (idx i = 0; i < Nsubsys; ++i)
1238  Dsubsys *= dims[subsys[i]];
1239  idx Dsubsys_bar = D / Dsubsys;
1240 
1241  idx Cdims[maxn];
1242  idx Csubsys[maxn];
1243  idx Cdimssubsys[maxn];
1244  idx Csubsys_bar[maxn];
1245  idx Cdimssubsys_bar[maxn];
1246 
1247  idx Cmidxcolsubsys_bar[maxn];
1248 
1249  std::vector<idx> subsys_bar = complement(subsys, N);
1250  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1251  std::begin(Csubsys_bar));
1252 
1253  for (idx i = 0; i < N; ++i)
1254  {
1255  Cdims[i] = dims[i];
1256  }
1257  for (idx i = 0; i < Nsubsys; ++i)
1258  {
1259  Csubsys[i] = subsys[i];
1260  Cdimssubsys[i] = dims[subsys[i]];
1261  }
1262  for (idx i = 0; i < Nsubsys_bar; ++i)
1263  {
1264  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1265  }
1266 
1268  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1269 
1270  //************ ket ************//
1271  if (internal::check_cvector(rA)) // we have a ket
1272  {
1273  // check that dims match the dimension of A
1274  if (!internal::check_dims_match_cvect(dims, rA))
1275  throw Exception("qpp::ptrace()",
1277 
1278  if (subsys.size() == dims.size())
1279  {
1280  result(0, 0) = (adjoint(rA) * rA).value();
1281  return result;
1282  }
1283 
1284  if (subsys.size() == 0)
1285  return rA * adjoint(rA);
1286 
1287  auto worker = [=, &Cmidxcolsubsys_bar](idx i) noexcept
1288  -> typename Derived::Scalar
1289  {
1290  // use static allocation for speed!
1291 
1292  idx Cmidxrow[maxn];
1293  idx Cmidxcol[maxn];
1294  idx Cmidxrowsubsys_bar[maxn];
1295  idx Cmidxsubsys[maxn];
1296 
1297  /* get the row multi-indexes of the complement */
1298  internal::n2multiidx(i, Nsubsys_bar,
1299  Cdimssubsys_bar, Cmidxrowsubsys_bar);
1300  /* write them in the global row/col multi-indexes */
1301  for (idx k = 0; k < Nsubsys_bar; ++k)
1302  {
1303  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1304  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1305  }
1306  typename Derived::Scalar sm = 0;
1307  for (idx a = 0; a < Dsubsys; ++a)
1308  {
1309  // get the multi-index over which we do the summation
1310  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1311  // write it into the global row/col multi-indexes
1312  for (idx k = 0; k < Nsubsys; ++k)
1313  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1314  = Cmidxsubsys[k];
1315 
1316  // now do the sum
1317  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims)) *
1318  std::conj(rA(internal::multiidx2n(Cmidxcol, N,
1319  Cdims)));
1320  }
1321 
1322  return sm;
1323  }; /* end worker */
1324 
1325  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1326  {
1327  // compute the column multi-indexes of the complement
1328  internal::n2multiidx(j, Nsubsys_bar,
1329  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1330 #ifdef WITH_OPENMP_
1331 #pragma omp parallel for
1332 #endif // WITH_OPENMP_
1333  for (idx i = 0; i < Dsubsys_bar; ++i)
1334  {
1335  result(i, j) = worker(i);
1336  }
1337  }
1338 
1339  return result;
1340  }
1341  //************ density matrix ************//
1342  else if (internal::check_square_mat(rA)) // we have a density operator
1343  {
1344  // check that dims match the dimension of A
1345  if (!internal::check_dims_match_mat(dims, rA))
1346  throw Exception("qpp::ptrace()",
1348 
1349  if (subsys.size() == dims.size())
1350  {
1351  result(0, 0) = rA.trace();
1352  return result;
1353  }
1354 
1355  if (subsys.size() == 0)
1356  return rA;
1357 
1358  auto worker = [=, &Cmidxcolsubsys_bar](idx i) noexcept
1359  -> typename Derived::Scalar
1360  {
1361  // use static allocation for speed!
1362 
1363  idx Cmidxrow[maxn];
1364  idx Cmidxcol[maxn];
1365  idx Cmidxrowsubsys_bar[maxn];
1366  idx Cmidxsubsys[maxn];
1367 
1368  /* get the row/col multi-indexes of the complement */
1369  internal::n2multiidx(i, Nsubsys_bar,
1370  Cdimssubsys_bar, Cmidxrowsubsys_bar);
1371  /* write them in the global row/col multi-indexes */
1372  for (idx k = 0; k < Nsubsys_bar; ++k)
1373  {
1374  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1375  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1376  }
1377  typename Derived::Scalar sm = 0;
1378  for (idx a = 0; a < Dsubsys; ++a)
1379  {
1380  // get the multi-index over which we do the summation
1381  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1382  // write it into the global row/col multi-indexes
1383  for (idx k = 0; k < Nsubsys; ++k)
1384  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1385  = Cmidxsubsys[k];
1386 
1387  // now do the sum
1388  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims),
1389  internal::multiidx2n(Cmidxcol, N, Cdims));
1390  }
1391 
1392  return sm;
1393  }; /* end worker */
1394 
1395  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1396  {
1397  // compute the column multi-indexes of the complement
1398  internal::n2multiidx(j, Nsubsys_bar,
1399  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1400 #ifdef WITH_OPENMP_
1401 #pragma omp parallel for
1402 #endif // WITH_OPENMP_
1403  for (idx i = 0; i < Dsubsys_bar; ++i)
1404  {
1405  result(i, j) = worker(i);
1406  }
1407  }
1408 
1409  return result;
1410  }
1411  //************ Exception: not ket nor density matrix ************//
1412  else
1413  throw Exception("qpp::ptrace()",
1415 }
1416 
1431 template<typename Derived>
1432 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1433  const std::vector<idx>& subsys,
1434  idx d = 2)
1435 {
1436  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1437 
1438  // EXCEPTION CHECKS
1439 
1440  // check zero size
1442  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1443 
1444  // check valid dims
1445  if (d == 0)
1446  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1447  // END EXCEPTION CHECKS
1448 
1449  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1450  std::vector<idx> dims(N, d); // local dimensions vector
1451 
1452  return ptrace(rA, subsys, dims);
1453 }
1454 
1468 template<typename Derived>
1470  const Eigen::MatrixBase<Derived>& A,
1471  const std::vector<idx>& subsys,
1472  const std::vector<idx>& dims)
1473 {
1474  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1475 
1476  // EXCEPTION CHECKS
1477 
1478  // check zero-size
1480  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1481 
1482  // check that dims is a valid dimension vector
1483  if (!internal::check_dims(dims))
1484  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1485 
1486  // check that subsys are valid
1487  if (!internal::check_subsys_match_dims(subsys, dims))
1488  throw Exception("qpp::ptranspose()",
1490  // END EXCEPTION CHECKS
1491 
1492  idx D = static_cast<idx>(rA.rows());
1493  idx N = dims.size();
1494  idx Nsubsys = subsys.size();
1495  idx Cdims[maxn];
1496  idx Cmidxcol[maxn];
1497  idx Csubsys[maxn];
1498 
1499  // copy dims in Cdims and subsys in Csubsys
1500  for (idx i = 0; i < N; ++i)
1501  Cdims[i] = dims[i];
1502  for (idx i = 0; i < Nsubsys; ++i)
1503  Csubsys[i] = subsys[i];
1504 
1505  dyn_mat<typename Derived::Scalar> result(D, D);
1506 
1507  //************ ket ************//
1508  if (internal::check_cvector(rA)) // we have a ket
1509  {
1510  // check that dims match the dimension of A
1511  if (!internal::check_dims_match_cvect(dims, rA))
1512  throw Exception("qpp::ptranspose()",
1514 
1515  if (subsys.size() == dims.size())
1516  return (rA * adjoint(rA)).transpose();
1517 
1518  if (subsys.size() == 0)
1519  return rA * adjoint(rA);
1520 
1521  auto worker = [=, &Cmidxcol](idx i) noexcept
1522  -> typename Derived::Scalar
1523  {
1524  // use static allocation for speed!
1525  idx midxcoltmp[maxn];
1526  idx midxrow[maxn];
1527 
1528  for (idx k = 0; k < N; ++k)
1529  midxcoltmp[k] = Cmidxcol[k];
1530 
1531  /* compute the row multi-index */
1532  internal::n2multiidx(i, N, Cdims, midxrow);
1533 
1534  for (idx k = 0; k < Nsubsys; ++k)
1535  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1536 
1537  /* writes the result */
1538  return rA(internal::multiidx2n(midxrow, N, Cdims)) *
1539  std::conj(rA(internal::multiidx2n(midxcoltmp, N,
1540  Cdims)));
1541  }; /* end worker */
1542 
1543  for (idx j = 0; j < D; ++j)
1544  {
1545  // compute the column multi-index
1546  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1547 
1548 #ifdef WITH_OPENMP_
1549 #pragma omp parallel for
1550 #endif // WITH_OPENMP_
1551  for (idx i = 0; i < D; ++i)
1552  result(i, j) = worker(i);
1553  }
1554 
1555  return result;
1556  }
1557  //************ density matrix ************//
1558  else if (internal::check_square_mat(rA)) // we have a density operator
1559  {
1560  // check that dims match the dimension of A
1561  if (!internal::check_dims_match_mat(dims, rA))
1562  throw Exception("qpp::ptranspose()",
1564 
1565  if (subsys.size() == dims.size())
1566  return rA.transpose();
1567 
1568  if (subsys.size() == 0)
1569  return rA;
1570 
1571  auto worker = [=, &Cmidxcol](idx i) noexcept
1572  -> typename Derived::Scalar
1573  {
1574  // use static allocation for speed!
1575  idx midxcoltmp[maxn];
1576  idx midxrow[maxn];
1577 
1578  for (idx k = 0; k < N; ++k)
1579  midxcoltmp[k] = Cmidxcol[k];
1580 
1581  /* compute the row multi-index */
1582  internal::n2multiidx(i, N, Cdims, midxrow);
1583 
1584  for (idx k = 0; k < Nsubsys; ++k)
1585  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1586 
1587  /* writes the result */
1588  return rA(internal::multiidx2n(midxrow, N, Cdims),
1589  internal::multiidx2n(midxcoltmp, N, Cdims));
1590  }; /* end worker */
1591 
1592  for (idx j = 0; j < D; ++j)
1593  {
1594  // compute the column multi-index
1595  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1596 
1597 #ifdef WITH_OPENMP_
1598 #pragma omp parallel for
1599 #endif // WITH_OPENMP_
1600  for (idx i = 0; i < D; ++i)
1601  result(i, j) = worker(i);
1602  }
1603 
1604  return result;
1605  }
1606  //************ Exception: not ket nor density matrix ************//
1607  else
1608  throw Exception("qpp::ptranspose()",
1610 }
1611 
1625 template<typename Derived>
1627  const Eigen::MatrixBase<Derived>& A,
1628  const std::vector<idx>& subsys,
1629  idx d = 2)
1630 {
1631  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1632 
1633  // EXCEPTION CHECKS
1634 
1635  // check zero size
1637  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1638 
1639  // check valid dims
1640  if (d == 0)
1641  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1642  // END EXCEPTION CHECKS
1643 
1644  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1645  std::vector<idx> dims(N, d); // local dimensions vector
1646 
1647  return ptranspose(rA, subsys, dims);
1648 }
1649 
1662 template<typename Derived>
1664  const Eigen::MatrixBase<Derived>& A,
1665  const std::vector<idx>& perm,
1666  const std::vector<idx>& dims)
1667 {
1668  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1669 
1670  // EXCEPTION CHECKS
1671 
1672  // check zero-size
1674  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1675 
1676  // check that dims is a valid dimension vector
1677  if (!internal::check_dims(dims))
1678  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1679 
1680  // check that we have a valid permutation
1681  if (!internal::check_perm(perm))
1682  throw Exception("qpp::syspermute()", Exception::Type::PERM_INVALID);
1683 
1684  // check that permutation match dimensions
1685  if (perm.size() != dims.size())
1686  throw Exception("qpp::syspermute()",
1688  // END EXCEPTION CHECKS
1689 
1690  idx D = static_cast<idx>(rA.rows());
1691  idx N = dims.size();
1692 
1694 
1695  //************ ket ************//
1696  if (internal::check_cvector(rA)) // we have a column vector
1697  {
1698  idx Cdims[maxn];
1699  idx Cperm[maxn];
1700 
1701  // check that dims match the dimension of rA
1702  if (!internal::check_dims_match_cvect(dims, rA))
1703  throw Exception("qpp::syspermute()",
1705 
1706  // copy dims in Cdims and perm in Cperm
1707  for (idx i = 0; i < N; ++i)
1708  {
1709  Cdims[i] = dims[i];
1710  Cperm[i] = perm[i];
1711  }
1712  result.resize(D, 1);
1713 
1714  auto worker = [&Cdims, &Cperm, N](idx i) noexcept
1715  -> idx
1716  {
1717  // use static allocation for speed,
1718  // double the size for matrices reshaped as vectors
1719  idx midx[maxn];
1720  idx midxtmp[maxn];
1721  idx permdims[maxn];
1722 
1723  /* compute the multi-index */
1724  internal::n2multiidx(i, N, Cdims, midx);
1725 
1726  for (idx k = 0; k < N; ++k)
1727  {
1728  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1729  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1730  }
1731 
1732  return internal::multiidx2n(midxtmp, N, permdims);
1733  }; /* end worker */
1734 
1735 #ifdef WITH_OPENMP_
1736 #pragma omp parallel for
1737 #endif // WITH_OPENMP_
1738  for (idx i = 0; i < D; ++i)
1739  result(worker(i)) = rA(i);
1740 
1741  return result;
1742  }
1743  //************ density matrix ************//
1744  else if (internal::check_square_mat(rA)) // we have a density operator
1745  {
1746  idx Cdims[2 * maxn];
1747  idx Cperm[2 * maxn];
1748 
1749  // check that dims match the dimension of rA
1750  if (!internal::check_dims_match_mat(dims, rA))
1751  throw Exception("qpp::syspermute()",
1753 
1754  // copy dims in Cdims and perm in Cperm
1755  for (idx i = 0; i < N; ++i)
1756  {
1757  Cdims[i] = dims[i];
1758  Cdims[i + N] = dims[i];
1759  Cperm[i] = perm[i];
1760  Cperm[i + N] = perm[i] + N;
1761  }
1762  result.resize(D * D, 1);
1763  // map A to a column vector
1765  Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1766  const_cast<typename Derived::Scalar*>(rA.data()), D *
1767  D,
1768  1);
1769 
1770  auto worker = [&Cdims, &Cperm, N](idx i) noexcept
1771  -> idx
1772  {
1773  // use static allocation for speed,
1774  // double the size for matrices reshaped as vectors
1775  idx midx[2 * maxn];
1776  idx midxtmp[2 * maxn];
1777  idx permdims[2 * maxn];
1778 
1779  /* compute the multi-index */
1780  internal::n2multiidx(i, 2 * N, Cdims, midx);
1781 
1782  for (idx k = 0; k < 2 * N; ++k)
1783  {
1784  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1785  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1786  }
1787 
1788  return internal::multiidx2n(midxtmp, 2 * N, permdims);
1789  }; /* end worker */
1790 
1791 #ifdef WITH_OPENMP_
1792 #pragma omp parallel for
1793 #endif // WITH_OPENMP_
1794  for (idx i = 0; i < D * D; ++i)
1795  result(worker(i)) = rA(i);
1796 
1797  return reshape(result, D, D);
1798  }
1799  //************ Exception: not ket nor density matrix ************//
1800  else
1801  throw Exception("qpp::syspermute()",
1803 }
1804 
1817 template<typename Derived>
1819  const Eigen::MatrixBase<Derived>& A,
1820  const std::vector<idx>& perm,
1821  idx d = 2)
1822 {
1823  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1824 
1825  // EXCEPTION CHECKS
1826 
1827  // check zero size
1829  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1830 
1831  // check valid dims
1832  if (d == 0)
1833  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1834  // END EXCEPTION CHECKS
1835 
1836  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1837  std::vector<idx> dims(N, d); // local dimensions vector
1838 
1839  return syspermute(rA, perm, dims);
1840 }
1841 
1842 } /* namespace qpp */
1843 
1844 #endif /* OPERATIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:112
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:74
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:899
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:66
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:67
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:56
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:1663
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:78
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:51
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:485
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:242
idx get_num_subsys(idx sz, idx d)
Definition: util.h:355
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:387
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:152
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:935
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:971
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1814
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1015
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:895
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1117
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:44
void n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:48
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:778
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:363
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:840
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:1209
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:127
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:762
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:54
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:1469
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:125
std::size_t idx
Non-negative integer index.
Definition: types.h:36
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:209
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1135
idx multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:61