Quantum++  v0.8.6
C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef OPERATIONS_H_
28 #define OPERATIONS_H_
29 
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 
66  // check types
67  if (!std::is_same<typename Derived1::Scalar,
68  typename Derived2::Scalar>::value)
69  throw Exception("qpp::applyCTRL()", Exception::Type::TYPE_MISMATCH);
70 
71  // check zero sizes
73  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
74 
75  // check zero sizes
76  if (!internal::_check_nonzero_size(rstate))
77  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
78 
79  // check square matrix for the gate
81  throw Exception("qpp::applyCTRL()",
83 
84  // check that all control subsystems have the same dimension
85  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
86  for (idx i = 1; i < ctrl.size(); ++i)
87  if (dims[ctrl[i]] != d)
88  throw Exception("qpp::applyCTRL()",
90 
91  // check that dimension is valid
92  if (!internal::_check_dims(dims))
93  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
94 
95  // check subsys is valid w.r.t. dims
96  if (!internal::_check_subsys_match_dims(subsys, dims))
97  throw Exception("qpp::applyCTRL()",
99 
100  // check that gate matches the dimensions of the subsys
101  std::vector<idx> subsys_dims(subsys.size());
102  for (idx i = 0; i < subsys.size(); ++i)
103  subsys_dims[i] = dims[subsys[i]];
104  if (!internal::_check_dims_match_mat(subsys_dims, rA))
105  throw Exception("qpp::applyCTRL()",
107 
108  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
109  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
110  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
111 
112  // check that ctrl + gate subsystem is valid
113  // with respect to local dimensions
114  if (!internal::_check_subsys_match_dims(ctrlgate, dims))
115  throw Exception("qpp::applyCTRL()",
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 = static_cast<idx>(rstate.rows()); // total dimension
129  idx n = dims.size(); // total number of subsystems
130  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
131  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
132  idx subsyssize = subsys.size(); // number of subsystems of the target
133  // dimension of ctrl subsystem
134  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
135  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
136 
137  idx Cdims[maxn]; // local dimensions
138  idx CdimsA[maxn]; // local dimensions
139  idx CdimsCTRL[maxn]; // local dimensions
140  idx CdimsCTRLAbar[maxn]; // local dimensions
141 
142  // compute the complementary subsystem of ctrlgate w.r.t. dims
143  std::vector<idx> ctrlgatebar = complement(ctrlgate, n);
144  // number of subsystems that are complementary to the ctrl+gate
145  idx ctrlgatebarsize = ctrlgatebar.size();
146 
147  idx DCTRLAbar = 1; // dimension of the rest
148  for (idx i = 0; i < ctrlgatebarsize; ++i)
149  DCTRLAbar *= dims[ctrlgatebar[i]];
150 
151  for (idx k = 0; k < n; ++k)
152  Cdims[k] = dims[k];
153  for (idx k = 0; k < subsyssize; ++k)
154  CdimsA[k] = dims[subsys[k]];
155  for (idx k = 0; k < ctrlsize; ++k)
156  CdimsCTRL[k] = d;
157  for (idx k = 0; k < ctrlgatebarsize; ++k)
158  CdimsCTRLAbar[k] = dims[ctrlgatebar[k]];
159 
160  // worker, computes the coefficient and the index for the ket case
161  // used in #pragma omp parallel for collapse
162  auto coeff_idx_ket = [=](idx _i, idx _m, idx _r) noexcept
163  -> std::pair<typename Derived1::Scalar, idx>
164  {
165  idx indx = 0;
166  typename Derived1::Scalar coeff = 0;
167 
168  idx Cmidx[maxn]; // the total multi-index
169  idx CmidxA[maxn]; // the gate part multi-index
170  idx CmidxCTRLAbar[maxn]; // the rest multi-index
171 
172  // compute the index
173 
174  // set the CTRL part
175  for (idx k = 0; k < ctrlsize; ++k)
176  {
177  Cmidx[ctrl[k]] = _i;
178  }
179 
180  // set the rest
181  internal::_n2multiidx(_r, n - ctrlgatesize,
182  CdimsCTRLAbar, CmidxCTRLAbar);
183  for (idx k = 0; k < n - ctrlgatesize; ++k)
184  {
185  Cmidx[ctrlgatebar[k]] = CmidxCTRLAbar[k];
186  }
187 
188  // set the A part
189  internal::_n2multiidx(_m, subsyssize, CdimsA, CmidxA);
190  for (idx k = 0; k < subsyssize; ++k)
191  {
192  Cmidx[subsys[k]] = CmidxA[k];
193  }
194 
195  // we now got the total index
196  indx = internal::_multiidx2n(Cmidx, n, Cdims);
197 
198  // compute the coefficient
199  for (idx _n = 0; _n < DA; ++_n)
200  {
201  internal::_n2multiidx(_n, subsyssize, CdimsA, CmidxA);
202  for (idx k = 0; k < subsyssize; ++k)
203  {
204  Cmidx[subsys[k]] = CmidxA[k];
205  }
206  coeff += Ai[_i](_m, _n) *
207  rstate(internal::_multiidx2n(Cmidx, n, Cdims));
208  }
209 
210  return std::make_pair(coeff, indx);
211  }; /* end coeff_idx_ket */
212 
213  // worker, computes the coefficient and the index
214  // for the density matrix case
215  // used in #pragma omp parallel for collapse
216  auto coeff_idx_rho = [&](idx _i1, idx _m1,
217  idx _r1, idx _i2, idx _m2, idx _r2) noexcept
218  -> std::tuple<typename Derived1::Scalar, idx, idx>
219  {
220  idx idxrow = 0;
221  idx idxcol = 0;
222  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
223 
224  idx Cmidxrow[maxn]; // the total row multi-index
225  idx Cmidxcol[maxn]; // the total col multi-index
226  idx CmidxArow[maxn]; // the gate part row multi-index
227  idx CmidxAcol[maxn]; // the gate part col multi-index
228  idx CmidxCTRLrow[maxn]; // the control row multi-index
229  idx CmidxCTRLcol[maxn]; // the control col multi-index
230  idx CmidxCTRLAbarrow[maxn]; // the rest row multi-index
231  idx CmidxCTRLAbarcol[maxn]; // the rest col multi-index
232 
233  // compute the ket/bra indexes
234 
235  // set the CTRL part
236  internal::_n2multiidx(_i1, ctrlsize,
237  CdimsCTRL, CmidxCTRLrow);
238  internal::_n2multiidx(_i2, ctrlsize,
239  CdimsCTRL, CmidxCTRLcol);
240 
241  for (idx k = 0; k < ctrlsize; ++k)
242  {
243  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
244  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
245  }
246 
247  // set the rest
248  internal::_n2multiidx(_r1, n - ctrlgatesize,
249  CdimsCTRLAbar, CmidxCTRLAbarrow);
250  internal::_n2multiidx(_r2, n - ctrlgatesize,
251  CdimsCTRLAbar, CmidxCTRLAbarcol);
252  for (idx k = 0; k < n - ctrlgatesize; ++k)
253  {
254  Cmidxrow[ctrlgatebar[k]] = CmidxCTRLAbarrow[k];
255  Cmidxcol[ctrlgatebar[k]] = CmidxCTRLAbarcol[k];
256  }
257 
258  // set the A part
259  internal::_n2multiidx(_m1, subsyssize, CdimsA, CmidxArow);
260  internal::_n2multiidx(_m2, subsyssize, CdimsA, CmidxAcol);
261  for (idx k = 0; k < subsys.size(); ++k)
262  {
263  Cmidxrow[subsys[k]] = CmidxArow[k];
264  Cmidxcol[subsys[k]] = CmidxAcol[k];
265  }
266 
267  // we now got the total row/col indexes
268  idxrow = internal::_multiidx2n(Cmidxrow, n, Cdims);
269  idxcol = internal::_multiidx2n(Cmidxcol, n, Cdims);
270 
271  // check whether all CTRL row and col multi indexes are equal
272  bool all_ctrl_rows_equal = true;
273  bool all_ctrl_cols_equal = true;
274 
275  idx first_ctrl_row, first_ctrl_col;
276  if (ctrlsize > 0)
277  {
278  first_ctrl_row = CmidxCTRLrow[0];
279  first_ctrl_col = CmidxCTRLcol[0];
280  }
281  else
282  {
283  first_ctrl_row = first_ctrl_col = 1;
284  }
285 
286  for (idx k = 1; k < ctrlsize; ++k)
287  {
288  if (CmidxCTRLrow[k] != first_ctrl_row)
289  {
290  all_ctrl_rows_equal = false;
291  break;
292  }
293  }
294  for (idx k = 1; k < ctrlsize; ++k)
295  {
296  if (CmidxCTRLcol[k] != first_ctrl_col)
297  {
298  all_ctrl_cols_equal = false;
299  break;
300  }
301  }
302 
303  // at least one control activated, compute the coefficient
304  for (idx _n1 = 0; _n1 < DA; ++_n1)
305  {
306  internal::_n2multiidx(_n1, subsyssize, CdimsA, CmidxArow);
307  for (idx k = 0; k < subsyssize; ++k)
308  {
309  Cmidxrow[subsys[k]] = CmidxArow[k];
310  }
311  idx idxrowtmp = internal::_multiidx2n(Cmidxrow, n, Cdims);
312 
313  if (all_ctrl_rows_equal)
314  {
315  lhs = Ai[first_ctrl_row](_m1, _n1);
316  }
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  }
334  else
335  {
336  rhs = (_n2 == _m2) ? 1 : 0; // identity matrix
337  }
338 
339  idx idxcoltmp = internal::_multiidx2n(Cmidxcol, n, Cdims);
340 
341  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
342  }
343  }
344 
345  return std::make_tuple(coeff, idxrow, idxcol);
346  }; /* end coeff_idx_rho */
347 
348  //************ ket ************//
349  if (internal::_check_cvector(rstate)) // we have a ket
350  {
351  // check that dims match state vector
352  if (!internal::_check_dims_match_cvect(dims, rstate))
353  throw Exception("qpp::applyCTRL()",
355  if (D == 1)
356  return rstate;
357 
358  dyn_mat<typename Derived1::Scalar> result = rstate;
359 
360 #pragma omp parallel for collapse(2)
361  for (idx m = 0; m < DA; ++m)
362  for (idx r = 0; r < DCTRLAbar; ++r)
363  {
364  if (ctrlsize == 0) // no control
365  {
366  result(coeff_idx_ket(1, m, r).second) =
367  coeff_idx_ket(1, m, r).first;
368  }
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 #pragma omp parallel for collapse(4)
393  for (idx m1 = 0; m1 < DA; ++m1)
394  for (idx r1 = 0; r1 < DCTRLAbar; ++r1)
395  for (idx m2 = 0; m2 < DA; ++m2)
396  for (idx r2 = 0; r2 < DCTRLAbar; ++r2)
397  if (ctrlsize == 0) // no control
398  {
399  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
400  1, m2, r2);
401  result(std::get<1>(coeff_idxes),
402  std::get<2>(coeff_idxes)) =
403  std::get<0>(coeff_idxes);
404  }
405  else
406  {
407  for (idx i1 = 0; i1 < Dctrl; ++i1)
408  for (idx i2 = 0; i2 < Dctrl; ++i2)
409  {
410  auto coeff_idxes = coeff_idx_rho(
411  i1, m1, r1,
412  i2, m2, r2);
413  result(std::get<1>(coeff_idxes),
414  std::get<2>(coeff_idxes)) =
415  std::get<0>(coeff_idxes);
416  }
417  }
418 
419  return result;
420  }
421  //************ Exception: not ket nor density matrix ************//
422  else
423  throw Exception("qpp::applyCTRL()",
425 }
426 
442 template<typename Derived1, typename Derived2>
444  const Eigen::MatrixBase<Derived1>& state,
445  const Eigen::MatrixBase<Derived2>& A,
446  const std::vector<idx>& ctrl,
447  const std::vector<idx>& subsys,
448  idx d = 2)
449 {
450  const dyn_mat<typename Derived1::Scalar>& rstate = state;
452 
453  // EXCEPTION CHECKS
454 
455  // check zero size
456  if (!internal::_check_nonzero_size(rstate))
457  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
458 
459  // check valid dims
460  if (d == 0)
461  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
462  // END EXCEPTION CHECKS
463 
464  idx n =
465  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
466  std::log2(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 dyn_mat<typename Derived1::Scalar>& rstate = state;
494 
495  // EXCEPTION CHECKS
496 
497  // check types
498  if (!std::is_same<typename Derived1::Scalar,
499  typename Derived2::Scalar>::value)
500  throw Exception("qpp::apply()", Exception::Type::TYPE_MISMATCH);
501 
502  // check zero sizes
504  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
505 
506  // check zero sizes
507  if (!internal::_check_nonzero_size(rstate))
508  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
509 
510  // check square matrix for the gate
512  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
513 
514  // check that dimension is valid
515  if (!internal::_check_dims(dims))
516  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
517 
518  // check subsys is valid w.r.t. dims
519  if (!internal::_check_subsys_match_dims(subsys, dims))
520  throw Exception("qpp::apply()", Exception::Type::SUBSYS_MISMATCH_DIMS);
521 
522  // check that gate matches the dimensions of the subsys
523  std::vector<idx> subsys_dims(subsys.size());
524  for (idx i = 0; i < subsys.size(); ++i)
525  subsys_dims[i] = dims[subsys[i]];
526  if (!internal::_check_dims_match_mat(subsys_dims, rA))
527  throw Exception("qpp::apply()",
529  // END EXCEPTION CHECKS
530 
531  //************ ket ************//
532  if (internal::_check_cvector(rstate)) // we have a ket
533  {
534  // check that dims match state vector
535  if (!internal::_check_dims_match_cvect(dims, rstate))
536  throw Exception("qpp::apply()",
538 
539  return applyCTRL(rstate, rA, {}, subsys, dims);
540  }
541  //************ density matrix ************//
542  else if (internal::_check_square_mat(rstate)) // we have a density operator
543  {
544 
545  // check that dims match state matrix
546  if (!internal::_check_dims_match_mat(dims, rstate))
547  throw Exception("qpp::apply()",
549 
550  return applyCTRL(rstate, rA, {}, subsys, dims);
551  }
552  //************ Exception: not ket nor density matrix ************//
553  else
554  throw Exception("qpp::apply()",
556 }
557 
571 template<typename Derived1, typename Derived2>
573  const Eigen::MatrixBase<Derived1>& state,
574  const Eigen::MatrixBase<Derived2>& A,
575  const std::vector<idx>& subsys,
576  idx d = 2)
577 {
578  const dyn_mat<typename Derived1::Scalar>& rstate = state;
580 
581  // EXCEPTION CHECKS
582 
583  // check zero size
584  if (!internal::_check_nonzero_size(rstate))
585  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
586 
587  // check valid dims
588  if (d == 0)
589  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
590  // END EXCEPTION CHECKS
591 
592  idx n =
593  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
594  std::log2(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>& rho,
610  const std::vector<cmat>& Ks)
611 {
612  const cmat& rrho = rho;
613 
614  // EXCEPTION CHECKS
615 
617  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
618  if (!internal::_check_square_mat(rrho))
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() != rrho.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(rrho.rows(), rrho.rows());
633 
634 #pragma omp parallel for
635  for (idx i = 0; i < Ks.size(); ++i)
636  {
637 #pragma omp critical
638  {
639  result += Ks[i] * rrho * adjoint(Ks[i]);
640  }
641  }
642 
643  return result;
644 }
645 
656 template<typename Derived>
657 cmat apply(const Eigen::MatrixBase<Derived>& rho,
658  const std::vector<cmat>& Ks,
659  const std::vector<idx>& subsys,
660  const std::vector<idx>& dims)
661 {
662  const cmat& rrho = rho;
663 
664  // EXCEPTION CHECKS
665 
666  // check zero sizes
668  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
669 
670  // check square matrix for the rho
671  if (!internal::_check_square_mat(rrho))
672  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
673 
674  // check that dimension is valid
675  if (!internal::_check_dims(dims))
676  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
677 
678  // check that dims match rho matrix
679  if (!internal::_check_dims_match_mat(dims, rrho))
680  throw Exception("qpp::apply()",
682 
683  // check subsys is valid w.r.t. dims
684  if (!internal::_check_subsys_match_dims(subsys, dims))
685  throw Exception("qpp::apply()",
687 
688  std::vector<idx> subsys_dims(subsys.size());
689  for (idx i = 0; i < subsys.size(); ++i)
690  subsys_dims[i] = dims[subsys[i]];
691 
692  // check the Kraus operators
693  if (Ks.size() == 0)
694  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
695  if (!internal::_check_square_mat(Ks[0]))
696  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
697  if (!internal::_check_dims_match_mat(subsys_dims, Ks[0]))
698  throw Exception("qpp::apply()",
700  for (auto&& it : Ks)
701  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
702  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
703  // END EXCEPTION CHECKS
704 
705  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
706 
707  for (idx i = 0; i < Ks.size(); ++i)
708  result += apply(rrho, Ks[i], subsys, dims);
709 
710  return result;
711 }
712 
723 template<typename Derived>
724 cmat apply(const Eigen::MatrixBase<Derived>& rho,
725  const std::vector<cmat>& Ks,
726  const std::vector<idx>& subsys,
727  idx d = 2)
728 {
729  const cmat& rrho = rho;
730 
731  // EXCEPTION CHECKS
732 
733  // check zero sizes
735  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
736 
737  // check valid dims
738  if (d == 0)
739  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
740  // END EXCEPTION CHECKS
741 
742  idx n =
743  static_cast<idx>(std::llround(std::log2(rrho.rows()) /
744  std::log2(d)));
745  std::vector<idx> dims(n, d); // local dimensions vector
746 
747  return apply(rrho, Ks, subsys, dims);
748 }
749 
761 inline cmat kraus2super(const std::vector<cmat>& Ks)
762 {
763  // EXCEPTION CHECKS
764 
765  if (Ks.size() == 0)
766  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
767  if (!internal::_check_nonzero_size(Ks[0]))
768  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
769  if (!internal::_check_square_mat(Ks[0]))
770  throw Exception("qpp::kraus2super()",
772  for (auto&& it : Ks)
773  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
774  throw Exception("qpp::kraus2super()",
776  // END EXCEPTION CHECKS
777 
778  idx D = static_cast<idx>(Ks[0].rows());
779 
780  cmat result(D * D, D * D);
781  cmat MN = cmat::Zero(D, D);
782  bra A = bra::Zero(D);
783  ket B = ket::Zero(D);
784  cmat EMN = cmat::Zero(D, D);
785 
786 #pragma omp parallel for collapse(2)
787  for (idx m = 0; m < D; ++m)
788  {
789  for (idx n = 0; n < D; ++n)
790  {
791 #pragma omp critical
792  {
793  // compute E(|m><n|)
794  MN(m, n) = 1;
795  for (idx i = 0; i < Ks.size(); ++i)
796  EMN += Ks[i] * MN * adjoint(Ks[i]);
797  MN(m, n) = 0;
798 
799  for (idx a = 0; a < D; ++a)
800  {
801  A(a) = 1;
802  for (idx b = 0; b < D; ++b)
803  {
804  // compute result(ab,mn)=<a|E(|m><n)|b>
805  B(b) = 1;
806  result(a * D + b, m * D + n) =
807  static_cast<cmat>(A * EMN * B).value();
808  B(b) = 0;
809  }
810  A(a) = 0;
811  }
812  EMN = cmat::Zero(D, D);
813  }
814  }
815  }
816 
817  return result;
818 }
819 
835 inline cmat kraus2choi(const std::vector<cmat>& Ks)
836 {
837  // EXCEPTION CHECKS
838 
839  if (Ks.size() == 0)
840  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
841  if (!internal::_check_nonzero_size(Ks[0]))
842  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
843  if (!internal::_check_square_mat(Ks[0]))
844  throw Exception("qpp::kraus2choi()",
846  for (auto&& it : Ks)
847  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
848  throw Exception("qpp::kraus2choi()",
850  // END EXCEPTION CHECKS
851 
852  idx D = static_cast<idx>(Ks[0].rows());
853 
854  // construct the D x D \sum |jj> vector
855  // (un-normalized maximally entangled state)
856  cmat MES = cmat::Zero(D * D, 1);
857  for (idx a = 0; a < D; ++a)
858  MES(a * D + a) = 1;
859 
860  cmat Omega = MES * adjoint(MES);
861 
862  cmat result = cmat::Zero(D * D, D * D);
863 
864 #pragma omp parallel for
865  for (idx i = 0; i < Ks.size(); ++i)
866  {
867 #pragma omp critical
868  {
869  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
870  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
871  }
872  }
873 
874  return result;
875 }
876 
890 inline std::vector<cmat> choi2kraus(const cmat& A)
891 {
892  // EXCEPTION CHECKS
893 
895  throw Exception("qpp::choi2kraus()", Exception::Type::ZERO_SIZE);
897  throw Exception("qpp::choi2kraus()",
899  idx D = static_cast<idx>(std::llround(
900  std::sqrt(static_cast<double>(A.rows()))));
901  if (D * D != static_cast<idx>(A.rows()))
902  throw Exception("qpp::choi2kraus()", Exception::Type::DIMS_INVALID);
903  // END EXCEPTION CHECKS
904 
905  dmat ev = hevals(A);
906  cmat evec = hevects(A);
907  std::vector<cmat> result;
908 
909  for (idx i = 0; i < D * D; ++i)
910  {
911  if (std::abs(ev(i)) > eps)
912  result.push_back(
913  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
914  }
915 
916  return result;
917 }
918 
926 inline cmat choi2super(const cmat& A)
927 {
928  // EXCEPTION CHECKS
929 
931  throw Exception("qpp::choi2super()", Exception::Type::ZERO_SIZE);
933  throw Exception("qpp::choi2super()",
935  idx D = static_cast<idx>(std::llround(
936  std::sqrt(static_cast<double>(A.rows()))));
937  if (D * D != static_cast<idx>(A.rows()))
938  throw Exception("qpp::choi2super()", Exception::Type::DIMS_INVALID);
939  // END EXCEPTION CHECKS
940 
941  cmat result(D * D, D * D);
942 
943 #pragma omp parallel for collapse(4)
944  for (idx a = 0; a < D; ++a)
945  for (idx b = 0; b < D; ++b)
946  for (idx m = 0; m < D; ++m)
947  for (idx n = 0; n < D; ++n)
948  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
949 
950  return result;
951 }
952 
960 inline cmat super2choi(const cmat& A)
961 {
962  // EXCEPTION CHECKS
963 
965  throw Exception("qpp::super2choi()", Exception::Type::ZERO_SIZE);
967  throw Exception("qpp::super2choi()",
969  idx D = static_cast<idx>(std::llround(
970  std::sqrt(static_cast<double>(A.rows()))));
971  if (D * D != static_cast<idx>(A.rows()))
972  throw Exception("qpp::super2choi()", Exception::Type::DIMS_INVALID);
973  // END EXCEPTION CHECKS
974 
975  cmat result(D * D, D * D);
976 
977 #pragma omp parallel for collapse(4)
978  for (idx a = 0; a < D; ++a)
979  for (idx b = 0; b < D; ++b)
980  for (idx m = 0; m < D; ++m)
981  for (idx n = 0; n < D; ++n)
982  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
983 
984  return result;
985 }
986 
1001 template<typename Derived>
1002 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1003  const std::vector<idx>& dims)
1004 {
1005  const dyn_mat<typename Derived::Scalar>& rA = A;
1006 
1007  // EXCEPTION CHECKS
1008 
1009  // check zero-size
1011  throw Exception("qpp::ptrace1()", Exception::Type::ZERO_SIZE);
1012 
1013  // check that dims is a valid dimension vector
1014  if (!internal::_check_dims(dims))
1015  throw Exception("qpp::ptrace1()", Exception::Type::DIMS_INVALID);
1016 
1017  // check dims has only 2 elements
1018  if (dims.size() != 2)
1019  throw Exception("qpp::ptrace1()", Exception::Type::NOT_BIPARTITE);
1020  // END EXCEPTION CHECKS
1021 
1022  idx DA = dims[0];
1023  idx DB = dims[1];
1024 
1027 
1028  //************ ket ************//
1029  if (internal::_check_cvector(rA)) // we have a ket
1030  {
1031  // check that dims match the dimension of A
1032  if (!internal::_check_dims_match_cvect(dims, rA))
1033  throw Exception("qpp::ptrace1()",
1035 
1036  auto worker = [=](idx i, idx j) noexcept
1037  -> typename Derived::Scalar
1038  {
1039  typename Derived::Scalar sum = 0;
1040  for (idx m = 0; m < DA; ++m)
1041  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
1042 
1043  return sum;
1044  }; /* end worker */
1045 
1046 #pragma omp parallel for collapse(2)
1047  for (idx j = 0; j < DB; ++j) // column major order for speed
1048  for (idx i = 0; i < DB; ++i)
1049  result(i, j) = worker(i, j);
1050 
1051  return result;
1052  }
1053  //************ density matrix ************//
1054  else if (internal::_check_square_mat(rA)) // we have a density operator
1055  {
1056  // check that dims match the dimension of A
1057  if (!internal::_check_dims_match_mat(dims, rA))
1058  throw Exception("qpp::ptrace1()",
1060 
1061  auto worker = [=](idx i, idx j) noexcept
1062  -> typename Derived::Scalar
1063  {
1064  typename Derived::Scalar sum = 0;
1065  for (idx m = 0; m < DA; ++m)
1066  sum += rA(m * DB + i, m * DB + j);
1067 
1068  return sum;
1069  }; /* end worker */
1070 
1071 #pragma omp parallel for collapse(2)
1072  for (idx j = 0; j < DB; ++j) // column major order for speed
1073  for (idx i = 0; i < DB; ++i)
1074  result(i, j) = worker(i, j);
1075 
1076  return result;
1077  }
1078  //************ Exception: not ket nor density matrix ************//
1079  else
1080  throw Exception("qpp::ptrace1()",
1082 }
1083 
1098 template<typename Derived>
1099 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1100  const std::vector<idx>& dims)
1101 {
1102  const dyn_mat<typename Derived::Scalar>& rA = A;
1103 
1104  // EXCEPTION CHECKS
1105 
1106  // check zero-size
1108  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1109 
1110  // check that dims is a valid dimension vector
1111  if (!internal::_check_dims(dims))
1112  throw Exception("qpp::ptrace2()", Exception::Type::DIMS_INVALID);
1113 
1114  // check dims has only 2 elements
1115  if (dims.size() != 2)
1116  throw Exception("qpp::ptrace2()", Exception::Type::NOT_BIPARTITE);
1117  // END EXCEPTION CHECKS
1118 
1119  idx DA = dims[0];
1120  idx DB = dims[1];
1121 
1124 
1125  //************ ket ************//
1126  if (internal::_check_cvector(rA)) // we have a ket
1127  {
1128  // check that dims match the dimension of A
1129  if (!internal::_check_dims_match_cvect(dims, rA))
1130  throw Exception("qpp::ptrace2()",
1132 
1133  auto worker = [=](idx i, idx j) noexcept
1134  -> typename Derived::Scalar
1135  {
1136  typename Derived::Scalar sum = 0;
1137  for (idx m = 0; m < DB; ++m)
1138  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1139 
1140  return sum;
1141  }; /* end worker */
1142 
1143 #pragma omp parallel for collapse(2)
1144  for (idx j = 0; j < DA; ++j) // column major order for speed
1145  for (idx i = 0; i < DA; ++i)
1146  result(i, j) = worker(i, j);
1147 
1148  return result;
1149  }
1150  //************ density matrix ************//
1151  else if (internal::_check_square_mat(rA)) // we have a density operator
1152  {
1153  // check that dims match the dimension of A
1154  if (!internal::_check_dims_match_mat(dims, rA))
1155  throw Exception("qpp::ptrace2()",
1157 
1158 #pragma omp parallel for collapse(2)
1159  for (idx j = 0; j < DA; ++j) // column major order for speed
1160  for (idx i = 0; i < DA; ++i)
1161  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1162 
1163  return result;
1164  }
1165  //************ Exception: not ket nor density matrix ************//
1166  else
1167  throw Exception("qpp::ptrace1()",
1169 }
1170 
1185 template<typename Derived>
1186 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1187  const std::vector<idx>& subsys,
1188  const std::vector<idx>& dims)
1189 {
1190  const dyn_mat<typename Derived::Scalar>& rA = A;
1191 
1192  // EXCEPTION CHECKS
1193 
1194  // check zero-size
1196  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1197 
1198  // check that dims is a valid dimension vector
1199  if (!internal::_check_dims(dims))
1200  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1201 
1202  // check that subsys are valid
1203  if (!internal::_check_subsys_match_dims(subsys, dims))
1204  throw Exception("qpp::ptrace()",
1206  // END EXCEPTION CHECKS
1207 
1208  idx D = static_cast<idx>(rA.rows());
1209  idx n = dims.size();
1210  idx nsubsys = subsys.size();
1211  idx nsubsysbar = n - nsubsys;
1212  idx Dsubsys = 1;
1213  for (idx i = 0; i < nsubsys; ++i)
1214  Dsubsys *= dims[subsys[i]];
1215  idx Dsubsysbar = D / Dsubsys;
1216 
1217  idx Cdims[maxn];
1218  idx Csubsys[maxn];
1219  idx Cdimssubsys[maxn];
1220  idx Csubsysbar[maxn];
1221  idx Cdimssubsysbar[maxn];
1222 
1223  idx Cmidxcolsubsysbar[maxn];
1224 
1225  std::vector<idx> subsys_bar = complement(subsys, n);
1226  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1227  std::begin(Csubsysbar));
1228 
1229  for (idx i = 0; i < n; ++i)
1230  {
1231  Cdims[i] = dims[i];
1232  }
1233  for (idx i = 0; i < nsubsys; ++i)
1234  {
1235  Csubsys[i] = subsys[i];
1236  Cdimssubsys[i] = dims[subsys[i]];
1237  }
1238  for (idx i = 0; i < nsubsysbar; ++i)
1239  {
1240  Cdimssubsysbar[i] = dims[subsys_bar[i]];
1241  }
1242 
1244  dyn_mat<typename Derived::Scalar>(Dsubsysbar, Dsubsysbar);
1245 
1246  //************ ket ************//
1247  if (internal::_check_cvector(rA)) // we have a ket
1248  {
1249  // check that dims match the dimension of A
1250  if (!internal::_check_dims_match_cvect(dims, rA))
1251  throw Exception("qpp::ptrace()",
1253 
1254  if (subsys.size() == dims.size())
1255  {
1256  result(0, 0) = (adjoint(rA) * rA).value();
1257  return result;
1258  }
1259 
1260  if (subsys.size() == 0)
1261  return rA * adjoint(rA);
1262 
1263  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1264  -> typename Derived::Scalar
1265  {
1266  // use static allocation for speed!
1267 
1268  idx Cmidxrow[maxn];
1269  idx Cmidxcol[maxn];
1270  idx Cmidxrowsubsysbar[maxn];
1271  idx Cmidxsubsys[maxn];
1272 
1273  /* get the row multi-indexes of the complement */
1274  internal::_n2multiidx(i, nsubsysbar,
1275  Cdimssubsysbar, Cmidxrowsubsysbar);
1276  /* write them in the global row/col multi-indexes */
1277  for (idx k = 0; k < nsubsysbar; ++k)
1278  {
1279  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1280  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1281  }
1282  typename Derived::Scalar sm = 0;
1283  for (idx a = 0; a < Dsubsys; ++a)
1284  {
1285  // get the multi-index over which we do the summation
1286  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1287  // write it into the global row/col multi-indexes
1288  for (idx k = 0; k < nsubsys; ++k)
1289  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1290  = Cmidxsubsys[k];
1291 
1292  // now do the sum
1293  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims)) *
1294  std::conj(rA(internal::_multiidx2n(Cmidxcol, n,
1295  Cdims)));
1296  }
1297 
1298  return sm;
1299  }; /* end worker */
1300 
1301  for (idx j = 0; j < Dsubsysbar; ++j) // column major order for speed
1302  {
1303  // compute the column multi-indexes of the complement
1304  internal::_n2multiidx(j, nsubsysbar,
1305  Cdimssubsysbar, Cmidxcolsubsysbar);
1306 #pragma omp parallel for
1307  for (idx i = 0; i < Dsubsysbar; ++i)
1308  {
1309  result(i, j) = worker(i);
1310  }
1311  }
1312 
1313  return result;
1314  }
1315  //************ density matrix ************//
1316  else if (internal::_check_square_mat(rA)) // we have a density operator
1317  {
1318  // check that dims match the dimension of A
1319  if (!internal::_check_dims_match_mat(dims, rA))
1320  throw Exception("qpp::ptrace()",
1322 
1323  if (subsys.size() == dims.size())
1324  {
1325  result(0, 0) = rA.trace();
1326  return result;
1327  }
1328 
1329  if (subsys.size() == 0)
1330  return rA;
1331 
1332  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1333  -> typename Derived::Scalar
1334  {
1335  // use static allocation for speed!
1336 
1337  idx Cmidxrow[maxn];
1338  idx Cmidxcol[maxn];
1339  idx Cmidxrowsubsysbar[maxn];
1340  idx Cmidxsubsys[maxn];
1341 
1342  /* get the row/col multi-indexes of the complement */
1343  internal::_n2multiidx(i, nsubsysbar,
1344  Cdimssubsysbar, Cmidxrowsubsysbar);
1345  /* write them in the global row/col multi-indexes */
1346  for (idx k = 0; k < nsubsysbar; ++k)
1347  {
1348  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1349  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1350  }
1351  typename Derived::Scalar sm = 0;
1352  for (idx a = 0; a < Dsubsys; ++a)
1353  {
1354  // get the multi-index over which we do the summation
1355  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1356  // write it into the global row/col multi-indexes
1357  for (idx k = 0; k < nsubsys; ++k)
1358  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1359  = Cmidxsubsys[k];
1360 
1361  // now do the sum
1362  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims),
1363  internal::_multiidx2n(Cmidxcol, n, Cdims));
1364  }
1365 
1366  return sm;
1367  }; /* end worker */
1368 
1369  for (idx j = 0; j < Dsubsysbar; ++j) // column major order for speed
1370  {
1371  // compute the column multi-indexes of the complement
1372  internal::_n2multiidx(j, nsubsysbar,
1373  Cdimssubsysbar, Cmidxcolsubsysbar);
1374 #pragma omp parallel for
1375  for (idx i = 0; i < Dsubsysbar; ++i)
1376  {
1377  result(i, j) = worker(i);
1378  }
1379  }
1380 
1381  return result;
1382  }
1383  //************ Exception: not ket nor density matrix ************//
1384  else
1385  throw Exception("qpp::ptrace()",
1387 }
1388 
1403 template<typename Derived>
1404 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1405  const std::vector<idx>& subsys,
1406  idx d = 2)
1407 {
1408  const dyn_mat<typename Derived::Scalar>& rA = A;
1409 
1410  // EXCEPTION CHECKS
1411 
1412  // check zero size
1414  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1415 
1416  // check valid dims
1417  if (d == 0)
1418  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1419  // END EXCEPTION CHECKS
1420 
1421  idx n =
1422  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1423  std::log2(d)));
1424  std::vector<idx> dims(n, d); // local dimensions vector
1425 
1426  return ptrace(rA, subsys, dims);
1427 }
1428 
1442 template<typename Derived>
1444  const Eigen::MatrixBase<Derived>& A,
1445  const std::vector<idx>& subsys,
1446  const std::vector<idx>& dims)
1447 {
1448  const dyn_mat<typename Derived::Scalar>& rA = A;
1449 
1450  // EXCEPTION CHECKS
1451 
1452  // check zero-size
1454  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1455 
1456  // check that dims is a valid dimension vector
1457  if (!internal::_check_dims(dims))
1458  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1459 
1460  // check that subsys are valid
1461  if (!internal::_check_subsys_match_dims(subsys, dims))
1462  throw Exception("qpp::ptranspose()",
1464  // END EXCEPTION CHECKS
1465 
1466  idx D = static_cast<idx>(rA.rows());
1467  idx numdims = dims.size();
1468  idx numsubsys = subsys.size();
1469  idx Cdims[maxn];
1470  idx Cmidxcol[maxn];
1471  idx Csubsys[maxn];
1472 
1473  // copy dims in Cdims and subsys in Csubsys
1474  for (idx i = 0; i < numdims; ++i)
1475  Cdims[i] = dims[i];
1476  for (idx i = 0; i < numsubsys; ++i)
1477  Csubsys[i] = subsys[i];
1478 
1479  dyn_mat<typename Derived::Scalar> result(D, D);
1480 
1481  //************ ket ************//
1482  if (internal::_check_cvector(rA)) // we have a ket
1483  {
1484  // check that dims match the dimension of A
1485  if (!internal::_check_dims_match_cvect(dims, rA))
1486  throw Exception("qpp::ptranspose()",
1488 
1489  if (subsys.size() == dims.size())
1490  return (rA * adjoint(rA)).transpose();
1491 
1492  if (subsys.size() == 0)
1493  return rA * adjoint(rA);
1494 
1495  auto worker = [=, &Cmidxcol](idx i) noexcept
1496  -> typename Derived::Scalar
1497  {
1498  // use static allocation for speed!
1499  idx midxcoltmp[maxn];
1500  idx midxrow[maxn];
1501 
1502  for (idx k = 0; k < numdims; ++k)
1503  midxcoltmp[k] = Cmidxcol[k];
1504 
1505  /* compute the row multi-index */
1506  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1507 
1508  for (idx k = 0; k < numsubsys; ++k)
1509  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1510 
1511  /* writes the result */
1512  return rA(internal::_multiidx2n(midxrow, numdims, Cdims)) *
1513  std::conj(rA(internal::_multiidx2n(midxcoltmp, numdims,
1514  Cdims)));
1515  }; /* end worker */
1516 
1517  for (idx j = 0; j < D; ++j)
1518  {
1519  // compute the column multi-index
1520  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1521 #pragma omp parallel for
1522  for (idx i = 0; i < D; ++i)
1523  result(i, j) = worker(i);
1524  }
1525 
1526  return result;
1527  }
1528  //************ density matrix ************//
1529  else if (internal::_check_square_mat(rA)) // we have a density operator
1530  {
1531  // check that dims match the dimension of A
1532  if (!internal::_check_dims_match_mat(dims, rA))
1533  throw Exception("qpp::ptranspose()",
1535 
1536  if (subsys.size() == dims.size())
1537  return rA.transpose();
1538 
1539  if (subsys.size() == 0)
1540  return rA;
1541 
1542  auto worker = [=, &Cmidxcol](idx i) noexcept
1543  -> typename Derived::Scalar
1544  {
1545  // use static allocation for speed!
1546  idx midxcoltmp[maxn];
1547  idx midxrow[maxn];
1548 
1549  for (idx k = 0; k < numdims; ++k)
1550  midxcoltmp[k] = Cmidxcol[k];
1551 
1552  /* compute the row multi-index */
1553  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1554 
1555  for (idx k = 0; k < numsubsys; ++k)
1556  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1557 
1558  /* writes the result */
1559  return rA(internal::_multiidx2n(midxrow, numdims, Cdims),
1560  internal::_multiidx2n(midxcoltmp, numdims, Cdims));
1561  }; /* end worker */
1562 
1563  for (idx j = 0; j < D; ++j)
1564  {
1565  // compute the column multi-index
1566  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1567 #pragma omp parallel for
1568  for (idx i = 0; i < D; ++i)
1569  result(i, j) = worker(i);
1570  }
1571 
1572  return result;
1573  }
1574  //************ Exception: not ket nor density matrix ************//
1575  else
1576  throw Exception("qpp::ptranspose()",
1578 }
1579 
1593 template<typename Derived>
1595  const Eigen::MatrixBase<Derived>& A,
1596  const std::vector<idx>& subsys,
1597  idx d = 2)
1598 {
1599  const dyn_mat<typename Derived::Scalar>& rA = A;
1600 
1601  // EXCEPTION CHECKS
1602 
1603  // check zero size
1605  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1606 
1607  // check valid dims
1608  if (d == 0)
1609  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1610  // END EXCEPTION CHECKS
1611 
1612  idx n =
1613  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1614  std::log2(d)));
1615  std::vector<idx> dims(n, d); // local dimensions vector
1616 
1617  return ptranspose(rA, subsys, dims);
1618 }
1619 
1632 template<typename Derived>
1634  const Eigen::MatrixBase<Derived>& A,
1635  const std::vector<idx>& perm,
1636  const std::vector<idx>& dims)
1637 {
1638  const dyn_mat<typename Derived::Scalar>& rA = A;
1639 
1640  // EXCEPTION CHECKS
1641 
1642  // check zero-size
1644  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1645 
1646  // check that dims is a valid dimension vector
1647  if (!internal::_check_dims(dims))
1648  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1649 
1650  // check that we have a valid permutation
1651  if (!internal::_check_perm(perm))
1652  throw Exception("qpp::syspermute()", Exception::Type::PERM_INVALID);
1653 
1654  // check that permutation match dimensions
1655  if (perm.size() != dims.size())
1656  throw Exception("qpp::syspermute()",
1658  // END EXCEPTION CHECKS
1659 
1660  idx D = static_cast<idx>(rA.rows());
1661  idx numdims = dims.size();
1662 
1664 
1665  //************ ket ************//
1666  if (internal::_check_cvector(rA)) // we have a column vector
1667  {
1668  idx Cdims[maxn];
1669  idx Cperm[maxn];
1670 
1671  // check that dims match the dimension of rA
1672  if (!internal::_check_dims_match_cvect(dims, rA))
1673  throw Exception("qpp::syspermute()",
1675 
1676  // copy dims in Cdims and perm in Cperm
1677  for (idx i = 0; i < numdims; ++i)
1678  {
1679  Cdims[i] = dims[i];
1680  Cperm[i] = perm[i];
1681  }
1682  result.resize(D, 1);
1683 
1684  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1685  -> idx
1686  {
1687  // use static allocation for speed,
1688  // double the size for matrices reshaped as vectors
1689  idx midx[maxn];
1690  idx midxtmp[maxn];
1691  idx permdims[maxn];
1692 
1693  /* compute the multi-index */
1694  internal::_n2multiidx(i, numdims, Cdims, midx);
1695 
1696  for (idx k = 0; k < numdims; ++k)
1697  {
1698  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1699  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1700  }
1701 
1702  return internal::_multiidx2n(midxtmp, numdims, permdims);
1703  }; /* end worker */
1704 
1705 #pragma omp parallel for
1706  for (idx i = 0; i < D; ++i)
1707  result(worker(i)) = rA(i);
1708 
1709  return result;
1710  }
1711  //************ density matrix ************//
1712  else if (internal::_check_square_mat(rA)) // we have a density operator
1713  {
1714  idx Cdims[2 * maxn];
1715  idx Cperm[2 * maxn];
1716 
1717  // check that dims match the dimension of rA
1718  if (!internal::_check_dims_match_mat(dims, rA))
1719  throw Exception("qpp::syspermute()",
1721 
1722  // copy dims in Cdims and perm in Cperm
1723  for (idx i = 0; i < numdims; ++i)
1724  {
1725  Cdims[i] = dims[i];
1726  Cdims[i + numdims] = dims[i];
1727  Cperm[i] = perm[i];
1728  Cperm[i + numdims] = perm[i] + numdims;
1729  }
1730  result.resize(D * D, 1);
1731  // map A to a column vector
1733  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1734  const_cast<typename Derived::Scalar*>(rA.data()), D * D,
1735  1);
1736 
1737  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1738  -> idx
1739  {
1740  // use static allocation for speed,
1741  // double the size for matrices reshaped as vectors
1742  idx midx[2 * maxn];
1743  idx midxtmp[2 * maxn];
1744  idx permdims[2 * maxn];
1745 
1746  /* compute the multi-index */
1747  internal::_n2multiidx(i, 2 * numdims, Cdims, midx);
1748 
1749  for (idx k = 0; k < 2 * numdims; ++k)
1750  {
1751  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1752  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1753  }
1754 
1755  return internal::_multiidx2n(midxtmp, 2 * numdims, permdims);
1756  }; /* end worker */
1757 
1758 #pragma omp parallel for
1759  for (idx i = 0; i < D * D; ++i)
1760  result(worker(i)) = rA(i);
1761 
1762  return reshape(result, D, D);
1763  }
1764  //************ Exception: not ket nor density matrix ************//
1765  else
1766  throw Exception("qpp::syspermute()",
1768 }
1769 
1782 template<typename Derived>
1784  const Eigen::MatrixBase<Derived>& A,
1785  const std::vector<idx>& perm,
1786  idx d = 2)
1787 {
1788  const dyn_mat<typename Derived::Scalar>& rA = A;
1789 
1790  // EXCEPTION CHECKS
1791 
1792  // check zero size
1794  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1795 
1796  // check valid dims
1797  if (d == 0)
1798  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1799  // END EXCEPTION CHECKS
1800 
1801  idx n =
1802  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1803  std::log2(d)));
1804  std::vector<idx> dims(n, d); // local dimensions vector
1805 
1806  return syspermute(rA, perm, dims);
1807 }
1808 
1809 } /* namespace qpp */
1810 
1811 #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:890
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:1633
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:486
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:391
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:926
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:960
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1806
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1002
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:890
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1099
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:783
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:417
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:835
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:1186
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:761
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:1443
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:209
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1130