Quantum++  v0.8.2
C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2015 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 
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 dyn_mat<typename Derived1::Scalar>& rstate = state;
63 
64  // EXCEPTION CHECKS
65  // check types
66  if (!std::is_same<typename Derived1::Scalar,
67  typename Derived2::Scalar>::value)
68  throw Exception("qpp::applyCTRL()", Exception::Type::TYPE_MISMATCH);
69 
70  // check zero sizes
72  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
73 
74  // check zero sizes
75  if (!internal::_check_nonzero_size(rstate))
76  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
77 
78  // check square matrix for the gate
80  throw Exception("qpp::applyCTRL()",
82 
83  // check that all control subsystems have the same dimension
84  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
85  for (idx i = 1; i < ctrl.size(); ++i)
86  if (dims[ctrl[i]] != d)
87  throw Exception("qpp::applyCTRL()",
89 
90  // check that dimension is valid
91  if (!internal::_check_dims(dims))
92  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
93 
94  // check subsys is valid w.r.t. dims
95  if (!internal::_check_subsys_match_dims(subsys, dims))
96  throw Exception("qpp::applyCTRL()",
98 
99  // check that gate matches the dimensions of the subsys
100  std::vector <idx> subsys_dims(subsys.size());
101  for (idx i = 0; i < subsys.size(); ++i)
102  subsys_dims[i] = dims[subsys[i]];
103  if (!internal::_check_dims_match_mat(subsys_dims, rA))
104  throw Exception("qpp::applyCTRL()",
106 
107  std::vector <idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
108  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
109  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
110 
111  // check that ctrl + gate subsystem is valid
112  // with respect to local dimensions
113  if (!internal::_check_subsys_match_dims(ctrlgate, dims))
114  throw Exception("qpp::applyCTRL()",
116 
117  // END EXCEPTION CHECKS
118 
119  // construct the table of A^i and (A^dagger)^i
120  std::vector <dyn_mat<typename Derived1::Scalar>> Ai;
121  std::vector <dyn_mat<typename Derived1::Scalar>> Aidagger;
122  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i)
123  {
124  Ai.push_back(powm(rA, i));
125  Aidagger.push_back(powm(adjoint(rA), i));
126  }
127 
128  idx D = rstate.rows(); // total dimension
129  idx n = dims.size(); // total number of subsystems
130  idx ctrlsize = ctrl.size(); // dimension of ctrl subsystem
131  idx DA = rA.rows(); // dimension of gate subsystem
132 
133  idx Cdims[maxn]; // local dimensions
134  idx CdimsA[maxn]; // local dimensions
135  idx CdimsCTRLAbar[maxn]; // local dimensions
136 
137  // compute the complementary subsystem of ctrlgate w.r.t. dims
138  std::vector <idx> ctrlgatebar = complement(ctrlgate, n);
139 
140  idx DCTRLAbar = 1; // dimension of the rest
141  for (idx i = 0; i < ctrlgatebar.size(); ++i)
142  DCTRLAbar *= dims[ctrlgatebar[i]];
143 
144  for (idx k = 0; k < n; ++k)
145  Cdims[k] = dims[k];
146  for (idx k = 0; k < subsys.size(); ++k)
147  CdimsA[k] = dims[subsys[k]];
148  for (idx k = 0; k < ctrlgatebar.size(); ++k)
149  CdimsCTRLAbar[k] = dims[ctrlgatebar[k]];
150 
151 
152  // worker, computes the coefficient and the index for the ket case
153  // used in #pragma omp parallel for collapse
154  auto coeff_idx_ket = [=](idx _i, idx _m, idx _r) noexcept
155  -> std::pair<typename Derived1::Scalar, idx>
156  {
157  idx indx = 0;
158  typename Derived1::Scalar coeff = 0;
159 
160  idx Cmidx[maxn]; // the total multi-index
161  idx CmidxA[maxn];// the gate part multi-index
162  idx CmidxCTRLAbar[maxn];// the rest multi-index
163 
164  // compute the index
165 
166  // set the CTRL part
167  for (idx k = 0; k < ctrl.size(); ++k)
168  {
169  Cmidx[ctrl[k]] = _i;
170  }
171 
172  // set the rest
173  internal::_n2multiidx(_r, n - ctrlgate.size(),
174  CdimsCTRLAbar, CmidxCTRLAbar);
175  for (idx k = 0; k < n - ctrlgate.size(); ++k)
176  {
177  Cmidx[ctrlgatebar[k]] = CmidxCTRLAbar[k];
178  }
179 
180  // set the A part
181  internal::_n2multiidx(_m, subsys.size(), CdimsA, CmidxA);
182  for (idx k = 0; k < subsys.size(); ++k)
183  {
184  Cmidx[subsys[k]] = CmidxA[k];
185  }
186 
187  // we now got the total index
188  indx = internal::_multiidx2n(Cmidx, n, Cdims);
189 
190  // compute the coefficient
191  for (idx _n = 0; _n < DA; ++_n)
192  {
193  internal::_n2multiidx(_n, subsys.size(), CdimsA, CmidxA);
194  for (idx k = 0; k < subsys.size(); ++k)
195  {
196  Cmidx[subsys[k]] = CmidxA[k];
197  }
198  coeff += Ai[_i](_m, _n) *
199  rstate(internal::_multiidx2n(Cmidx, n, Cdims));
200  }
201 
202  return std::make_pair(coeff, indx);
203  };
204 
205  // worker, computes the coefficient and the index
206  // for the density matrix case
207  // used in #pragma omp parallel for collapse
208  auto coeff_idx_rho = [=](idx _i1, idx _m1,
209  idx _r1, idx _i2, idx _m2, idx _r2) noexcept
210  -> std::tuple<typename Derived1::Scalar, idx, idx>
211  {
212  idx idxrow = 0;
213  idx idxcol = 0;
214  typename Derived1::Scalar coeff = 0;
215 
216  idx Cmidxrow[maxn]; // the total row multi-index
217  idx Cmidxcol[maxn];// the total col multi-index
218  idx CmidxArow[maxn];// the gate part row multi-index
219  idx CmidxAcol[maxn];// the gate part col multi-index
220  idx CmidxCTRLAbarrow[maxn];// the rest row multi-index
221  idx CmidxCTRLAbarcol[maxn];// the rest col multi-index
222 
223  // compute the ket/bra indexes
224 
225  // set the CTRL part
226  for (idx k = 0; k < ctrl.size(); ++k)
227  {
228  Cmidxrow[ctrl[k]] = _i1;
229  Cmidxcol[ctrl[k]] = _i2;
230  }
231 
232  // set the rest
233  internal::_n2multiidx(_r1, n - ctrlgate.size(),
234  CdimsCTRLAbar, CmidxCTRLAbarrow);
235  internal::_n2multiidx(_r2, n - ctrlgate.size(),
236  CdimsCTRLAbar, CmidxCTRLAbarcol);
237  for (idx k = 0; k < n - ctrlgate.size(); ++k)
238  {
239  Cmidxrow[ctrlgatebar[k]] = CmidxCTRLAbarrow[k];
240  Cmidxcol[ctrlgatebar[k]] = CmidxCTRLAbarcol[k];
241  }
242 
243  // set the A part
244  internal::_n2multiidx(_m1, subsys.size(), CdimsA, CmidxArow);
245  internal::_n2multiidx(_m2, subsys.size(), CdimsA, CmidxAcol);
246  for (idx k = 0; k < subsys.size(); ++k)
247  {
248  Cmidxrow[subsys[k]] = CmidxArow[k];
249  Cmidxcol[subsys[k]] = CmidxAcol[k];
250  }
251 
252  // we now got the total row/col indexes
253  idxrow = internal::_multiidx2n(Cmidxrow, n, Cdims);
254  idxcol = internal::_multiidx2n(Cmidxcol, n, Cdims);
255 
256  // compute the coefficient
257  for (idx _n1 = 0; _n1 < DA; ++_n1)
258  {
259  internal::_n2multiidx(_n1, subsys.size(), CdimsA, CmidxArow);
260  for (idx k = 0; k < subsys.size(); ++k)
261  {
262  Cmidxrow[subsys[k]] = CmidxArow[k];
263  }
264  for (idx _n2 = 0; _n2 < DA; ++_n2)
265  {
266  internal::_n2multiidx(_n2, subsys.size(), CdimsA, CmidxAcol);
267  for (idx k = 0; k < subsys.size(); ++k)
268  {
269  Cmidxcol[subsys[k]] = CmidxAcol[k];
270  }
271  coeff += Ai[_i1](_m1, _n1) *
272  rstate(internal::_multiidx2n(Cmidxrow, n, Cdims),
273  internal::_multiidx2n(Cmidxcol, n, Cdims)) *
274  Aidagger[_i2](_n2, _m2);
275  }
276  }
277 
278  return std::make_tuple(coeff, idxrow, idxcol);
279  };
280 
281  //************ ket ************//
282  if (internal::_check_cvector(rstate)) // we have a ket
283  {
284  // check that dims match state vector
285  if (!internal::_check_dims_match_cvect(dims, rstate))
286  throw Exception("qpp::applyCTRL()",
288  if (D == 1)
289  return rstate;
290 
291  dyn_mat<typename Derived1::Scalar> result = rstate;
292 
293 #pragma omp parallel for collapse(2)
294  for (idx m = 0; m < DA; ++m)
295  for (idx r = 0; r < DCTRLAbar; ++r)
296  {
297  if (ctrlsize == 0) // no control
298  {
299  result(coeff_idx_ket(1, m, r).second) =
300  coeff_idx_ket(1, m, r).first;
301  }
302  else
303  for (idx i = 0; i < d; ++i)
304  {
305  result(coeff_idx_ket(i, m, r).second) =
306  coeff_idx_ket(i, m, r).first;
307  }
308  }
309 
310  return result;
311  }
312  //************ density matrix ************//
313  else if (internal::_check_square_mat(rstate)) // we have a density operator
314  {
315  // check that dims match state matrix
316  if (!internal::_check_dims_match_mat(dims, rstate))
317  throw Exception("qpp::applyCTRL()",
319 
320  if (D == 1)
321  return rstate;
322 
323  dyn_mat<typename Derived1::Scalar> result = rstate;
324 
325 #pragma omp parallel for collapse(4)
326  for (idx m1 = 0; m1 < DA; ++m1)
327  for (idx r1 = 0; r1 < DCTRLAbar; ++r1)
328  for (idx m2 = 0; m2 < DA; ++m2)
329  for (idx r2 = 0; r2 < DCTRLAbar; ++r2)
330  if (ctrlsize == 0) // no control
331  {
332  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
333  1, m2, r2);
334  result(std::get<1>(coeff_idxes),
335  std::get<2>(coeff_idxes)) =
336  std::get<0>(coeff_idxes);
337  }
338  else
339  {
340  for (idx i1 = 0; i1 < d; ++i1)
341  for (idx i2 = 0; i2 < d; ++i2)
342  {
343  auto coeff_idxes = coeff_idx_rho(
344  i1, m1, r1,
345  i2, m2, r2);
346  result(std::get<1>(coeff_idxes),
347  std::get<2>(coeff_idxes)) =
348  std::get<0>(coeff_idxes);
349  }
350  }
351 
352  return result;
353  }
354  //************ Exception: not ket nor density matrix ************//
355  else
356  throw Exception("qpp::applyCTRL()",
358 }
359 
375 template<typename Derived1, typename Derived2>
377  const Eigen::MatrixBase<Derived1>& state,
378  const Eigen::MatrixBase<Derived2>& A,
379  const std::vector <idx>& ctrl,
380  const std::vector <idx>& subsys,
381  idx d = 2)
382 {
383  const dyn_mat<typename Derived1::Scalar>& rstate = state;
385 
386  // check zero size
387  if (!internal::_check_nonzero_size(rstate))
388  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
389 
390  idx n =
391  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
392  std::log2(d)));
393  std::vector <idx> dims(n, d); // local dimensions vector
394 
395  return applyCTRL(rstate, rA, ctrl, subsys, dims);
396 }
397 
411 template<typename Derived1, typename Derived2>
413  const Eigen::MatrixBase<Derived1>& state,
414  const Eigen::MatrixBase<Derived2>& A,
415  const std::vector <idx>& subsys,
416  const std::vector <idx>& dims)
417 {
418  const dyn_mat<typename Derived1::Scalar>& rstate = state;
420 
421  // EXCEPTION CHECKS
422 
423  // check types
424  if (!std::is_same<typename Derived1::Scalar,
425  typename Derived2::Scalar>::value)
426  throw Exception("qpp::apply()", Exception::Type::TYPE_MISMATCH);
427 
428  // check zero sizes
430  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
431 
432  // check zero sizes
433  if (!internal::_check_nonzero_size(rstate))
434  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
435 
436  // check square matrix for the gate
438  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
439 
440  // check that dimension is valid
441  if (!internal::_check_dims(dims))
442  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
443 
444  // check subsys is valid w.r.t. dims
445  if (!internal::_check_subsys_match_dims(subsys, dims))
446  throw Exception("qpp::apply()", Exception::Type::SUBSYS_MISMATCH_DIMS);
447 
448  // check that gate matches the dimensions of the subsys
449  std::vector <idx> subsys_dims(subsys.size());
450  for (idx i = 0; i < subsys.size(); ++i)
451  subsys_dims[i] = dims[subsys[i]];
452  if (!internal::_check_dims_match_mat(subsys_dims, rA))
453  throw Exception("qpp::apply()",
455 
456  //************ ket ************//
457  if (internal::_check_cvector(rstate)) // we have a ket
458  {
459  // check that dims match state vector
460  if (!internal::_check_dims_match_cvect(dims, rstate))
461  throw Exception("qpp::apply()",
463 
464  return applyCTRL(rstate, rA, {}, subsys, dims);
465  }
466  //************ density matrix ************//
467  else if (internal::_check_square_mat(rstate)) // we have a density operator
468  {
469 
470  // check that dims match state matrix
471  if (!internal::_check_dims_match_mat(dims, rstate))
472  throw Exception("qpp::apply()",
474 
475  return applyCTRL(rstate, rA, {}, subsys, dims);
476  }
477  //************ Exception: not ket nor density matrix ************//
478  else
479  throw Exception("qpp::apply()",
481 }
482 
496 template<typename Derived1, typename Derived2>
498  const Eigen::MatrixBase<Derived1>& state,
499  const Eigen::MatrixBase<Derived2>& A,
500  const std::vector <idx>& subsys,
501  idx d = 2)
502 {
503  const dyn_mat<typename Derived1::Scalar>& rstate = state;
505 
506  // check zero size
507  if (!internal::_check_nonzero_size(rstate))
508  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
509 
510  idx n =
511  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
512  std::log2(d)));
513  std::vector <idx> dims(n, d); // local dimensions vector
514 
515  return apply(rstate, rA, subsys, dims);
516 }
517 
526 template<typename Derived>
527 cmat apply(const Eigen::MatrixBase<Derived>& rho,
528  const std::vector <cmat>& Ks)
529 {
530  const cmat& rrho = rho;
531 
532  // EXCEPTION CHECKS
534  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
535  if (!internal::_check_square_mat(rrho))
536  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
537  if (Ks.size() == 0)
538  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
539  if (!internal::_check_square_mat(Ks[0]))
540  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
541  if (Ks[0].rows() != rrho.rows())
542  throw Exception("qpp::apply()",
544  for (auto&& it : Ks)
545  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
546  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
547 
548  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
549 
550 #pragma omp parallel for
551  for (idx i = 0; i < Ks.size(); ++i)
552  {
553 #pragma omp critical
554  {
555  result += Ks[i] * rrho * adjoint(Ks[i]);
556  }
557  }
558 
559  return result;
560 }
561 
572 template<typename Derived>
573 cmat apply(const Eigen::MatrixBase<Derived>& rho,
574  const std::vector <cmat>& Ks,
575  const std::vector <idx>& subsys,
576  const std::vector <idx>& dims)
577 {
578  const cmat& rrho = rho;
579 
580  // EXCEPTION CHECKS
581  // check zero sizes
583  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
584 
585  // check square matrix for the rho
586  if (!internal::_check_square_mat(rrho))
587  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
588 
589  // check that dimension is valid
590  if (!internal::_check_dims(dims))
591  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
592 
593  // check that dims match rho matrix
594  if (!internal::_check_dims_match_mat(dims, rrho))
595  throw Exception("qpp::apply()",
597 
598  // check subsys is valid w.r.t. dims
599  if (!internal::_check_subsys_match_dims(subsys, dims))
600  throw Exception("qpp::apply()",
602 
603  std::vector <idx> subsys_dims(subsys.size());
604  for (idx i = 0; i < subsys.size(); ++i)
605  subsys_dims[i] = dims[subsys[i]];
606 
607  // check the Kraus operators
608  if (Ks.size() == 0)
609  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
610  if (!internal::_check_square_mat(Ks[0]))
611  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
612  if (!internal::_check_dims_match_mat(subsys_dims, Ks[0]))
613  throw Exception("qpp::apply()",
615  for (auto&& it : Ks)
616  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
617  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
618 
619  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
620 
621  for (idx i = 0; i < Ks.size(); ++i)
622  result += apply(rrho, Ks[i], subsys, dims);
623 
624  return result;
625 }
626 
637 template<typename Derived>
638 cmat apply(const Eigen::MatrixBase<Derived>& rho,
639  const std::vector <cmat>& Ks,
640  const std::vector <idx>& subsys,
641  idx d = 2)
642 {
643  const cmat& rrho = rho;
644 
645  // check zero sizes
647  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
648 
649  idx n =
650  static_cast<idx>(std::llround(std::log2(rrho.rows()) /
651  std::log2(d)));
652  std::vector <idx> dims(n, d); // local dimensions vector
653 
654  return apply(rrho, Ks, subsys, dims);
655 }
656 
668 inline cmat kraus2super(const std::vector <cmat>& Ks)
669 {
670  // EXCEPTION CHECKS
671  if (Ks.size() == 0)
672  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
673  if (!internal::_check_nonzero_size(Ks[0]))
674  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
675  if (!internal::_check_square_mat(Ks[0]))
676  throw Exception("qpp::kraus2super()",
678  for (auto&& it : Ks)
679  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
680  throw Exception("qpp::kraus2super()",
682  idx D = static_cast<idx>(Ks[0].rows());
683 
684  cmat result(D * D, D * D);
685  cmat MN = cmat::Zero(D, D);
686  bra A = bra::Zero(D);
687  ket B = ket::Zero(D);
688  cmat EMN = cmat::Zero(D, D);
689 
690 #pragma omp parallel for collapse(2)
691  for (idx m = 0; m < D; ++m)
692  {
693  for (idx n = 0; n < D; ++n)
694  {
695 #pragma omp critical
696  {
697  // compute E(|m><n|)
698  MN(m, n) = 1;
699  for (idx i = 0; i < Ks.size(); ++i)
700  EMN += Ks[i] * MN * adjoint(Ks[i]);
701  MN(m, n) = 0;
702 
703  for (idx a = 0; a < D; ++a)
704  {
705  A(a) = 1;
706  for (idx b = 0; b < D; ++b)
707  {
708  // compute result(ab,mn)=<a|E(|m><n)|b>
709  B(b) = 1;
710  result(a * D + b, m * D + n) =
711  static_cast<cmat>(A * EMN * B).value();
712  B(b) = 0;
713  }
714  A(a) = 0;
715  }
716  EMN = cmat::Zero(D, D);
717  }
718  }
719  }
720 
721  return result;
722 }
723 
739 inline cmat kraus2choi(const std::vector <cmat>& Ks)
740 {
741  // EXCEPTION CHECKS
742  if (Ks.size() == 0)
743  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
744  if (!internal::_check_nonzero_size(Ks[0]))
745  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
746  if (!internal::_check_square_mat(Ks[0]))
747  throw Exception("qpp::kraus2choi()",
749  for (auto&& it : Ks)
750  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
751  throw Exception("qpp::kraus2choi()",
753  idx D = static_cast<idx>(Ks[0].rows());
754 
755  // construct the D x D \sum |jj> vector
756  // (un-normalized maximally entangled state)
757  cmat MES = cmat::Zero(D * D, 1);
758  for (idx a = 0; a < D; ++a)
759  MES(a * D + a) = 1;
760 
761  cmat Omega = MES * adjoint(MES);
762 
763  cmat result = cmat::Zero(D * D, D * D);
764 
765 #pragma omp parallel for
766  for (idx i = 0; i < Ks.size(); ++i)
767  {
768 #pragma omp critical
769  {
770  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
771  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
772  }
773  }
774 
775  return result;
776 }
777 
791 inline std::vector <cmat> choi2kraus(const cmat& A)
792 {
793  // EXCEPTION CHECKS
795  throw Exception("qpp::choi2kraus()", Exception::Type::ZERO_SIZE);
797  throw Exception("qpp::choi2kraus()",
799  idx D = static_cast<idx>(std::llround(
800  std::sqrt(static_cast<double>(A.rows()))));
801  if (D * D != static_cast<idx>(A.rows()))
802  throw Exception("qpp::choi2kraus()", Exception::Type::DIMS_INVALID);
803 
804  dmat ev = hevals(A);
805  cmat evec = hevects(A);
806  std::vector <cmat> result;
807 
808  for (idx i = 0; i < D * D; ++i)
809  {
810  if (std::abs(ev(i)) > eps)
811  result.push_back(
812  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
813  }
814 
815  return result;
816 }
817 
825 inline cmat choi2super(const cmat& A)
826 {
827  // EXCEPTION CHECKS
829  throw Exception("qpp::choi2super()", Exception::Type::ZERO_SIZE);
831  throw Exception("qpp::choi2super()",
833  idx D = static_cast<idx>(std::llround(
834  std::sqrt(static_cast<double>(A.rows()))));
835  if (D * D != static_cast<idx>(A.rows()))
836  throw Exception("qpp::choi2super()", Exception::Type::DIMS_INVALID);
837 
838  cmat result(D * D, D * D);
839 
840 #pragma omp parallel for collapse(4)
841  for (idx a = 0; a < D; ++a)
842  for (idx b = 0; b < D; ++b)
843  for (idx m = 0; m < D; ++m)
844  for (idx n = 0; n < D; ++n)
845  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
846 
847  return result;
848 }
849 
857 inline cmat super2choi(const cmat& A)
858 {
859  // EXCEPTION CHECKS
861  throw Exception("qpp::super2choi()", Exception::Type::ZERO_SIZE);
863  throw Exception("qpp::super2choi()",
865  idx D = static_cast<idx>(std::llround(
866  std::sqrt(static_cast<double>(A.rows()))));
867  if (D * D != static_cast<idx>(A.rows()))
868  throw Exception("qpp::super2choi()", Exception::Type::DIMS_INVALID);
869 
870  cmat result(D * D, D * D);
871 
872 #pragma omp parallel for collapse(4)
873  for (idx a = 0; a < D; ++a)
874  for (idx b = 0; b < D; ++b)
875  for (idx m = 0; m < D; ++m)
876  for (idx n = 0; n < D; ++n)
877  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
878 
879  return result;
880 }
881 
896 template<typename Derived>
897 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
898  const std::vector <idx>& dims)
899 {
900  const dyn_mat<typename Derived::Scalar>& rA = A;
901 
902  // Error checks
903 
904  // check zero-size
906  throw Exception("qpp::ptrace1()", Exception::Type::ZERO_SIZE);
907 
908  // check that dims is a valid dimension vector
909  if (!internal::_check_dims(dims))
910  throw Exception("qpp::ptrace1()", Exception::Type::DIMS_INVALID);
911 
912  // check dims has only 2 elements
913  if (dims.size() != 2)
914  throw Exception("qpp::ptrace1()", Exception::Type::NOT_BIPARTITE);
915 
916  idx DA = dims[0];
917  idx DB = dims[1];
918 
921 
922  //************ ket ************//
923  if (internal::_check_cvector(rA)) // we have a ket
924  {
925  // check that dims match the dimension of A
926  if (!internal::_check_dims_match_cvect(dims, rA))
927  throw Exception("qpp::ptrace1()",
929 
930  auto worker = [=](idx i, idx j) noexcept
931  -> typename Derived::Scalar
932  {
933  typename Derived::Scalar sum = 0;
934  for (idx m = 0; m < DA; ++m)
935  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
936 
937  return sum;
938  };
939 
940 #pragma omp parallel for collapse(2)
941  for (idx j = 0; j < DB; ++j) // column major order for speed
942  for (idx i = 0; i < DB; ++i)
943  result(i, j) = worker(i, j);
944 
945  return result;
946  }
947  //************ density matrix ************//
948  else if (internal::_check_square_mat(rA)) // we have a density operator
949  {
950  // check that dims match the dimension of A
951  if (!internal::_check_dims_match_mat(dims, rA))
952  throw Exception("qpp::ptrace1()",
954 
955  auto worker = [=](idx i, idx j) noexcept
956  -> typename Derived::Scalar
957  {
958  typename Derived::Scalar sum = 0;
959  for (idx m = 0; m < DA; ++m)
960  sum += rA(m * DB + i, m * DB + j);
961 
962  return sum;
963  };
964 
965 #pragma omp parallel for collapse(2)
966  for (idx j = 0; j < DB; ++j) // column major order for speed
967  for (idx i = 0; i < DB; ++i)
968  result(i, j) = worker(i, j);
969 
970  return result;
971  }
972  //************ Exception: not ket nor density matrix ************//
973  else
974  throw Exception("qpp::ptrace1()",
976 }
977 
992 template<typename Derived>
993 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
994  const std::vector <idx>& dims)
995 {
996  const dyn_mat<typename Derived::Scalar>& rA = A;
997 
998  // Error checks
999 
1000  // check zero-size
1002  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1003 
1004  // check that dims is a valid dimension vector
1005  if (!internal::_check_dims(dims))
1006  throw Exception("qpp::ptrace2()", Exception::Type::DIMS_INVALID);
1007 
1008  // check dims has only 2 elements
1009  if (dims.size() != 2)
1010  throw Exception("qpp::ptrace2()", Exception::Type::NOT_BIPARTITE);
1011 
1012  idx DA = dims[0];
1013  idx DB = dims[1];
1014 
1016  dyn_mat < typename Derived::Scalar > ::Zero(DA, DA);
1017 
1018  //************ ket ************//
1019  if (internal::_check_cvector(rA)) // we have a ket
1020  {
1021  // check that dims match the dimension of A
1022  if (!internal::_check_dims_match_cvect(dims, rA))
1023  throw Exception("qpp::ptrace2()",
1025 
1026  auto worker = [=](idx i, idx j) noexcept
1027  -> typename Derived::Scalar
1028  {
1029  typename Derived::Scalar sum = 0;
1030  for (idx m = 0; m < DB; ++m)
1031  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1032 
1033  return sum;
1034  };
1035 
1036 #pragma omp parallel for collapse(2)
1037  for (idx j = 0; j < DA; ++j) // column major order for speed
1038  for (idx i = 0; i < DA; ++i)
1039  result(i, j) = worker(i, j);
1040 
1041  return result;
1042  }
1043  //************ density matrix ************//
1044  else if (internal::_check_square_mat(rA)) // we have a density operator
1045  {
1046  // check that dims match the dimension of A
1047  if (!internal::_check_dims_match_mat(dims, rA))
1048  throw Exception("qpp::ptrace2()",
1050 
1051 #pragma omp parallel for collapse(2)
1052  for (idx j = 0; j < DA; ++j) // column major order for speed
1053  for (idx i = 0; i < DA; ++i)
1054  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1055 
1056  return result;
1057  }
1058  //************ Exception: not ket nor density matrix ************//
1059  else
1060  throw Exception("qpp::ptrace1()",
1062 }
1063 
1078 template<typename Derived>
1079 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1080  const std::vector <idx>& subsys,
1081  const std::vector <idx>& dims)
1082 {
1083  const dyn_mat<typename Derived::Scalar>& rA = A;
1084 
1085  // error checks
1086 
1087  // check zero-size
1089  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1090 
1091  // check that dims is a valid dimension vector
1092  if (!internal::_check_dims(dims))
1093  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1094 
1095  // check that subsys are valid
1096  if (!internal::_check_subsys_match_dims(subsys, dims))
1097  throw Exception("qpp::ptrace()",
1099 
1100  idx D = static_cast<idx>(rA.rows());
1101  idx n = dims.size();
1102  idx nsubsys = subsys.size();
1103  idx nsubsysbar = n - nsubsys;
1104  idx dimsubsys = 1;
1105  for (idx i = 0; i < nsubsys; ++i)
1106  dimsubsys *= dims[subsys[i]];
1107  idx dimsubsysbar = D / dimsubsys;
1108 
1109  idx Cdims[maxn];
1110  idx Csubsys[maxn];
1111  idx Cdimssubsys[maxn];
1112  idx Csubsysbar[maxn];
1113  idx Cdimssubsysbar[maxn];
1114 
1115  idx Cmidxcolsubsysbar[maxn];
1116 
1117  std::vector <idx> subsys_bar = complement(subsys, n);
1118  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1119  std::begin(Csubsysbar));
1120 
1121  for (idx i = 0; i < n; ++i)
1122  {
1123  Cdims[i] = dims[i];
1124  }
1125  for (idx i = 0; i < nsubsys; ++i)
1126  {
1127  Csubsys[i] = subsys[i];
1128  Cdimssubsys[i] = dims[subsys[i]];
1129  }
1130  for (idx i = 0; i < nsubsysbar; ++i)
1131  {
1132  Cdimssubsysbar[i] = dims[subsys_bar[i]];
1133  }
1134 
1136  dyn_mat < typename Derived::Scalar > (dimsubsysbar, dimsubsysbar);
1137 
1138  //************ ket ************//
1139  if (internal::_check_cvector(rA)) // we have a ket
1140  {
1141  // check that dims match the dimension of A
1142  if (!internal::_check_dims_match_cvect(dims, rA))
1143  throw Exception("qpp::ptrace()",
1145 
1146  if (subsys.size() == dims.size())
1147  {
1148  result(0, 0) = (adjoint(rA) * rA).value();
1149  return result;
1150  }
1151 
1152  if (subsys.size() == 0)
1153  return rA * adjoint(rA);
1154 
1155  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1156  -> typename Derived::Scalar
1157  {
1158  // use static allocation for speed!
1159 
1160  idx Cmidxrow[maxn];
1161  idx Cmidxcol[maxn];
1162  idx Cmidxrowsubsysbar[maxn];
1163  idx Cmidxsubsys[maxn];
1164 
1165  /* get the row multi-indexes of the complement */
1166  internal::_n2multiidx(i, nsubsysbar,
1167  Cdimssubsysbar, Cmidxrowsubsysbar);
1168  /* write them in the global row/col multi-indexes */
1169  for (idx k = 0; k < nsubsysbar; ++k)
1170  {
1171  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1172  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1173  }
1174  typename Derived::Scalar sm = 0;
1175  for (idx a = 0; a < dimsubsys; ++a)
1176  {
1177  // get the multi-index over which we do the summation
1178  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1179  // write it into the global row/col multi-indexes
1180  for (idx k = 0; k < nsubsys; ++k)
1181  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1182  = Cmidxsubsys[k];
1183 
1184  // now do the sum
1185  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims)) *
1186  std::conj(rA(internal::_multiidx2n(Cmidxcol, n,
1187  Cdims)));
1188  }
1189 
1190  return sm;
1191  };
1192 
1193  for (idx j = 0; j < dimsubsysbar; ++j) // column major order for speed
1194  {
1195  // compute the column multi-indexes of the complement
1196  internal::_n2multiidx(j, nsubsysbar,
1197  Cdimssubsysbar, Cmidxcolsubsysbar);
1198 #pragma omp parallel for
1199  for (idx i = 0; i < dimsubsysbar; ++i)
1200  {
1201  result(i, j) = worker(i);
1202  }
1203  }
1204 
1205  return result;
1206  }
1207  //************ density matrix ************//
1208  else if (internal::_check_square_mat(rA)) // we have a density operator
1209  {
1210  // check that dims match the dimension of A
1211  if (!internal::_check_dims_match_mat(dims, rA))
1212  throw Exception("qpp::ptrace()",
1214 
1215  if (subsys.size() == dims.size())
1216  {
1217  result(0, 0) = rA.trace();
1218  return result;
1219  }
1220 
1221  if (subsys.size() == 0)
1222  return rA;
1223 
1224  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1225  -> typename Derived::Scalar
1226  {
1227  // use static allocation for speed!
1228 
1229  idx Cmidxrow[maxn];
1230  idx Cmidxcol[maxn];
1231  idx Cmidxrowsubsysbar[maxn];
1232  idx Cmidxsubsys[maxn];
1233 
1234  /* get the row/col multi-indexes of the complement */
1235  internal::_n2multiidx(i, nsubsysbar,
1236  Cdimssubsysbar, Cmidxrowsubsysbar);
1237  /* write them in the global row/col multi-indexes */
1238  for (idx k = 0; k < nsubsysbar; ++k)
1239  {
1240  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1241  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1242  }
1243  typename Derived::Scalar sm = 0;
1244  for (idx a = 0; a < dimsubsys; ++a)
1245  {
1246  // get the multi-index over which we do the summation
1247  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1248  // write it into the global row/col multi-indexes
1249  for (idx k = 0; k < nsubsys; ++k)
1250  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1251  = Cmidxsubsys[k];
1252 
1253  // now do the sum
1254  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims),
1255  internal::_multiidx2n(Cmidxcol, n, Cdims));
1256  }
1257 
1258  return sm;
1259  };
1260 
1261  for (idx j = 0; j < dimsubsysbar; ++j) // column major order for speed
1262  {
1263  // compute the column multi-indexes of the complement
1264  internal::_n2multiidx(j, nsubsysbar,
1265  Cdimssubsysbar, Cmidxcolsubsysbar);
1266 #pragma omp parallel for
1267  for (idx i = 0; i < dimsubsysbar; ++i)
1268  {
1269  result(i, j) = worker(i);
1270  }
1271  }
1272 
1273  return result;
1274  }
1275  //************ Exception: not ket nor density matrix ************//
1276  else
1277  throw Exception("qpp::ptrace()",
1279 }
1280 
1295 template<typename Derived>
1296 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1297  const std::vector <idx>& subsys,
1298  idx d = 2)
1299 {
1300  const dyn_mat<typename Derived::Scalar>& rA = A;
1301 
1302  // check zero size
1304  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1305 
1306  idx n =
1307  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1308  std::log2(d)));
1309  std::vector <idx> dims(n, d); // local dimensions vector
1310 
1311  return ptrace(rA, subsys, dims);
1312 }
1313 
1327 template<typename Derived>
1329  const Eigen::MatrixBase<Derived>& A,
1330  const std::vector <idx>& subsys,
1331  const std::vector <idx>& dims)
1332 {
1333  const dyn_mat<typename Derived::Scalar>& rA = A;
1334 
1335  // error checks
1336 
1337  // check zero-size
1339  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1340 
1341  // check that dims is a valid dimension vector
1342  if (!internal::_check_dims(dims))
1343  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1344 
1345  // check that subsys are valid
1346  if (!internal::_check_subsys_match_dims(subsys, dims))
1347  throw Exception("qpp::ptranspose()",
1349 
1350  idx D = static_cast<idx>(rA.rows());
1351  idx numdims = dims.size();
1352  idx numsubsys = subsys.size();
1353  idx Cdims[maxn];
1354  idx Cmidxcol[maxn];
1355  idx Csubsys[maxn];
1356 
1357  // copy dims in Cdims and subsys in Csubsys
1358  for (idx i = 0; i < numdims; ++i)
1359  Cdims[i] = dims[i];
1360  for (idx i = 0; i < numsubsys; ++i)
1361  Csubsys[i] = subsys[i];
1362 
1363  dyn_mat<typename Derived::Scalar> result(D, D);
1364 
1365  //************ ket ************//
1366  if (internal::_check_cvector(rA)) // we have a ket
1367  {
1368  // check that dims match the dimension of A
1369  if (!internal::_check_dims_match_cvect(dims, rA))
1370  throw Exception("qpp::ptranspose()",
1372 
1373  if (subsys.size() == dims.size())
1374  return (rA * adjoint(rA)).transpose();
1375 
1376  if (subsys.size() == 0)
1377  return rA * adjoint(rA);
1378 
1379  auto worker = [=, &Cmidxcol](idx i) noexcept
1380  -> typename Derived::Scalar
1381  {
1382  // use static allocation for speed!
1383  idx midxcoltmp[maxn];
1384  idx midxrow[maxn];
1385 
1386  for (idx k = 0; k < numdims; ++k)
1387  midxcoltmp[k] = Cmidxcol[k];
1388 
1389  /* compute the row multi-index */
1390  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1391 
1392  for (idx k = 0; k < numsubsys; ++k)
1393  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1394 
1395  /* writes the result */
1396  return rA(internal::_multiidx2n(midxrow, numdims, Cdims)) *
1397  std::conj(rA(internal::_multiidx2n(midxcoltmp, numdims,
1398  Cdims)));
1399  };
1400 
1401  for (idx j = 0; j < D; ++j)
1402  {
1403  // compute the column multi-index
1404  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1405 #pragma omp parallel for
1406  for (idx i = 0; i < D; ++i)
1407  result(i, j) = worker(i);
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::ptranspose()",
1419 
1420  if (subsys.size() == dims.size())
1421  return rA.transpose();
1422 
1423  if (subsys.size() == 0)
1424  return rA;
1425 
1426  auto worker = [=, &Cmidxcol](idx i) noexcept
1427  -> typename Derived::Scalar
1428  {
1429  // use static allocation for speed!
1430  idx midxcoltmp[maxn];
1431  idx midxrow[maxn];
1432 
1433  for (idx k = 0; k < numdims; ++k)
1434  midxcoltmp[k] = Cmidxcol[k];
1435 
1436  /* compute the row multi-index */
1437  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1438 
1439  for (idx k = 0; k < numsubsys; ++k)
1440  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1441 
1442  /* writes the result */
1443  return rA(internal::_multiidx2n(midxrow, numdims, Cdims),
1444  internal::_multiidx2n(midxcoltmp, numdims, Cdims));
1445  };
1446 
1447  for (idx j = 0; j < D; ++j)
1448  {
1449  // compute the column multi-index
1450  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1451 #pragma omp parallel for
1452  for (idx i = 0; i < D; ++i)
1453  result(i, j) = worker(i);
1454  }
1455 
1456  return result;
1457  }
1458  //************ Exception: not ket nor density matrix ************//
1459  else
1460  throw Exception("qpp::ptranspose()",
1462 }
1463 
1477 template<typename Derived>
1479  const Eigen::MatrixBase<Derived>& A,
1480  const std::vector <idx>& subsys,
1481  idx d = 2)
1482 {
1483  const dyn_mat<typename Derived::Scalar>& rA = A;
1484 
1485  // check zero size
1487  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1488 
1489  idx n =
1490  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1491  std::log2(d)));
1492  std::vector <idx> dims(n, d); // local dimensions vector
1493 
1494  return ptranspose(rA, subsys, dims);
1495 }
1496 
1509 template<typename Derived>
1511  const Eigen::MatrixBase<Derived>& A,
1512  const std::vector <idx>& perm,
1513  const std::vector <idx>& dims)
1514 {
1515  const dyn_mat<typename Derived::Scalar>& rA = A;
1516 
1517  // Error checks
1518 
1519  // check zero-size
1521  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1522 
1523  // check that dims is a valid dimension vector
1524  if (!internal::_check_dims(dims))
1525  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1526 
1527  // check that we have a valid permutation
1528  if (!internal::_check_perm(perm))
1529  throw Exception("qpp::syspermute()", Exception::Type::PERM_INVALID);
1530 
1531  // check that permutation match dimensions
1532  if (perm.size() != dims.size())
1533  throw Exception("qpp::syspermute()",
1535 
1536  idx D = static_cast<idx>(rA.rows());
1537  idx numdims = dims.size();
1538 
1540 
1541  //************ ket ************//
1542  if (internal::_check_cvector(rA)) // we have a column vector
1543  {
1544  idx Cdims[maxn];
1545  idx Cperm[maxn];
1546 
1547  // check that dims match the dimension of rA
1548  if (!internal::_check_dims_match_cvect(dims, rA))
1549  throw Exception("qpp::syspermute()",
1551 
1552  // copy dims in Cdims and perm in Cperm
1553  for (idx i = 0; i < numdims; ++i)
1554  {
1555  Cdims[i] = dims[i];
1556  Cperm[i] = perm[i];
1557  }
1558  result.resize(D, 1);
1559 
1560  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1561  -> idx
1562  {
1563  // use static allocation for speed,
1564  // double the size for matrices reshaped as vectors
1565  idx midx[maxn];
1566  idx midxtmp[maxn];
1567  idx permdims[maxn];
1568 
1569  /* compute the multi-index */
1570  internal::_n2multiidx(i, numdims, Cdims, midx);
1571 
1572  for (idx k = 0; k < numdims; ++k)
1573  {
1574  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1575  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1576  }
1577 
1578  return internal::_multiidx2n(midxtmp, numdims, permdims);
1579  };
1580 
1581 #pragma omp parallel for
1582  for (idx i = 0; i < D; ++i)
1583  result(worker(i)) = rA(i);
1584 
1585  return result;
1586  }
1587  //************ density matrix ************//
1588  else if (internal::_check_square_mat(rA)) // we have a density operator
1589  {
1590  idx Cdims[2 * maxn];
1591  idx Cperm[2 * maxn];
1592 
1593  // check that dims match the dimension of rA
1594  if (!internal::_check_dims_match_mat(dims, rA))
1595  throw Exception("qpp::syspermute()",
1597 
1598  // copy dims in Cdims and perm in Cperm
1599  for (idx i = 0; i < numdims; ++i)
1600  {
1601  Cdims[i] = dims[i];
1602  Cdims[i + numdims] = dims[i];
1603  Cperm[i] = perm[i];
1604  Cperm[i + numdims] = perm[i] + numdims;
1605  }
1606  result.resize(D * D, 1);
1607  // map A to a column vector
1609  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1610  const_cast<typename Derived::Scalar*>(rA.data()), D * D,
1611  1);
1612 
1613  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1614  -> idx
1615  {
1616  // use static allocation for speed,
1617  // double the size for matrices reshaped as vectors
1618  idx midx[2 * maxn];
1619  idx midxtmp[2 * maxn];
1620  idx permdims[2 * maxn];
1621 
1622  /* compute the multi-index */
1623  internal::_n2multiidx(i, 2 * numdims, Cdims, midx);
1624 
1625  for (idx k = 0; k < 2 * numdims; ++k)
1626  {
1627  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1628  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1629  }
1630 
1631  return internal::_multiidx2n(midxtmp, 2 * numdims, permdims);
1632  };
1633 
1634 #pragma omp parallel for
1635  for (idx i = 0; i < D * D; ++i)
1636  result(worker(i)) = rA(i);
1637 
1638  return reshape(result, D, D);
1639  }
1640  //************ Exception: not ket nor density matrix ************//
1641  else
1642  throw Exception("qpp::syspermute()",
1644 }
1645 
1658 template<typename Derived>
1660  const Eigen::MatrixBase<Derived>& A,
1661  const std::vector <idx>& perm,
1662  idx d = 2)
1663 {
1664  const dyn_mat<typename Derived::Scalar>& rA = A;
1665 
1666  // check zero size
1668  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1669 
1670  idx n =
1671  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1672  std::log2(d)));
1673  std::vector <idx> dims(n, d); // local dimensions vector
1674 
1675  return syspermute(rA, perm, dims);
1676 }
1677 
1678 } /* namespace qpp */
1679 
1680 #endif /* OPERATIONS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:80
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:791
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:190
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:142
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:71
bool _check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &V)
Definition: util.h:155
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:73
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:61
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:1510
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:83
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:56
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:412
Quantum++ main namespace.
Definition: codes.h:30
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:352
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:825
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:857
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1663
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:897
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:798
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:993
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:257
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:54
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)
Matrix power.
Definition: functions.h:701
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:375
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:739
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:1079
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:121
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:668
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:1328
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:66
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:41
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:126
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:194
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1026