Quantum++  v1.0.0-beta2
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 // 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 {
37 
54 template<typename Derived1, typename Derived2>
56  const Eigen::MatrixBase<Derived1>& state,
57  const Eigen::MatrixBase<Derived2>& A,
58  const std::vector<idx>& ctrl,
59  const std::vector<idx>& subsys,
60  const std::vector<idx>& dims)
61 {
62  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
63  = state.derived();
64  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
65 
66  // EXCEPTION CHECKS
67 
68  // check types
69  if (!std::is_same<typename Derived1::Scalar,
70  typename Derived2::Scalar>::value)
71  throw Exception("qpp::applyCTRL()", Exception::Type::TYPE_MISMATCH);
72 
73  // check zero sizes
75  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
76 
77  // check zero sizes
78  if (!internal::check_nonzero_size(rstate))
79  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
80 
81  // check square matrix for the gate
83  throw Exception("qpp::applyCTRL()",
85 
86  // check that all control subsystems have the same dimension
87  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
88  for (idx i = 1; i < ctrl.size(); ++i)
89  if (dims[ctrl[i]] != d)
90  throw Exception("qpp::applyCTRL()",
92 
93  // check that dimension is valid
94  if (!internal::check_dims(dims))
95  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
96 
97  // check subsys is valid w.r.t. dims
98  if (!internal::check_subsys_match_dims(subsys, dims))
99  throw Exception("qpp::applyCTRL()",
101 
102  // check that gate matches the dimensions of the subsys
103  std::vector<idx> subsys_dims(subsys.size());
104  for (idx i = 0; i < subsys.size(); ++i)
105  subsys_dims[i] = dims[subsys[i]];
106  if (!internal::check_dims_match_mat(subsys_dims, rA))
107  throw Exception("qpp::applyCTRL()",
109 
110  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
111  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
112  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
113 
114  // check that ctrl + gate subsystem is valid
115  // with respect to local dimensions
116  if (!internal::check_subsys_match_dims(ctrlgate, dims))
117  throw Exception("qpp::applyCTRL()",
119  // END EXCEPTION CHECKS
120 
121  // construct the table of A^i and (A^dagger)^i
122  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
123  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
124  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i)
125  {
126  Ai.push_back(powm(rA, i));
127  Aidagger.push_back(powm(adjoint(rA), i));
128  }
129 
130  idx D = static_cast<idx>(rstate.rows()); // total dimension
131  idx N = dims.size(); // total number of subsystems
132  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
133  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
134  idx subsyssize = subsys.size(); // number of subsystems of the target
135  // dimension of ctrl subsystem
136  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
137  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
138 
139  idx Cdims[maxn]; // local dimensions
140  idx CdimsA[maxn]; // local dimensions
141  idx CdimsCTRL[maxn]; // local dimensions
142  idx CdimsCTRLA_bar[maxn]; // local dimensions
143 
144  // compute the complementary subsystem of ctrlgate w.r.t. dims
145  std::vector<idx> ctrlgate_bar = complement(ctrlgate, N);
146  // number of subsystems that are complementary to the ctrl+gate
147  idx ctrlgate_barsize = ctrlgate_bar.size();
148 
149  idx DCTRLA_bar = 1; // dimension of the rest
150  for (idx i = 0; i < ctrlgate_barsize; ++i)
151  DCTRLA_bar *= dims[ctrlgate_bar[i]];
152 
153  for (idx k = 0; k < N; ++k)
154  Cdims[k] = dims[k];
155  for (idx k = 0; k < subsyssize; ++k)
156  CdimsA[k] = dims[subsys[k]];
157  for (idx k = 0; k < ctrlsize; ++k)
158  CdimsCTRL[k] = d;
159  for (idx k = 0; k < ctrlgate_barsize; ++k)
160  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
161 
162  // worker, computes the coefficient and the index for the ket case
163  // used in #pragma omp parallel for collapse
164  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
165  -> std::pair<typename Derived1::Scalar, idx>
166  {
167  idx indx = 0;
168  typename Derived1::Scalar coeff = 0;
169 
170  idx Cmidx[maxn]; // the total multi-index
171  idx CmidxA[maxn]; // the gate part multi-index
172  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
173 
174  // compute the index
175 
176  // set the CTRL part
177  for (idx k = 0; k < ctrlsize; ++k)
178  {
179  Cmidx[ctrl[k]] = i_;
180  }
181 
182  // set the rest
183  internal::n2multiidx(r_, N - ctrlgatesize,
184  CdimsCTRLA_bar, CmidxCTRLA_bar);
185  for (idx k = 0; k < N - ctrlgatesize; ++k)
186  {
187  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
188  }
189 
190  // set the A part
191  internal::n2multiidx(m_, subsyssize, CdimsA, CmidxA);
192  for (idx k = 0; k < subsyssize; ++k)
193  {
194  Cmidx[subsys[k]] = CmidxA[k];
195  }
196 
197  // we now got the total index
198  indx = internal::multiidx2n(Cmidx, N, Cdims);
199 
200  // compute the coefficient
201  for (idx n_ = 0; n_ < DA; ++n_)
202  {
203  internal::n2multiidx(n_, subsyssize, CdimsA, CmidxA);
204  for (idx k = 0; k < subsyssize; ++k)
205  {
206  Cmidx[subsys[k]] = CmidxA[k];
207  }
208  coeff += Ai[i_](m_, n_) *
209  rstate(internal::multiidx2n(Cmidx, N, Cdims));
210  }
211 
212  return std::make_pair(coeff, indx);
213  }; /* end coeff_idx_ket */
214 
215  // worker, computes the coefficient and the index
216  // for the density matrix case
217  // used in #pragma omp parallel for collapse
218  auto coeff_idx_rho = [&](idx i1_, idx m1_,
219  idx r1_, idx i2_, idx m2_, idx r2_) noexcept
220  -> std::tuple<typename Derived1::Scalar, idx, idx>
221  {
222  idx idxrow = 0;
223  idx idxcol = 0;
224  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
225 
226  idx Cmidxrow[maxn]; // the total row multi-index
227  idx Cmidxcol[maxn]; // the total col multi-index
228  idx CmidxArow[maxn]; // the gate part row multi-index
229  idx CmidxAcol[maxn]; // the gate part col multi-index
230  idx CmidxCTRLrow[maxn]; // the control row multi-index
231  idx CmidxCTRLcol[maxn]; // the control col multi-index
232  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
233  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
234 
235  // compute the ket/bra indexes
236 
237  // set the CTRL part
238  internal::n2multiidx(i1_, ctrlsize,
239  CdimsCTRL, CmidxCTRLrow);
240  internal::n2multiidx(i2_, ctrlsize,
241  CdimsCTRL, CmidxCTRLcol);
242 
243  for (idx k = 0; k < ctrlsize; ++k)
244  {
245  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
246  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
247  }
248 
249  // set the rest
250  internal::n2multiidx(r1_, N - ctrlgatesize,
251  CdimsCTRLA_bar, CmidxCTRLA_barrow);
252  internal::n2multiidx(r2_, N - ctrlgatesize,
253  CdimsCTRLA_bar, CmidxCTRLA_barcol);
254  for (idx k = 0; k < N - ctrlgatesize; ++k)
255  {
256  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
257  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
258  }
259 
260  // set the A part
261  internal::n2multiidx(m1_, subsyssize, CdimsA, CmidxArow);
262  internal::n2multiidx(m2_, subsyssize, CdimsA, CmidxAcol);
263  for (idx k = 0; k < subsys.size(); ++k)
264  {
265  Cmidxrow[subsys[k]] = CmidxArow[k];
266  Cmidxcol[subsys[k]] = CmidxAcol[k];
267  }
268 
269  // we now got the total row/col indexes
270  idxrow = internal::multiidx2n(Cmidxrow, N, Cdims);
271  idxcol = internal::multiidx2n(Cmidxcol, N, Cdims);
272 
273  // check whether all CTRL row and col multi indexes are equal
274  bool all_ctrl_rows_equal = true;
275  bool all_ctrl_cols_equal = true;
276 
277  idx first_ctrl_row, first_ctrl_col;
278  if (ctrlsize > 0)
279  {
280  first_ctrl_row = CmidxCTRLrow[0];
281  first_ctrl_col = CmidxCTRLcol[0];
282  } else
283  {
284  first_ctrl_row = first_ctrl_col = 1;
285  }
286 
287  for (idx k = 1; k < ctrlsize; ++k)
288  {
289  if (CmidxCTRLrow[k] != first_ctrl_row)
290  {
291  all_ctrl_rows_equal = false;
292  break;
293  }
294  }
295  for (idx k = 1; k < ctrlsize; ++k)
296  {
297  if (CmidxCTRLcol[k] != first_ctrl_col)
298  {
299  all_ctrl_cols_equal = false;
300  break;
301  }
302  }
303 
304  // at least one control activated, compute the coefficient
305  for (idx n1_ = 0; n1_ < DA; ++n1_)
306  {
307  internal::n2multiidx(n1_, subsyssize, CdimsA, CmidxArow);
308  for (idx k = 0; k < subsyssize; ++k)
309  {
310  Cmidxrow[subsys[k]] = CmidxArow[k];
311  }
312  idx idxrowtmp = internal::multiidx2n(Cmidxrow, N, Cdims);
313 
314  if (all_ctrl_rows_equal)
315  {
316  lhs = Ai[first_ctrl_row](m1_, n1_);
317  } else
318  {
319  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
320  }
321 
322  for (idx n2_ = 0; n2_ < DA; ++n2_)
323  {
324  internal::n2multiidx(n2_, subsyssize, CdimsA, CmidxAcol);
325  for (idx k = 0; k < subsyssize; ++k)
326  {
327  Cmidxcol[subsys[k]] = CmidxAcol[k];
328  }
329 
330  if (all_ctrl_cols_equal)
331  {
332  rhs = Aidagger[first_ctrl_col](n2_, m2_);
333  } else
334  {
335  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
336  }
337 
338  idx idxcoltmp = internal::multiidx2n(Cmidxcol, N, Cdims);
339 
340  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
341  }
342  }
343 
344  return std::make_tuple(coeff, idxrow, idxcol);
345  }; /* end coeff_idx_rho */
346 
347  //************ ket ************//
348  if (internal::check_cvector(rstate)) // we have a ket
349  {
350  // check that dims match state vector
351  if (!internal::check_dims_match_cvect(dims, rstate))
352  throw Exception("qpp::applyCTRL()",
354  if (D == 1)
355  return rstate;
356 
357  dyn_mat<typename Derived1::Scalar> result = rstate;
358 
359 #ifdef WITH_OPENMP_
360 #pragma omp parallel for collapse(2)
361 #endif // WITH_OPENMP_
362  for (idx m = 0; m < DA; ++m)
363  for (idx r = 0; r < DCTRLA_bar; ++r)
364  {
365  if (ctrlsize == 0) // no control
366  {
367  result(coeff_idx_ket(1, m, r).second) =
368  coeff_idx_ket(1, m, r).first;
369  } else
370  for (idx i = 0; i < d; ++i)
371  {
372  result(coeff_idx_ket(i, m, r).second) =
373  coeff_idx_ket(i, m, r).first;
374  }
375  }
376 
377  return result;
378  }
379  //************ density matrix ************//
380  else if (internal::check_square_mat(rstate)) // we have a density operator
381  {
382  // check that dims match state matrix
383  if (!internal::check_dims_match_mat(dims, rstate))
384  throw Exception("qpp::applyCTRL()",
386 
387  if (D == 1)
388  return rstate;
389 
390  dyn_mat<typename Derived1::Scalar> result = rstate;
391 
392 #ifdef WITH_OPENMP_
393 #pragma omp parallel for collapse(4)
394 #endif // WITH_OPENMP_
395  for (idx m1 = 0; m1 < DA; ++m1)
396  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
397  for (idx m2 = 0; m2 < DA; ++m2)
398  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
399  if (ctrlsize == 0) // no control
400  {
401  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
402  1, m2, r2);
403  result(std::get<1>(coeff_idxes),
404  std::get<2>(coeff_idxes)) =
405  std::get<0>(coeff_idxes);
406  } else
407  {
408  for (idx i1 = 0; i1 < Dctrl; ++i1)
409  for (idx i2 = 0; i2 < Dctrl; ++i2)
410  {
411  auto coeff_idxes = coeff_idx_rho(
412  i1, m1, r1,
413  i2, m2, r2);
414  result(std::get<1>(coeff_idxes),
415  std::get<2>(coeff_idxes)) =
416  std::get<0>(coeff_idxes);
417  }
418  }
419 
420  return result;
421  }
422  //************ Exception: not ket nor density matrix ************//
423  else
424  throw Exception("qpp::applyCTRL()",
426 }
427 
443 template<typename Derived1, typename Derived2>
445  const Eigen::MatrixBase<Derived1>& state,
446  const Eigen::MatrixBase<Derived2>& A,
447  const std::vector<idx>& ctrl,
448  const std::vector<idx>& subsys,
449  idx d = 2)
450 {
451  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
452  = state.derived();
453  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
454 
455  // EXCEPTION CHECKS
456 
457  // check zero size
458  if (!internal::check_nonzero_size(rstate))
459  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
460 
461  // check valid dims
462  if (d == 0)
463  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
464  // END EXCEPTION CHECKS
465 
466  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
467  std::vector<idx> dims(N, d); // local dimensions vector
468 
469  return applyCTRL(rstate, rA, ctrl, subsys, dims);
470 }
471 
485 template<typename Derived1, typename Derived2>
487  const Eigen::MatrixBase<Derived1>& state,
488  const Eigen::MatrixBase<Derived2>& A,
489  const std::vector<idx>& subsys,
490  const std::vector<idx>& dims)
491 {
492  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
493  = state.derived();
494  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
495 
496  // EXCEPTION CHECKS
497 
498  // check types
499  if (!std::is_same<typename Derived1::Scalar,
500  typename Derived2::Scalar>::value)
501  throw Exception("qpp::apply()", Exception::Type::TYPE_MISMATCH);
502 
503  // check zero sizes
505  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
506 
507  // check zero sizes
508  if (!internal::check_nonzero_size(rstate))
509  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
510 
511  // check square matrix for the gate
513  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
514 
515  // check that dimension is valid
516  if (!internal::check_dims(dims))
517  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
518 
519  // check subsys is valid w.r.t. dims
520  if (!internal::check_subsys_match_dims(subsys, dims))
521  throw Exception("qpp::apply()", Exception::Type::SUBSYS_MISMATCH_DIMS);
522 
523  // check that gate matches the dimensions of the subsys
524  std::vector<idx> subsys_dims(subsys.size());
525  for (idx i = 0; i < subsys.size(); ++i)
526  subsys_dims[i] = dims[subsys[i]];
527  if (!internal::check_dims_match_mat(subsys_dims, rA))
528  throw Exception("qpp::apply()",
530  // END EXCEPTION CHECKS
531 
532  //************ ket ************//
533  if (internal::check_cvector(rstate)) // we have a ket
534  {
535  // check that dims match state vector
536  if (!internal::check_dims_match_cvect(dims, rstate))
537  throw Exception("qpp::apply()",
539 
540  return applyCTRL(rstate, rA, {}, subsys, dims);
541  }
542  //************ density matrix ************//
543  else if (internal::check_square_mat(rstate)) // we have a density operator
544  {
545 
546  // check that dims match state matrix
547  if (!internal::check_dims_match_mat(dims, rstate))
548  throw Exception("qpp::apply()",
550 
551  return applyCTRL(rstate, rA, {}, subsys, dims);
552  }
553  //************ Exception: not ket nor density matrix ************//
554  else
555  throw Exception("qpp::apply()",
557 }
558 
572 template<typename Derived1, typename Derived2>
574  const Eigen::MatrixBase<Derived1>& state,
575  const Eigen::MatrixBase<Derived2>& A,
576  const std::vector<idx>& subsys,
577  idx d = 2)
578 {
579  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
580  = state.derived();
581  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
582 
583  // EXCEPTION CHECKS
584 
585  // check zero size
586  if (!internal::check_nonzero_size(rstate))
587  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
588 
589  // check valid dims
590  if (d == 0)
591  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
592  // END EXCEPTION CHECKS
593 
594  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
595  std::vector<idx> dims(N, d); // local dimensions vector
596 
597  return apply(rstate, rA, subsys, dims);
598 }
599 
608 template<typename Derived>
609 cmat apply(const Eigen::MatrixBase<Derived>& A,
610  const std::vector<cmat>& Ks)
611 {
612  const cmat& rA = A.derived();
613 
614  // EXCEPTION CHECKS
615 
617  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
619  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
620  if (Ks.size() == 0)
621  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
622  if (!internal::check_square_mat(Ks[0]))
623  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
624  if (Ks[0].rows() != rA.rows())
625  throw Exception("qpp::apply()",
627  for (auto&& it : Ks)
628  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
629  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
630  // END EXCEPTION CHECKS
631 
632  cmat result = cmat::Zero(rA.rows(), rA.rows());
633 
634 #ifdef WITH_OPENMP_
635 #pragma omp parallel for
636 #endif // WITH_OPENMP_
637  for (idx i = 0; i < Ks.size(); ++i)
638  {
639 #ifdef WITH_OPENMP_
640 #pragma omp critical
641 #endif // WITH_OPENMP_
642  {
643  result += Ks[i] * rA * adjoint(Ks[i]);
644  }
645  }
646 
647  return result;
648 }
649 
660 template<typename Derived>
661 cmat apply(const Eigen::MatrixBase<Derived>& A,
662  const std::vector<cmat>& Ks,
663  const std::vector<idx>& subsys,
664  const std::vector<idx>& dims)
665 {
666  const cmat& rA = A.derived();
667 
668  // EXCEPTION CHECKS
669 
670  // check zero sizes
672  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
673 
674  // check square matrix for the A
676  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
677 
678  // check that dimension is valid
679  if (!internal::check_dims(dims))
680  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
681 
682  // check that dims match A matrix
683  if (!internal::check_dims_match_mat(dims, rA))
684  throw Exception("qpp::apply()",
686 
687  // check subsys is valid w.r.t. dims
688  if (!internal::check_subsys_match_dims(subsys, dims))
689  throw Exception("qpp::apply()",
691 
692  std::vector<idx> subsys_dims(subsys.size());
693  for (idx i = 0; i < subsys.size(); ++i)
694  subsys_dims[i] = dims[subsys[i]];
695 
696  // check the Kraus operators
697  if (Ks.size() == 0)
698  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
699  if (!internal::check_square_mat(Ks[0]))
700  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
701  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
702  throw Exception("qpp::apply()",
704  for (auto&& it : Ks)
705  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
706  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
707  // END EXCEPTION CHECKS
708 
709  cmat result = cmat::Zero(rA.rows(), rA.rows());
710 
711  for (idx i = 0; i < Ks.size(); ++i)
712  result += apply(rA, Ks[i], subsys, dims);
713 
714  return result;
715 }
716 
727 template<typename Derived>
728 cmat apply(const Eigen::MatrixBase<Derived>& A,
729  const std::vector<cmat>& Ks,
730  const std::vector<idx>& subsys,
731  idx d = 2)
732 {
733  const cmat& rA = A.derived();
734 
735  // EXCEPTION CHECKS
736 
737  // check zero sizes
739  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
740 
741  // check valid dims
742  if (d == 0)
743  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
744  // END EXCEPTION CHECKS
745 
746  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
747  std::vector<idx> dims(N, d); // local dimensions vector
748 
749  return apply(rA, Ks, subsys, dims);
750 }
751 
763 inline cmat kraus2super(const std::vector<cmat>& Ks)
764 {
765  // EXCEPTION CHECKS
766 
767  if (Ks.size() == 0)
768  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
769  if (!internal::check_nonzero_size(Ks[0]))
770  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
771  if (!internal::check_square_mat(Ks[0]))
772  throw Exception("qpp::kraus2super()",
774  for (auto&& it : Ks)
775  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
776  throw Exception("qpp::kraus2super()",
778  // END EXCEPTION CHECKS
779 
780  idx D = static_cast<idx>(Ks[0].rows());
781 
782  cmat result(D * D, D * D);
783  cmat MN = cmat::Zero(D, D);
784  bra A = bra::Zero(D);
785  ket B = ket::Zero(D);
786  cmat EMN = cmat::Zero(D, D);
787 
788 #ifdef WITH_OPENMP_
789 #pragma omp parallel for collapse(2)
790 #endif // WITH_OPENMP_
791  for (idx m = 0; m < D; ++m)
792  {
793  for (idx n = 0; n < D; ++n)
794  {
795 #ifdef WITH_OPENMP_
796 #pragma omp critical
797 #endif // WITH_OPENMP_
798  {
799  // compute E(|m><n|)
800  MN(m, n) = 1;
801  for (idx i = 0; i < Ks.size(); ++i)
802  EMN += Ks[i] * MN * adjoint(Ks[i]);
803  MN(m, n) = 0;
804 
805  for (idx a = 0; a < D; ++a)
806  {
807  A(a) = 1;
808  for (idx b = 0; b < D; ++b)
809  {
810  // compute result(ab,mn)=<a|E(|m><n)|b>
811  B(b) = 1;
812  result(a * D + b, m * D + n) =
813  static_cast<cmat>(A * EMN * B).value();
814  B(b) = 0;
815  }
816  A(a) = 0;
817  }
818  EMN = cmat::Zero(D, D);
819  }
820  }
821  }
822 
823  return result;
824 }
825 
841 inline cmat kraus2choi(const std::vector<cmat>& Ks)
842 {
843  // EXCEPTION CHECKS
844 
845  if (Ks.size() == 0)
846  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
847  if (!internal::check_nonzero_size(Ks[0]))
848  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
849  if (!internal::check_square_mat(Ks[0]))
850  throw Exception("qpp::kraus2choi()",
852  for (auto&& it : Ks)
853  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
854  throw Exception("qpp::kraus2choi()",
856  // END EXCEPTION CHECKS
857 
858  idx D = static_cast<idx>(Ks[0].rows());
859 
860  // construct the D x D \sum |jj> vector
861  // (un-normalized maximally entangled state)
862  cmat MES = cmat::Zero(D * D, 1);
863  for (idx a = 0; a < D; ++a)
864  MES(a * D + a) = 1;
865 
866  cmat Omega = MES * adjoint(MES);
867 
868  cmat result = cmat::Zero(D * D, D * D);
869 
870 #ifdef WITH_OPENMP_
871 #pragma omp parallel for
872 #endif // WITH_OPENMP_
873  for (idx i = 0; i < Ks.size(); ++i)
874  {
875 #ifdef WITH_OPENMP_
876 #pragma omp critical
877 #endif // WITH_OPENMP_
878  {
879  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
880  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
881  }
882  }
883 
884  return result;
885 }
886 
900 inline std::vector<cmat> choi2kraus(const cmat& A)
901 {
902  // EXCEPTION CHECKS
903 
905  throw Exception("qpp::choi2kraus()", Exception::Type::ZERO_SIZE);
907  throw Exception("qpp::choi2kraus()",
909  idx D = internal::get_dim_subsys(A.rows(), 2);
910  // check equal dimensions
911  if (D * D != static_cast<idx>(A.rows()))
912  throw Exception("qpp::choi2kraus()", Exception::Type::DIMS_INVALID);
913  // END EXCEPTION CHECKS
914 
915  dmat ev = hevals(A);
916  cmat evec = hevects(A);
917  std::vector<cmat> result;
918 
919  for (idx i = 0; i < D * D; ++i)
920  {
921  if (std::abs(ev(i)) > eps)
922  result.push_back(
923  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
924  }
925 
926  return result;
927 }
928 
936 inline cmat choi2super(const cmat& A)
937 {
938  // EXCEPTION CHECKS
939 
941  throw Exception("qpp::choi2super()", Exception::Type::ZERO_SIZE);
943  throw Exception("qpp::choi2super()",
945  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
946  // check equal dimensions
947  if (D * D != static_cast<idx>(A.rows()))
948  throw Exception("qpp::choi2super()", Exception::Type::DIMS_INVALID);
949  // END EXCEPTION CHECKS
950 
951  cmat result(D * D, D * D);
952 
953 #ifdef WITH_OPENMP_
954 #pragma omp parallel for collapse(4)
955 #endif // WITH_OPENMP_
956  for (idx a = 0; a < D; ++a)
957  for (idx b = 0; b < D; ++b)
958  for (idx m = 0; m < D; ++m)
959  for (idx n = 0; n < D; ++n)
960  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
961 
962  return result;
963 }
964 
972 inline cmat super2choi(const cmat& A)
973 {
974  // EXCEPTION CHECKS
975 
977  throw Exception("qpp::super2choi()", Exception::Type::ZERO_SIZE);
979  throw Exception("qpp::super2choi()",
981  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()),2);
982  // check equal dimensions
983  if (D * D != static_cast<idx>(A.rows()))
984  throw Exception("qpp::super2choi()", Exception::Type::DIMS_INVALID);
985  // END EXCEPTION CHECKS
986 
987  cmat result(D * D, D * D);
988 
989 #ifdef WITH_OPENMP_
990 #pragma omp parallel for collapse(4)
991 #endif // WITH_OPENMP_
992  for (idx a = 0; a < D; ++a)
993  for (idx b = 0; b < D; ++b)
994  for (idx m = 0; m < D; ++m)
995  for (idx n = 0; n < D; ++n)
996  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
997 
998  return result;
999 }
1000 
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 
1115 template<typename Derived>
1116 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1117  idx d = 2)
1118 {
1119  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1120 
1121  // EXCEPTION CHECKS
1122 
1123  // check zero size
1125  throw Exception("qpp::ptrace1()", Exception::Type::ZERO_SIZE);
1126 
1127  // check valid dims
1128  if (d == 0)
1129  throw Exception("qpp::ptrace1()",
1131  // END EXCEPTION CHECKS
1132 
1133  std::vector<idx> dims(2, d); // local dimensions vector
1134 
1135  return ptrace1(A, dims);
1136 }
1137 
1151 template<typename Derived>
1152 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1153  const std::vector<idx>& dims)
1154 {
1155  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1156  = A.derived();
1157 
1158  // EXCEPTION CHECKS
1159 
1160  // check zero-size
1162  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1163 
1164  // check that dims is a valid dimension vector
1165  if (!internal::check_dims(dims))
1166  throw Exception("qpp::ptrace2()", Exception::Type::DIMS_INVALID);
1167 
1168  // check dims has only 2 elements
1169  if (dims.size() != 2)
1170  throw Exception("qpp::ptrace2()", Exception::Type::NOT_BIPARTITE);
1171  // END EXCEPTION CHECKS
1172 
1173  idx DA = dims[0];
1174  idx DB = dims[1];
1175 
1178 
1179  //************ ket ************//
1180  if (internal::check_cvector(rA)) // we have a ket
1181  {
1182  // check that dims match the dimension of A
1183  if (!internal::check_dims_match_cvect(dims, rA))
1184  throw Exception("qpp::ptrace2()",
1186 
1187  auto worker = [=](idx i, idx j) noexcept
1188  -> typename Derived::Scalar
1189  {
1190  typename Derived::Scalar sum = 0;
1191  for (idx m = 0; m < DB; ++m)
1192  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1193 
1194  return sum;
1195  }; /* end worker */
1196 
1197 #ifdef WITH_OPENMP_
1198 #pragma omp parallel for collapse(2)
1199 #endif // WITH_OPENMP_
1200  for (idx j = 0; j < DA; ++j) // column major order for speed
1201  for (idx i = 0; i < DA; ++i)
1202  result(i, j) = worker(i, j);
1203 
1204  return result;
1205  }
1206  //************ density matrix ************//
1207  else if (internal::check_square_mat(rA)) // we have a density operator
1208  {
1209  // check that dims match the dimension of A
1210  if (!internal::check_dims_match_mat(dims, rA))
1211  throw Exception("qpp::ptrace2()",
1213 
1214 #ifdef WITH_OPENMP_
1215 #pragma omp parallel for collapse(2)
1216 #endif // WITH_OPENMP_
1217  for (idx j = 0; j < DA; ++j) // column major order for speed
1218  for (idx i = 0; i < DA; ++i)
1219  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1220 
1221  return result;
1222  }
1223  //************ Exception: not ket nor density matrix ************//
1224  else
1225  throw Exception("qpp::ptrace1()",
1227 }
1228 
1242 template<typename Derived>
1243 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1244  idx d = 2)
1245 {
1246  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1247 
1248  // EXCEPTION CHECKS
1249 
1250  // check zero size
1252  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1253 
1254  // check valid dims
1255  if (d == 0)
1256  throw Exception("qpp::ptrace2()",
1258  // END EXCEPTION CHECKS
1259 
1260  std::vector<idx> dims(2, d); // local dimensions vector
1261 
1262  return ptrace2(A, dims);
1263 }
1264 
1279 template<typename Derived>
1280 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1281  const std::vector<idx>& subsys,
1282  const std::vector<idx>& dims)
1283 {
1284  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1285  = A.derived();
1286 
1287  // EXCEPTION CHECKS
1288 
1289  // check zero-size
1291  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1292 
1293  // check that dims is a valid dimension vector
1294  if (!internal::check_dims(dims))
1295  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1296 
1297  // check that subsys are valid
1298  if (!internal::check_subsys_match_dims(subsys, dims))
1299  throw Exception("qpp::ptrace()",
1301  // END EXCEPTION CHECKS
1302 
1303  idx D = static_cast<idx>(rA.rows());
1304  idx N = dims.size();
1305  idx Nsubsys = subsys.size();
1306  idx Nsubsys_bar = N - Nsubsys;
1307  idx Dsubsys = 1;
1308  for (idx i = 0; i < Nsubsys; ++i)
1309  Dsubsys *= dims[subsys[i]];
1310  idx Dsubsys_bar = D / Dsubsys;
1311 
1312  idx Cdims[maxn];
1313  idx Csubsys[maxn];
1314  idx Cdimssubsys[maxn];
1315  idx Csubsys_bar[maxn];
1316  idx Cdimssubsys_bar[maxn];
1317 
1318  idx Cmidxcolsubsys_bar[maxn];
1319 
1320  std::vector<idx> subsys_bar = complement(subsys, N);
1321  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1322  std::begin(Csubsys_bar));
1323 
1324  for (idx i = 0; i < N; ++i)
1325  {
1326  Cdims[i] = dims[i];
1327  }
1328  for (idx i = 0; i < Nsubsys; ++i)
1329  {
1330  Csubsys[i] = subsys[i];
1331  Cdimssubsys[i] = dims[subsys[i]];
1332  }
1333  for (idx i = 0; i < Nsubsys_bar; ++i)
1334  {
1335  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1336  }
1337 
1339  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1340 
1341  //************ ket ************//
1342  if (internal::check_cvector(rA)) // we have a ket
1343  {
1344  // check that dims match the dimension of A
1345  if (!internal::check_dims_match_cvect(dims, rA))
1346  throw Exception("qpp::ptrace()",
1348 
1349  if (subsys.size() == dims.size())
1350  {
1351  result(0, 0) = (adjoint(rA) * rA).value();
1352  return result;
1353  }
1354 
1355  if (subsys.size() == 0)
1356  return rA * adjoint(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 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  std::conj(rA(internal::multiidx2n(Cmidxcol, N,
1390  Cdims)));
1391  }
1392 
1393  return sm;
1394  }; /* end worker */
1395 
1396  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1397  {
1398  // compute the column multi-indexes of the complement
1399  internal::n2multiidx(j, Nsubsys_bar,
1400  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1401 #ifdef WITH_OPENMP_
1402 #pragma omp parallel for
1403 #endif // WITH_OPENMP_
1404  for (idx i = 0; i < Dsubsys_bar; ++i)
1405  {
1406  result(i, j) = worker(i);
1407  }
1408  }
1409 
1410  return result;
1411  }
1412  //************ density matrix ************//
1413  else if (internal::check_square_mat(rA)) // we have a density operator
1414  {
1415  // check that dims match the dimension of A
1416  if (!internal::check_dims_match_mat(dims, rA))
1417  throw Exception("qpp::ptrace()",
1419 
1420  if (subsys.size() == dims.size())
1421  {
1422  result(0, 0) = rA.trace();
1423  return result;
1424  }
1425 
1426  if (subsys.size() == 0)
1427  return rA;
1428 
1429  auto worker = [=, &Cmidxcolsubsys_bar](idx i) noexcept
1430  -> typename Derived::Scalar
1431  {
1432  // use static allocation for speed!
1433 
1434  idx Cmidxrow[maxn];
1435  idx Cmidxcol[maxn];
1436  idx Cmidxrowsubsys_bar[maxn];
1437  idx Cmidxsubsys[maxn];
1438 
1439  /* get the row/col multi-indexes of the complement */
1440  internal::n2multiidx(i, Nsubsys_bar,
1441  Cdimssubsys_bar, Cmidxrowsubsys_bar);
1442  /* write them in the global row/col multi-indexes */
1443  for (idx k = 0; k < Nsubsys_bar; ++k)
1444  {
1445  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1446  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1447  }
1448  typename Derived::Scalar sm = 0;
1449  for (idx a = 0; a < Dsubsys; ++a)
1450  {
1451  // get the multi-index over which we do the summation
1452  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1453  // write it into the global row/col multi-indexes
1454  for (idx k = 0; k < Nsubsys; ++k)
1455  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1456  = Cmidxsubsys[k];
1457 
1458  // now do the sum
1459  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims),
1460  internal::multiidx2n(Cmidxcol, N, Cdims));
1461  }
1462 
1463  return sm;
1464  }; /* end worker */
1465 
1466  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1467  {
1468  // compute the column multi-indexes of the complement
1469  internal::n2multiidx(j, Nsubsys_bar,
1470  Cdimssubsys_bar, Cmidxcolsubsys_bar);
1471 #ifdef WITH_OPENMP_
1472 #pragma omp parallel for
1473 #endif // WITH_OPENMP_
1474  for (idx i = 0; i < Dsubsys_bar; ++i)
1475  {
1476  result(i, j) = worker(i);
1477  }
1478  }
1479 
1480  return result;
1481  }
1482  //************ Exception: not ket nor density matrix ************//
1483  else
1484  throw Exception("qpp::ptrace()",
1486 }
1487 
1502 template<typename Derived>
1503 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1504  const std::vector<idx>& subsys,
1505  idx d = 2)
1506 {
1507  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1508 
1509  // EXCEPTION CHECKS
1510 
1511  // check zero size
1513  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1514 
1515  // check valid dims
1516  if (d == 0)
1517  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1518  // END EXCEPTION CHECKS
1519 
1520  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1521  std::vector<idx> dims(N, d); // local dimensions vector
1522 
1523  return ptrace(rA, subsys, dims);
1524 }
1525 
1539 template<typename Derived>
1541  const Eigen::MatrixBase<Derived>& A,
1542  const std::vector<idx>& subsys,
1543  const std::vector<idx>& dims)
1544 {
1545  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1546 
1547  // EXCEPTION CHECKS
1548 
1549  // check zero-size
1551  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1552 
1553  // check that dims is a valid dimension vector
1554  if (!internal::check_dims(dims))
1555  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1556 
1557  // check that subsys are valid
1558  if (!internal::check_subsys_match_dims(subsys, dims))
1559  throw Exception("qpp::ptranspose()",
1561  // END EXCEPTION CHECKS
1562 
1563  idx D = static_cast<idx>(rA.rows());
1564  idx N = dims.size();
1565  idx Nsubsys = subsys.size();
1566  idx Cdims[maxn];
1567  idx Cmidxcol[maxn];
1568  idx Csubsys[maxn];
1569 
1570  // copy dims in Cdims and subsys in Csubsys
1571  for (idx i = 0; i < N; ++i)
1572  Cdims[i] = dims[i];
1573  for (idx i = 0; i < Nsubsys; ++i)
1574  Csubsys[i] = subsys[i];
1575 
1576  dyn_mat<typename Derived::Scalar> result(D, D);
1577 
1578  //************ ket ************//
1579  if (internal::check_cvector(rA)) // we have a ket
1580  {
1581  // check that dims match the dimension of A
1582  if (!internal::check_dims_match_cvect(dims, rA))
1583  throw Exception("qpp::ptranspose()",
1585 
1586  if (subsys.size() == dims.size())
1587  return (rA * adjoint(rA)).transpose();
1588 
1589  if (subsys.size() == 0)
1590  return rA * adjoint(rA);
1591 
1592  auto worker = [=, &Cmidxcol](idx i) noexcept
1593  -> 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  std::conj(rA(internal::multiidx2n(midxcoltmp, N,
1611  Cdims)));
1612  }; /* end worker */
1613 
1614  for (idx j = 0; j < D; ++j)
1615  {
1616  // compute the column multi-index
1617  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1618 
1619 #ifdef WITH_OPENMP_
1620 #pragma omp parallel for
1621 #endif // WITH_OPENMP_
1622  for (idx i = 0; i < D; ++i)
1623  result(i, j) = worker(i);
1624  }
1625 
1626  return result;
1627  }
1628  //************ density matrix ************//
1629  else if (internal::check_square_mat(rA)) // we have a density operator
1630  {
1631  // check that dims match the dimension of A
1632  if (!internal::check_dims_match_mat(dims, rA))
1633  throw Exception("qpp::ptranspose()",
1635 
1636  if (subsys.size() == dims.size())
1637  return rA.transpose();
1638 
1639  if (subsys.size() == 0)
1640  return rA;
1641 
1642  auto worker = [=, &Cmidxcol](idx i) noexcept
1643  -> typename Derived::Scalar
1644  {
1645  // use static allocation for speed!
1646  idx midxcoltmp[maxn];
1647  idx midxrow[maxn];
1648 
1649  for (idx k = 0; k < N; ++k)
1650  midxcoltmp[k] = Cmidxcol[k];
1651 
1652  /* compute the row multi-index */
1653  internal::n2multiidx(i, N, Cdims, midxrow);
1654 
1655  for (idx k = 0; k < Nsubsys; ++k)
1656  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1657 
1658  /* writes the result */
1659  return rA(internal::multiidx2n(midxrow, N, Cdims),
1660  internal::multiidx2n(midxcoltmp, N, Cdims));
1661  }; /* end worker */
1662 
1663  for (idx j = 0; j < D; ++j)
1664  {
1665  // compute the column multi-index
1666  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1667 
1668 #ifdef WITH_OPENMP_
1669 #pragma omp parallel for
1670 #endif // WITH_OPENMP_
1671  for (idx i = 0; i < D; ++i)
1672  result(i, j) = worker(i);
1673  }
1674 
1675  return result;
1676  }
1677  //************ Exception: not ket nor density matrix ************//
1678  else
1679  throw Exception("qpp::ptranspose()",
1681 }
1682 
1696 template<typename Derived>
1698  const Eigen::MatrixBase<Derived>& A,
1699  const std::vector<idx>& subsys,
1700  idx d = 2)
1701 {
1702  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1703 
1704  // EXCEPTION CHECKS
1705 
1706  // check zero size
1708  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1709 
1710  // check valid dims
1711  if (d == 0)
1712  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1713  // END EXCEPTION CHECKS
1714 
1715  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1716  std::vector<idx> dims(N, d); // local dimensions vector
1717 
1718  return ptranspose(rA, subsys, dims);
1719 }
1720 
1733 template<typename Derived>
1735  const Eigen::MatrixBase<Derived>& A,
1736  const std::vector<idx>& perm,
1737  const std::vector<idx>& dims)
1738 {
1739  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1740 
1741  // EXCEPTION CHECKS
1742 
1743  // check zero-size
1745  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1746 
1747  // check that dims is a valid dimension vector
1748  if (!internal::check_dims(dims))
1749  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1750 
1751  // check that we have a valid permutation
1752  if (!internal::check_perm(perm))
1753  throw Exception("qpp::syspermute()", Exception::Type::PERM_INVALID);
1754 
1755  // check that permutation match dimensions
1756  if (perm.size() != dims.size())
1757  throw Exception("qpp::syspermute()",
1759  // END EXCEPTION CHECKS
1760 
1761  idx D = static_cast<idx>(rA.rows());
1762  idx N = dims.size();
1763 
1765 
1766  //************ ket ************//
1767  if (internal::check_cvector(rA)) // we have a column vector
1768  {
1769  idx Cdims[maxn];
1770  idx Cperm[maxn];
1771 
1772  // check that dims match the dimension of rA
1773  if (!internal::check_dims_match_cvect(dims, rA))
1774  throw Exception("qpp::syspermute()",
1776 
1777  // copy dims in Cdims and perm in Cperm
1778  for (idx i = 0; i < N; ++i)
1779  {
1780  Cdims[i] = dims[i];
1781  Cperm[i] = perm[i];
1782  }
1783  result.resize(D, 1);
1784 
1785  auto worker = [&Cdims, &Cperm, N](idx i) noexcept
1786  -> idx
1787  {
1788  // use static allocation for speed,
1789  // double the size for matrices reshaped as vectors
1790  idx midx[maxn];
1791  idx midxtmp[maxn];
1792  idx permdims[maxn];
1793 
1794  /* compute the multi-index */
1795  internal::n2multiidx(i, N, Cdims, midx);
1796 
1797  for (idx k = 0; k < N; ++k)
1798  {
1799  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1800  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1801  }
1802 
1803  return internal::multiidx2n(midxtmp, N, permdims);
1804  }; /* end worker */
1805 
1806 #ifdef WITH_OPENMP_
1807 #pragma omp parallel for
1808 #endif // WITH_OPENMP_
1809  for (idx i = 0; i < D; ++i)
1810  result(worker(i)) = rA(i);
1811 
1812  return result;
1813  }
1814  //************ density matrix ************//
1815  else if (internal::check_square_mat(rA)) // we have a density operator
1816  {
1817  idx Cdims[2 * maxn];
1818  idx Cperm[2 * maxn];
1819 
1820  // check that dims match the dimension of rA
1821  if (!internal::check_dims_match_mat(dims, rA))
1822  throw Exception("qpp::syspermute()",
1824 
1825  // copy dims in Cdims and perm in Cperm
1826  for (idx i = 0; i < N; ++i)
1827  {
1828  Cdims[i] = dims[i];
1829  Cdims[i + N] = dims[i];
1830  Cperm[i] = perm[i];
1831  Cperm[i + N] = perm[i] + N;
1832  }
1833  result.resize(D * D, 1);
1834  // map A to a column vector
1836  Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1837  const_cast<typename Derived::Scalar*>(rA.data()), D *
1838  D,
1839  1);
1840 
1841  auto worker = [&Cdims, &Cperm, N](idx i) noexcept
1842  -> idx
1843  {
1844  // use static allocation for speed,
1845  // double the size for matrices reshaped as vectors
1846  idx midx[2 * maxn];
1847  idx midxtmp[2 * maxn];
1848  idx permdims[2 * maxn];
1849 
1850  /* compute the multi-index */
1851  internal::n2multiidx(i, 2 * N, Cdims, midx);
1852 
1853  for (idx k = 0; k < 2 * N; ++k)
1854  {
1855  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1856  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1857  }
1858 
1859  return internal::multiidx2n(midxtmp, 2 * N, permdims);
1860  }; /* end worker */
1861 
1862 #ifdef WITH_OPENMP_
1863 #pragma omp parallel for
1864 #endif // WITH_OPENMP_
1865  for (idx i = 0; i < D * D; ++i)
1866  result(worker(i)) = rA(i);
1867 
1868  return reshape(result, D, D);
1869  }
1870  //************ Exception: not ket nor density matrix ************//
1871  else
1872  throw Exception("qpp::syspermute()",
1874 }
1875 
1888 template<typename Derived>
1890  const Eigen::MatrixBase<Derived>& A,
1891  const std::vector<idx>& perm,
1892  idx d = 2)
1893 {
1894  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1895 
1896  // EXCEPTION CHECKS
1897 
1898  // check zero size
1900  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1901 
1902  // check valid dims
1903  if (d == 0)
1904  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1905  // END EXCEPTION CHECKS
1906 
1907  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1908  std::vector<idx> dims(N, d); // local dimensions vector
1909 
1910  return syspermute(rA, perm, dims);
1911 }
1912 
1913 } /* namespace qpp */
1914 
1915 #endif /* OPERATIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:127
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:74
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:198
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:120
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:900
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:1734
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:156
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:48
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:486
Quantum++ main namespace.
Definition: codes.h:30
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:257
idx get_num_subsys(idx sz, idx d)
Definition: util.h:370
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:387
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:936
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:972
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1814
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:167
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:1152
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
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:378
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:841
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:1280
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:127
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:763
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Applies the controlled-gate A to the part subsys of the multi-partite state vector or density matrix ...
Definition: operations.h:55
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:1540
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:140
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