Quantum++  v1.0-rc3
A modern C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * This file is part of Quantum++.
3  *
4  * MIT License
5  *
6  * Copyright (c) 2013 - 2018 Vlad Gheorghiu (vgheorgh@gmail.com)
7  *
8  * Permission is hereby granted, free of charge, to any person obtaining a copy
9  * of this software and associated documentation files (the "Software"), to deal
10  * in the Software without restriction, including without limitation the rights
11  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12  * copies of the Software, and to permit persons to whom the Software is
13  * furnished to do so, subject to the following conditions:
14  *
15  * The above copyright notice and this permission notice shall be included in
16  * all copies or substantial portions of the Software.
17  *
18  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24  * SOFTWARE.
25  */
26 
32 #ifndef OPERATIONS_H_
33 #define OPERATIONS_H_
34 
35 namespace qpp {
36 
53 template <typename Derived1, typename Derived2>
54 dyn_mat<typename Derived1::Scalar>
55 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
56  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
57  const std::vector<idx>& subsys, const std::vector<idx>& dims) {
58  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
59  state.derived();
60  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
61 
62  // EXCEPTION CHECKS
63 
64  // check types
65  if (!std::is_same<typename Derived1::Scalar,
66  typename Derived2::Scalar>::value)
67  throw exception::TypeMismatch("qpp::applyCTRL()");
68 
69  // check zero sizes
71  throw exception::ZeroSize("qpp::applyCTRL()");
72 
73  // check zero sizes
74  if (!internal::check_nonzero_size(rstate))
75  throw exception::ZeroSize("qpp::applyCTRL()");
76 
77  // check square matrix for the gate
79  throw exception::MatrixNotSquare("qpp::applyCTRL()");
80 
81  // check that all control subsystems have the same dimension
82  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
83  for (idx i = 1; i < ctrl.size(); ++i)
84  if (dims[ctrl[i]] != d)
85  throw exception::DimsNotEqual("qpp::applyCTRL()");
86 
87  // check that dimension is valid
88  if (!internal::check_dims(dims))
89  throw exception::DimsInvalid("qpp::applyCTRL()");
90 
91  // check subsys is valid w.r.t. dims
92  if (!internal::check_subsys_match_dims(subsys, dims))
93  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
94 
95  // check that gate matches the dimensions of the subsys
96  std::vector<idx> subsys_dims(subsys.size());
97  for (idx i = 0; i < subsys.size(); ++i)
98  subsys_dims[i] = dims[subsys[i]];
99  if (!internal::check_dims_match_mat(subsys_dims, rA))
100  throw exception::MatrixMismatchSubsys("qpp::applyCTRL()");
101 
102  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
103  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
104  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
105 
106  // check that ctrl + gate subsystem is valid
107  // with respect to local dimensions
108  if (!internal::check_subsys_match_dims(ctrlgate, dims))
109  throw exception::SubsysMismatchDims("qpp::applyCTRL()");
110  // END EXCEPTION CHECKS
111 
112  // construct the table of A^i and (A^dagger)^i
113  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
114  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
115  for (idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i) {
116  Ai.push_back(powm(rA, i));
117  Aidagger.push_back(powm(adjoint(rA), i));
118  }
119 
120  idx D = static_cast<idx>(rstate.rows()); // total dimension
121  idx N = dims.size(); // total number of subsystems
122  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
123  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
124  idx subsyssize = subsys.size(); // number of subsystems of the target
125  // dimension of ctrl subsystem
126  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
127  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
128 
129  idx Cdims[maxn]; // local dimensions
130  idx CdimsA[maxn]; // local dimensions
131  idx CdimsCTRL[maxn]; // local dimensions
132  idx CdimsCTRLA_bar[maxn]; // local dimensions
133 
134  // compute the complementary subsystem of ctrlgate w.r.t. dims
135  std::vector<idx> ctrlgate_bar = complement(ctrlgate, N);
136  // number of subsystems that are complementary to the ctrl+gate
137  idx ctrlgate_barsize = ctrlgate_bar.size();
138 
139  idx DCTRLA_bar = 1; // dimension of the rest
140  for (idx i = 0; i < ctrlgate_barsize; ++i)
141  DCTRLA_bar *= dims[ctrlgate_bar[i]];
142 
143  for (idx k = 0; k < N; ++k)
144  Cdims[k] = dims[k];
145  for (idx k = 0; k < subsyssize; ++k)
146  CdimsA[k] = dims[subsys[k]];
147  for (idx k = 0; k < ctrlsize; ++k)
148  CdimsCTRL[k] = d;
149  for (idx k = 0; k < ctrlgate_barsize; ++k)
150  CdimsCTRLA_bar[k] = dims[ctrlgate_bar[k]];
151 
152  // worker, computes the coefficient and the index for the ket case
153  // used in #pragma omp parallel for collapse
154  auto coeff_idx_ket = [&](idx i_, idx m_, idx r_) noexcept
155  ->std::pair<typename Derived1::Scalar, idx> {
156  idx indx = 0;
157  typename Derived1::Scalar coeff = 0;
158 
159  idx Cmidx[maxn]; // the total multi-index
160  idx CmidxA[maxn]; // the gate part multi-index
161  idx CmidxCTRLA_bar[maxn]; // the rest multi-index
162 
163  // compute the index
164 
165  // set the CTRL part
166  for (idx k = 0; k < ctrlsize; ++k) {
167  Cmidx[ctrl[k]] = i_;
168  }
169 
170  // set the rest
171  internal::n2multiidx(r_, N - ctrlgatesize, CdimsCTRLA_bar,
172  CmidxCTRLA_bar);
173  for (idx k = 0; k < N - ctrlgatesize; ++k) {
174  Cmidx[ctrlgate_bar[k]] = CmidxCTRLA_bar[k];
175  }
176 
177  // set the A part
178  internal::n2multiidx(m_, subsyssize, CdimsA, CmidxA);
179  for (idx k = 0; k < subsyssize; ++k) {
180  Cmidx[subsys[k]] = CmidxA[k];
181  }
182 
183  // we now got the total index
184  indx = internal::multiidx2n(Cmidx, N, Cdims);
185 
186  // compute the coefficient
187  for (idx n_ = 0; n_ < DA; ++n_) {
188  internal::n2multiidx(n_, subsyssize, CdimsA, CmidxA);
189  for (idx k = 0; k < subsyssize; ++k) {
190  Cmidx[subsys[k]] = CmidxA[k];
191  }
192  coeff +=
193  Ai[i_](m_, n_) * rstate(internal::multiidx2n(Cmidx, N, Cdims));
194  }
195 
196  return std::make_pair(coeff, indx);
197  }; /* end coeff_idx_ket */
198 
199  // worker, computes the coefficient and the index
200  // for the density matrix case
201  // used in #pragma omp parallel for collapse
202  auto coeff_idx_rho = [&](idx i1_, idx m1_, idx r1_, idx i2_, idx m2_,
203  idx r2_) noexcept
204  ->std::tuple<typename Derived1::Scalar, idx, idx> {
205  idx idxrow = 0;
206  idx idxcol = 0;
207  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
208 
209  idx Cmidxrow[maxn]; // the total row multi-index
210  idx Cmidxcol[maxn]; // the total col multi-index
211  idx CmidxArow[maxn]; // the gate part row multi-index
212  idx CmidxAcol[maxn]; // the gate part col multi-index
213  idx CmidxCTRLrow[maxn]; // the control row multi-index
214  idx CmidxCTRLcol[maxn]; // the control col multi-index
215  idx CmidxCTRLA_barrow[maxn]; // the rest row multi-index
216  idx CmidxCTRLA_barcol[maxn]; // the rest col multi-index
217 
218  // compute the ket/bra indexes
219 
220  // set the CTRL part
221  internal::n2multiidx(i1_, ctrlsize, CdimsCTRL, CmidxCTRLrow);
222  internal::n2multiidx(i2_, ctrlsize, CdimsCTRL, CmidxCTRLcol);
223 
224  for (idx k = 0; k < ctrlsize; ++k) {
225  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
226  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
227  }
228 
229  // set the rest
230  internal::n2multiidx(r1_, N - ctrlgatesize, CdimsCTRLA_bar,
231  CmidxCTRLA_barrow);
232  internal::n2multiidx(r2_, N - ctrlgatesize, CdimsCTRLA_bar,
233  CmidxCTRLA_barcol);
234  for (idx k = 0; k < N - ctrlgatesize; ++k) {
235  Cmidxrow[ctrlgate_bar[k]] = CmidxCTRLA_barrow[k];
236  Cmidxcol[ctrlgate_bar[k]] = CmidxCTRLA_barcol[k];
237  }
238 
239  // set the A part
240  internal::n2multiidx(m1_, subsyssize, CdimsA, CmidxArow);
241  internal::n2multiidx(m2_, subsyssize, CdimsA, CmidxAcol);
242  for (idx k = 0; k < subsys.size(); ++k) {
243  Cmidxrow[subsys[k]] = CmidxArow[k];
244  Cmidxcol[subsys[k]] = CmidxAcol[k];
245  }
246 
247  // we now got the total row/col indexes
248  idxrow = internal::multiidx2n(Cmidxrow, N, Cdims);
249  idxcol = internal::multiidx2n(Cmidxcol, N, Cdims);
250 
251  // check whether all CTRL row and col multi indexes are equal
252  bool all_ctrl_rows_equal = true;
253  bool all_ctrl_cols_equal = true;
254 
255  idx first_ctrl_row, first_ctrl_col;
256  if (ctrlsize > 0) {
257  first_ctrl_row = CmidxCTRLrow[0];
258  first_ctrl_col = CmidxCTRLcol[0];
259  } else {
260  first_ctrl_row = first_ctrl_col = 1;
261  }
262 
263  for (idx k = 1; k < ctrlsize; ++k) {
264  if (CmidxCTRLrow[k] != first_ctrl_row) {
265  all_ctrl_rows_equal = false;
266  break;
267  }
268  }
269  for (idx k = 1; k < ctrlsize; ++k) {
270  if (CmidxCTRLcol[k] != first_ctrl_col) {
271  all_ctrl_cols_equal = false;
272  break;
273  }
274  }
275 
276  // at least one control activated, compute the coefficient
277  for (idx n1_ = 0; n1_ < DA; ++n1_) {
278  internal::n2multiidx(n1_, subsyssize, CdimsA, CmidxArow);
279  for (idx k = 0; k < subsyssize; ++k) {
280  Cmidxrow[subsys[k]] = CmidxArow[k];
281  }
282  idx idxrowtmp = internal::multiidx2n(Cmidxrow, N, Cdims);
283 
284  if (all_ctrl_rows_equal) {
285  lhs = Ai[first_ctrl_row](m1_, n1_);
286  } else {
287  lhs = (m1_ == n1_) ? 1 : 0; // identity matrix
288  }
289 
290  for (idx n2_ = 0; n2_ < DA; ++n2_) {
291  internal::n2multiidx(n2_, subsyssize, CdimsA, CmidxAcol);
292  for (idx k = 0; k < subsyssize; ++k) {
293  Cmidxcol[subsys[k]] = CmidxAcol[k];
294  }
295 
296  if (all_ctrl_cols_equal) {
297  rhs = Aidagger[first_ctrl_col](n2_, m2_);
298  } else {
299  rhs = (n2_ == m2_) ? 1 : 0; // identity matrix
300  }
301 
302  idx idxcoltmp = internal::multiidx2n(Cmidxcol, N, Cdims);
303 
304  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
305  }
306  }
307 
308  return std::make_tuple(coeff, idxrow, idxcol);
309  }; /* end coeff_idx_rho */
310 
311  //************ ket ************//
312  if (internal::check_cvector(rstate)) // we have a ket
313  {
314  // check that dims match state vector
315  if (!internal::check_dims_match_cvect(dims, rstate))
316  throw exception::DimsMismatchCvector("qpp::applyCTRL()");
317  if (D == 1)
318  return rstate;
319 
320  dyn_mat<typename Derived1::Scalar> result = rstate;
321 
322 #ifdef WITH_OPENMP_
323 #pragma omp parallel for collapse(2)
324 #endif // WITH_OPENMP_
325  for (idx m = 0; m < DA; ++m)
326  for (idx r = 0; r < DCTRLA_bar; ++r) {
327  if (ctrlsize == 0) // no control
328  {
329  result(coeff_idx_ket(1, m, r).second) =
330  coeff_idx_ket(1, m, r).first;
331  } else
332  for (idx i = 0; i < d; ++i) {
333  result(coeff_idx_ket(i, m, r).second) =
334  coeff_idx_ket(i, m, r).first;
335  }
336  }
337 
338  return result;
339  }
340  //************ density matrix ************//
341  else if (internal::check_square_mat(rstate)) // we have a density operator
342  {
343  // check that dims match state matrix
344  if (!internal::check_dims_match_mat(dims, rstate))
345  throw exception::DimsMismatchMatrix("qpp::applyCTRL()");
346 
347  if (D == 1)
348  return rstate;
349 
350  dyn_mat<typename Derived1::Scalar> result = rstate;
351 
352 #ifdef WITH_OPENMP_
353 #pragma omp parallel for collapse(4)
354 #endif // WITH_OPENMP_
355  for (idx m1 = 0; m1 < DA; ++m1)
356  for (idx r1 = 0; r1 < DCTRLA_bar; ++r1)
357  for (idx m2 = 0; m2 < DA; ++m2)
358  for (idx r2 = 0; r2 < DCTRLA_bar; ++r2)
359  if (ctrlsize == 0) // no control
360  {
361  auto coeff_idxes =
362  coeff_idx_rho(1, m1, r1, 1, m2, r2);
363  result(std::get<1>(coeff_idxes),
364  std::get<2>(coeff_idxes)) =
365  std::get<0>(coeff_idxes);
366  } else {
367  for (idx i1 = 0; i1 < Dctrl; ++i1)
368  for (idx i2 = 0; i2 < Dctrl; ++i2) {
369  auto coeff_idxes =
370  coeff_idx_rho(i1, m1, r1, i2, m2, r2);
371  result(std::get<1>(coeff_idxes),
372  std::get<2>(coeff_idxes)) =
373  std::get<0>(coeff_idxes);
374  }
375  }
376 
377  return result;
378  }
379  //************ Exception: not ket nor density matrix ************//
380  else
381  throw exception::MatrixNotSquareNorCvector("qpp::applyCTRL()");
382 }
383 
399 template <typename Derived1, typename Derived2>
401 applyCTRL(const Eigen::MatrixBase<Derived1>& state,
402  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& ctrl,
403  const std::vector<idx>& subsys, idx d = 2) {
404  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
405  state.derived();
406  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
407 
408  // EXCEPTION CHECKS
409 
410  // check zero size
411  if (!internal::check_nonzero_size(rstate))
412  throw exception::ZeroSize("qpp::applyCTRL()");
413 
414  // check valid dims
415  if (d < 2)
416  throw exception::DimsInvalid("qpp::applyCTRL()");
417  // END EXCEPTION CHECKS
418 
419  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
420  std::vector<idx> dims(N, d); // local dimensions vector
421 
422  return applyCTRL(rstate, rA, ctrl, subsys, dims);
423 }
424 
438 template <typename Derived1, typename Derived2>
440 apply(const Eigen::MatrixBase<Derived1>& state,
441  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& subsys,
442  const std::vector<idx>& dims) {
443  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
444  state.derived();
445  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
446 
447  // EXCEPTION CHECKS
448 
449  // check types
450  if (!std::is_same<typename Derived1::Scalar,
451  typename Derived2::Scalar>::value)
452  throw exception::TypeMismatch("qpp::apply()");
453 
454  // check zero sizes
456  throw exception::ZeroSize("qpp::apply()");
457 
458  // check zero sizes
459  if (!internal::check_nonzero_size(rstate))
460  throw exception::ZeroSize("qpp::apply()");
461 
462  // check square matrix for the gate
464  throw exception::MatrixNotSquare("qpp::apply()");
465 
466  // check that dimension is valid
467  if (!internal::check_dims(dims))
468  throw exception::DimsInvalid("qpp::apply()");
469 
470  // check subsys is valid w.r.t. dims
471  if (!internal::check_subsys_match_dims(subsys, dims))
472  throw exception::SubsysMismatchDims("qpp::apply()");
473 
474  // check that gate matches the dimensions of the subsys
475  std::vector<idx> subsys_dims(subsys.size());
476  for (idx i = 0; i < subsys.size(); ++i)
477  subsys_dims[i] = dims[subsys[i]];
478  if (!internal::check_dims_match_mat(subsys_dims, rA))
479  throw exception::MatrixMismatchSubsys("qpp::apply()");
480  // END EXCEPTION CHECKS
481 
482  //************ ket ************//
483  if (internal::check_cvector(rstate)) // we have a ket
484  {
485  // check that dims match state vector
486  if (!internal::check_dims_match_cvect(dims, rstate))
487  throw exception::DimsMismatchCvector("qpp::apply()");
488 
489  return applyCTRL(rstate, rA, {}, subsys, dims);
490  }
491  //************ density matrix ************//
492  else if (internal::check_square_mat(rstate)) // we have a density operator
493  {
494 
495  // check that dims match state matrix
496  if (!internal::check_dims_match_mat(dims, rstate))
497  throw exception::DimsMismatchMatrix("qpp::apply()");
498 
499  return applyCTRL(rstate, rA, {}, subsys, dims);
500  }
501  //************ Exception: not ket nor density matrix ************//
502  else
503  throw exception::MatrixNotSquareNorCvector("qpp::apply()");
504 }
505 
519 template <typename Derived1, typename Derived2>
521 apply(const Eigen::MatrixBase<Derived1>& state,
522  const Eigen::MatrixBase<Derived2>& A, const std::vector<idx>& subsys,
523  idx d = 2) {
524  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate =
525  state.derived();
526  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
527 
528  // EXCEPTION CHECKS
529 
530  // check zero size
531  if (!internal::check_nonzero_size(rstate))
532  throw exception::ZeroSize("qpp::apply()");
533 
534  // check valid dims
535  if (d < 2)
536  throw exception::DimsInvalid("qpp::apply()");
537  // END EXCEPTION CHECKS
538 
539  idx N = internal::get_num_subsys(static_cast<idx>(rstate.rows()), d);
540  std::vector<idx> dims(N, d); // local dimensions vector
541 
542  return apply(rstate, rA, subsys, dims);
543 }
544 
553 template <typename Derived>
554 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks) {
555  const cmat& rA = A.derived();
556 
557  // EXCEPTION CHECKS
558 
560  throw exception::ZeroSize("qpp::apply()");
562  throw exception::MatrixNotSquare("qpp::apply()");
563  if (Ks.size() == 0)
564  throw exception::ZeroSize("qpp::apply()");
565  if (!internal::check_square_mat(Ks[0]))
566  throw exception::MatrixNotSquare("qpp::apply()");
567  if (Ks[0].rows() != rA.rows())
568  throw exception::DimsMismatchMatrix("qpp::apply()");
569  for (auto&& it : Ks)
570  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
571  throw exception::DimsNotEqual("qpp::apply()");
572  // END EXCEPTION CHECKS
573 
574  cmat result = cmat::Zero(rA.rows(), rA.rows());
575 
576 #ifdef WITH_OPENMP_
577 #pragma omp parallel for
578 #endif // WITH_OPENMP_
579  for (idx i = 0; i < Ks.size(); ++i) {
580 #ifdef WITH_OPENMP_
581 #pragma omp critical
582 #endif // WITH_OPENMP_
583  { result += Ks[i] * rA * adjoint(Ks[i]); }
584  }
585 
586  return result;
587 }
588 
599 template <typename Derived>
600 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
601  const std::vector<idx>& subsys, const std::vector<idx>& dims) {
602  const cmat& rA = A.derived();
603 
604  // EXCEPTION CHECKS
605 
606  // check zero sizes
608  throw exception::ZeroSize("qpp::apply()");
609 
610  // check square matrix for the A
612  throw exception::MatrixNotSquare("qpp::apply()");
613 
614  // check that dimension is valid
615  if (!internal::check_dims(dims))
616  throw exception::DimsInvalid("qpp::apply()");
617 
618  // check that dims match A matrix
619  if (!internal::check_dims_match_mat(dims, rA))
620  throw exception::DimsMismatchMatrix("qpp::apply()");
621 
622  // check subsys is valid w.r.t. dims
623  if (!internal::check_subsys_match_dims(subsys, dims))
624  throw exception::SubsysMismatchDims("qpp::apply()");
625 
626  std::vector<idx> subsys_dims(subsys.size());
627  for (idx i = 0; i < subsys.size(); ++i)
628  subsys_dims[i] = dims[subsys[i]];
629 
630  // check the Kraus operators
631  if (Ks.size() == 0)
632  throw exception::ZeroSize("qpp::apply()");
633  if (!internal::check_square_mat(Ks[0]))
634  throw exception::MatrixNotSquare("qpp::apply()");
635  if (!internal::check_dims_match_mat(subsys_dims, Ks[0]))
636  throw exception::MatrixMismatchSubsys("qpp::apply()");
637  for (auto&& it : Ks)
638  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
639  throw exception::DimsNotEqual("qpp::apply()");
640  // END EXCEPTION CHECKS
641 
642  cmat result = cmat::Zero(rA.rows(), rA.rows());
643 
644  for (idx i = 0; i < Ks.size(); ++i)
645  result += apply(rA, Ks[i], subsys, dims);
646 
647  return result;
648 }
649 
660 template <typename Derived>
661 cmat apply(const Eigen::MatrixBase<Derived>& A, const std::vector<cmat>& Ks,
662  const std::vector<idx>& subsys, idx d = 2) {
663  const cmat& rA = A.derived();
664 
665  // EXCEPTION CHECKS
666 
667  // check zero sizes
669  throw exception::ZeroSize("qpp::apply()");
670 
671  // check valid dims
672  if (d < 2)
673  throw exception::DimsInvalid("qpp::apply()");
674  // END EXCEPTION CHECKS
675 
676  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
677  std::vector<idx> dims(N, d); // local dimensions vector
678 
679  return apply(rA, Ks, subsys, dims);
680 }
681 
693 inline cmat kraus2super(const std::vector<cmat>& Ks) {
694  // EXCEPTION CHECKS
695 
696  if (Ks.size() == 0)
697  throw exception::ZeroSize("qpp::kraus2super()");
698  if (!internal::check_nonzero_size(Ks[0]))
699  throw exception::ZeroSize("qpp::kraus2super()");
700  if (!internal::check_square_mat(Ks[0]))
701  throw exception::MatrixNotSquare("qpp::kraus2super()");
702  for (auto&& it : Ks)
703  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
704  throw exception::DimsNotEqual("qpp::kraus2super()");
705  // END EXCEPTION CHECKS
706 
707  idx D = static_cast<idx>(Ks[0].rows());
708 
709  cmat result(D * D, D * D);
710  cmat MN = cmat::Zero(D, D);
711  bra A = bra::Zero(D);
712  ket B = ket::Zero(D);
713  cmat EMN = cmat::Zero(D, D);
714 
715 #ifdef WITH_OPENMP_
716 #pragma omp parallel for collapse(2)
717 #endif // WITH_OPENMP_
718  for (idx m = 0; m < D; ++m) {
719  for (idx n = 0; n < D; ++n) {
720 #ifdef WITH_OPENMP_
721 #pragma omp critical
722 #endif // WITH_OPENMP_
723  {
724  // compute E(|m><n|)
725  MN(m, n) = 1;
726  for (idx i = 0; i < Ks.size(); ++i)
727  EMN += Ks[i] * MN * adjoint(Ks[i]);
728  MN(m, n) = 0;
729 
730  for (idx a = 0; a < D; ++a) {
731  A(a) = 1;
732  for (idx b = 0; b < D; ++b) {
733  // compute result(ab,mn)=<a|E(|m><n)|b>
734  B(b) = 1;
735  result(a * D + b, m * D + n) =
736  static_cast<cmat>(A * EMN * B).value();
737  B(b) = 0;
738  }
739  A(a) = 0;
740  }
741  EMN = cmat::Zero(D, D);
742  }
743  }
744  }
745 
746  return result;
747 }
748 
764 inline cmat kraus2choi(const std::vector<cmat>& Ks) {
765  // EXCEPTION CHECKS
766 
767  if (Ks.size() == 0)
768  throw exception::ZeroSize("qpp::kraus2choi()");
769  if (!internal::check_nonzero_size(Ks[0]))
770  throw exception::ZeroSize("qpp::kraus2choi()");
771  if (!internal::check_square_mat(Ks[0]))
772  throw exception::MatrixNotSquare("qpp::kraus2choi()");
773  for (auto&& it : Ks)
774  if (it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
775  throw exception::DimsNotEqual("qpp::kraus2choi()");
776  // END EXCEPTION CHECKS
777 
778  idx D = static_cast<idx>(Ks[0].rows());
779 
780  // construct the D x D \sum |jj> vector
781  // (un-normalized maximally entangled state)
782  cmat MES = cmat::Zero(D * D, 1);
783  for (idx a = 0; a < D; ++a)
784  MES(a * D + a) = 1;
785 
786  cmat Omega = MES * adjoint(MES);
787 
788  cmat result = cmat::Zero(D * D, D * D);
789 
790 #ifdef WITH_OPENMP_
791 #pragma omp parallel for
792 #endif // WITH_OPENMP_
793  for (idx i = 0; i < Ks.size(); ++i) {
794 #ifdef WITH_OPENMP_
795 #pragma omp critical
796 #endif // WITH_OPENMP_
797  {
798  result += kron(cmat::Identity(D, D), Ks[i]) * Omega *
799  adjoint(kron(cmat::Identity(D, D), Ks[i]));
800  }
801  }
802 
803  return result;
804 }
805 
819 inline std::vector<cmat> choi2kraus(const cmat& A) {
820  // EXCEPTION CHECKS
821 
823  throw exception::ZeroSize("qpp::choi2kraus()");
825  throw exception::MatrixNotSquare("qpp::choi2kraus()");
826  idx D = internal::get_dim_subsys(A.rows(), 2);
827  // check equal dimensions
828  if (D * D != static_cast<idx>(A.rows()))
829  throw exception::DimsInvalid("qpp::choi2kraus()");
830  // END EXCEPTION CHECKS
831 
832  dmat ev = hevals(A);
833  cmat evec = hevects(A);
834  std::vector<cmat> result;
835 
836  for (idx i = 0; i < D * D; ++i) {
837  if (std::abs(ev(i)) > eps)
838  result.push_back(std::sqrt(std::abs(ev(i))) *
839  reshape(evec.col(i), D, D));
840  }
841 
842  return result;
843 }
844 
852 inline cmat choi2super(const cmat& A) {
853  // EXCEPTION CHECKS
854 
856  throw exception::ZeroSize("qpp::choi2super()");
858  throw exception::MatrixNotSquare("qpp::choi2super()");
859  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
860  // check equal dimensions
861  if (D * D != static_cast<idx>(A.rows()))
862  throw exception::DimsInvalid("qpp::choi2super()");
863  // END EXCEPTION CHECKS
864 
865  cmat result(D * D, D * D);
866 
867 #ifdef WITH_OPENMP_
868 #pragma omp parallel for collapse(4)
869 #endif // WITH_OPENMP_
870  for (idx a = 0; a < D; ++a)
871  for (idx b = 0; b < D; ++b)
872  for (idx m = 0; m < D; ++m)
873  for (idx n = 0; n < D; ++n)
874  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
875 
876  return result;
877 }
878 
886 inline cmat super2choi(const cmat& A) {
887  // EXCEPTION CHECKS
888 
890  throw exception::ZeroSize("qpp::super2choi()");
892  throw exception::MatrixNotSquare("qpp::super2choi()");
893  idx D = internal::get_dim_subsys(static_cast<idx>(A.rows()), 2);
894  // check equal dimensions
895  if (D * D != static_cast<idx>(A.rows()))
896  throw exception::DimsInvalid("qpp::super2choi()");
897  // END EXCEPTION CHECKS
898 
899  cmat result(D * D, D * D);
900 
901 #ifdef WITH_OPENMP_
902 #pragma omp parallel for collapse(4)
903 #endif // WITH_OPENMP_
904  for (idx a = 0; a < D; ++a)
905  for (idx b = 0; b < D; ++b)
906  for (idx m = 0; m < D; ++m)
907  for (idx n = 0; n < D; ++n)
908  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
909 
910  return result;
911 }
912 
926 template <typename Derived>
927 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
928  const std::vector<idx>& dims) {
929  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
930 
931  // EXCEPTION CHECKS
932 
933  // check zero-size
935  throw exception::ZeroSize("qpp::ptrace1()");
936 
937  // check that dims is a valid dimension vector
938  if (!internal::check_dims(dims))
939  throw exception::DimsInvalid("qpp::ptrace1()");
940 
941  // check dims has only 2 elements
942  if (dims.size() != 2)
943  throw exception::NotBipartite("qpp::ptrace1()");
944  // END EXCEPTION CHECKS
945 
946  idx DA = dims[0];
947  idx DB = dims[1];
948 
951 
952  //************ ket ************//
953  if (internal::check_cvector(rA)) // we have a ket
954  {
955  // check that dims match the dimension of A
956  if (!internal::check_dims_match_cvect(dims, rA))
957  throw exception::DimsMismatchCvector("qpp::ptrace1()");
958 
959  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
960  typename Derived::Scalar sum = 0;
961  for (idx m = 0; m < DA; ++m)
962  sum += rA(m * DB + i) * std::conj(rA(m * DB + j));
963 
964  return sum;
965  }; /* end worker */
966 
967 #ifdef WITH_OPENMP_
968 #pragma omp parallel for collapse(2)
969 #endif // WITH_OPENMP_
970  // column major order for speed
971  for (idx j = 0; j < DB; ++j)
972  for (idx i = 0; i < DB; ++i)
973  result(i, j) = worker(i, j);
974 
975  return result;
976  }
977  //************ density matrix ************//
978  else if (internal::check_square_mat(rA)) // we have a density operator
979  {
980  // check that dims match the dimension of A
981  if (!internal::check_dims_match_mat(dims, rA))
982  throw exception::DimsMismatchMatrix("qpp::ptrace1()");
983 
984  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
985  typename Derived::Scalar sum = 0;
986  for (idx m = 0; m < DA; ++m)
987  sum += rA(m * DB + i, m * DB + j);
988 
989  return sum;
990  }; /* end worker */
991 
992 #ifdef WITH_OPENMP_
993 #pragma omp parallel for collapse(2)
994 #endif // WITH_OPENMP_
995  // column major order for speed
996  for (idx j = 0; j < DB; ++j)
997  for (idx i = 0; i < DB; ++i)
998  result(i, j) = worker(i, j);
999 
1000  return result;
1001  }
1002  //************ Exception: not ket nor density matrix ************//
1003  else
1004  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1005 }
1006 
1020 template <typename Derived>
1021 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1022  idx d = 2) {
1023  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1024 
1025  // EXCEPTION CHECKS
1026 
1027  // check zero size
1029  throw exception::ZeroSize("qpp::ptrace1()");
1030 
1031  // check valid dims
1032  if (d == 0)
1033  throw exception::DimsInvalid("qpp::ptrace1()");
1034  // END EXCEPTION CHECKS
1035 
1036  std::vector<idx> dims(2, d); // local dimensions vector
1037 
1038  return ptrace1(A, dims);
1039 }
1040 
1054 template <typename Derived>
1055 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1056  const std::vector<idx>& dims) {
1057  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1058 
1059  // EXCEPTION CHECKS
1060 
1061  // check zero-size
1063  throw exception::ZeroSize("qpp::ptrace2()");
1064 
1065  // check that dims is a valid dimension vector
1066  if (!internal::check_dims(dims))
1067  throw exception::DimsInvalid("qpp::ptrace2()");
1068 
1069  // check dims has only 2 elements
1070  if (dims.size() != 2)
1071  throw exception::NotBipartite("qpp::ptrace2()");
1072  // END EXCEPTION CHECKS
1073 
1074  idx DA = dims[0];
1075  idx DB = dims[1];
1076 
1079 
1080  //************ ket ************//
1081  if (internal::check_cvector(rA)) // we have a ket
1082  {
1083  // check that dims match the dimension of A
1084  if (!internal::check_dims_match_cvect(dims, rA))
1085  throw exception::DimsMismatchCvector("qpp::ptrace2()");
1086 
1087  auto worker = [&](idx i, idx j) noexcept->typename Derived::Scalar {
1088  typename Derived::Scalar sum = 0;
1089  for (idx m = 0; m < DB; ++m)
1090  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1091 
1092  return sum;
1093  }; /* end worker */
1094 
1095 #ifdef WITH_OPENMP_
1096 #pragma omp parallel for collapse(2)
1097 #endif // WITH_OPENMP_
1098  // column major order for speed
1099  for (idx j = 0; j < DA; ++j)
1100  for (idx i = 0; i < DA; ++i)
1101  result(i, j) = worker(i, j);
1102 
1103  return result;
1104  }
1105  //************ density matrix ************//
1106  else if (internal::check_square_mat(rA)) // we have a density operator
1107  {
1108  // check that dims match the dimension of A
1109  if (!internal::check_dims_match_mat(dims, rA))
1110  throw exception::DimsMismatchMatrix("qpp::ptrace2()");
1111 
1112 #ifdef WITH_OPENMP_
1113 #pragma omp parallel for collapse(2)
1114 #endif // WITH_OPENMP_
1115  // column major order for speed
1116  for (idx j = 0; j < DA; ++j)
1117  for (idx i = 0; i < DA; ++i)
1118  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1119 
1120  return result;
1121  }
1122  //************ Exception: not ket nor density matrix ************//
1123  else
1124  throw exception::MatrixNotSquareNorCvector("qpp::ptrace1()");
1125 }
1126 
1140 template <typename Derived>
1141 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1142  idx d = 2) {
1143  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1144 
1145  // EXCEPTION CHECKS
1146 
1147  // check zero size
1149  throw exception::ZeroSize("qpp::ptrace2()");
1150 
1151  // check valid dims
1152  if (d == 0)
1153  throw exception::DimsInvalid("qpp::ptrace2()");
1154  // END EXCEPTION CHECKS
1155 
1156  std::vector<idx> dims(2, d); // local dimensions vector
1157 
1158  return ptrace2(A, dims);
1159 }
1160 
1175 template <typename Derived>
1176 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1177  const std::vector<idx>& subsys,
1178  const std::vector<idx>& dims) {
1179  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1180 
1181  // EXCEPTION CHECKS
1182 
1183  // check zero-size
1185  throw exception::ZeroSize("qpp::ptrace()");
1186 
1187  // check that dims is a valid dimension vector
1188  if (!internal::check_dims(dims))
1189  throw exception::DimsInvalid("qpp::ptrace()");
1190 
1191  // check that subsys are valid
1192  if (!internal::check_subsys_match_dims(subsys, dims))
1193  throw exception::SubsysMismatchDims("qpp::ptrace()");
1194  // END EXCEPTION CHECKS
1195 
1196  idx D = static_cast<idx>(rA.rows());
1197  idx N = dims.size();
1198  idx Nsubsys = subsys.size();
1199  idx Nsubsys_bar = N - Nsubsys;
1200  idx Dsubsys = 1;
1201  for (idx i = 0; i < Nsubsys; ++i)
1202  Dsubsys *= dims[subsys[i]];
1203  idx Dsubsys_bar = D / Dsubsys;
1204 
1205  idx Cdims[maxn];
1206  idx Csubsys[maxn];
1207  idx Cdimssubsys[maxn];
1208  idx Csubsys_bar[maxn];
1209  idx Cdimssubsys_bar[maxn];
1210 
1211  idx Cmidxcolsubsys_bar[maxn];
1212 
1213  std::vector<idx> subsys_bar = complement(subsys, N);
1214  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1215  std::begin(Csubsys_bar));
1216 
1217  for (idx i = 0; i < N; ++i) {
1218  Cdims[i] = dims[i];
1219  }
1220  for (idx i = 0; i < Nsubsys; ++i) {
1221  Csubsys[i] = subsys[i];
1222  Cdimssubsys[i] = dims[subsys[i]];
1223  }
1224  for (idx i = 0; i < Nsubsys_bar; ++i) {
1225  Cdimssubsys_bar[i] = dims[subsys_bar[i]];
1226  }
1227 
1229  dyn_mat<typename Derived::Scalar>(Dsubsys_bar, Dsubsys_bar);
1230 
1231  //************ ket ************//
1232  if (internal::check_cvector(rA)) // we have a ket
1233  {
1234  // check that dims match the dimension of A
1235  if (!internal::check_dims_match_cvect(dims, rA))
1236  throw exception::DimsMismatchCvector("qpp::ptrace()");
1237 
1238  if (subsys.size() == dims.size()) {
1239  result(0, 0) = (adjoint(rA) * rA).value();
1240  return result;
1241  }
1242 
1243  if (subsys.size() == 0)
1244  return rA * adjoint(rA);
1245 
1246  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1247  // use static allocation for speed!
1248 
1249  idx Cmidxrow[maxn];
1250  idx Cmidxcol[maxn];
1251  idx Cmidxrowsubsys_bar[maxn];
1252  idx Cmidxsubsys[maxn];
1253 
1254  /* get the row multi-indexes of the complement */
1255  internal::n2multiidx(i, Nsubsys_bar, Cdimssubsys_bar,
1256  Cmidxrowsubsys_bar);
1257  /* write them in the global row/col multi-indexes */
1258  for (idx k = 0; k < Nsubsys_bar; ++k) {
1259  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1260  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1261  }
1262  typename Derived::Scalar sm = 0;
1263  for (idx a = 0; a < Dsubsys; ++a) {
1264  // get the multi-index over which we do the summation
1265  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1266  // write it into the global row/col multi-indexes
1267  for (idx k = 0; k < Nsubsys; ++k)
1268  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1269  Cmidxsubsys[k];
1270 
1271  // now do the sum
1272  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims)) *
1273  std::conj(rA(internal::multiidx2n(Cmidxcol, N, Cdims)));
1274  }
1275 
1276  return sm;
1277  }; /* end worker */
1278 
1279  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1280  {
1281  // compute the column multi-indexes of the complement
1282  internal::n2multiidx(j, Nsubsys_bar, Cdimssubsys_bar,
1283  Cmidxcolsubsys_bar);
1284 #ifdef WITH_OPENMP_
1285 #pragma omp parallel for
1286 #endif // WITH_OPENMP_
1287  for (idx i = 0; i < Dsubsys_bar; ++i) {
1288  result(i, j) = worker(i);
1289  }
1290  }
1291 
1292  return result;
1293  }
1294  //************ density matrix ************//
1295  else if (internal::check_square_mat(rA)) // we have a density operator
1296  {
1297  // check that dims match the dimension of A
1298  if (!internal::check_dims_match_mat(dims, rA))
1299  throw exception::DimsMismatchMatrix("qpp::ptrace()");
1300 
1301  if (subsys.size() == dims.size()) {
1302  result(0, 0) = rA.trace();
1303  return result;
1304  }
1305 
1306  if (subsys.size() == 0)
1307  return rA;
1308 
1309  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1310  // use static allocation for speed!
1311 
1312  idx Cmidxrow[maxn];
1313  idx Cmidxcol[maxn];
1314  idx Cmidxrowsubsys_bar[maxn];
1315  idx Cmidxsubsys[maxn];
1316 
1317  /* get the row/col multi-indexes of the complement */
1318  internal::n2multiidx(i, Nsubsys_bar, Cdimssubsys_bar,
1319  Cmidxrowsubsys_bar);
1320  /* write them in the global row/col multi-indexes */
1321  for (idx k = 0; k < Nsubsys_bar; ++k) {
1322  Cmidxrow[Csubsys_bar[k]] = Cmidxrowsubsys_bar[k];
1323  Cmidxcol[Csubsys_bar[k]] = Cmidxcolsubsys_bar[k];
1324  }
1325  typename Derived::Scalar sm = 0;
1326  for (idx a = 0; a < Dsubsys; ++a) {
1327  // get the multi-index over which we do the summation
1328  internal::n2multiidx(a, Nsubsys, Cdimssubsys, Cmidxsubsys);
1329  // write it into the global row/col multi-indexes
1330  for (idx k = 0; k < Nsubsys; ++k)
1331  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]] =
1332  Cmidxsubsys[k];
1333 
1334  // now do the sum
1335  sm += rA(internal::multiidx2n(Cmidxrow, N, Cdims),
1336  internal::multiidx2n(Cmidxcol, N, Cdims));
1337  }
1338 
1339  return sm;
1340  }; /* end worker */
1341 
1342  for (idx j = 0; j < Dsubsys_bar; ++j) // column major order for speed
1343  {
1344  // compute the column multi-indexes of the complement
1345  internal::n2multiidx(j, Nsubsys_bar, Cdimssubsys_bar,
1346  Cmidxcolsubsys_bar);
1347 #ifdef WITH_OPENMP_
1348 #pragma omp parallel for
1349 #endif // WITH_OPENMP_
1350  for (idx i = 0; i < Dsubsys_bar; ++i) {
1351  result(i, j) = worker(i);
1352  }
1353  }
1354 
1355  return result;
1356  }
1357  //************ Exception: not ket nor density matrix ************//
1358  else
1359  throw exception::MatrixNotSquareNorCvector("qpp::ptrace()");
1360 }
1361 
1376 template <typename Derived>
1377 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1378  const std::vector<idx>& subsys,
1379  idx d = 2) {
1380  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1381 
1382  // EXCEPTION CHECKS
1383 
1384  // check zero size
1386  throw exception::ZeroSize("qpp::ptrace()");
1387 
1388  // check valid dims
1389  if (d < 2)
1390  throw exception::DimsInvalid("qpp::ptrace()");
1391  // END EXCEPTION CHECKS
1392 
1393  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1394  std::vector<idx> dims(N, d); // local dimensions vector
1395 
1396  return ptrace(rA, subsys, dims);
1397 }
1398 
1412 template <typename Derived>
1414 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& subsys,
1415  const std::vector<idx>& dims) {
1416  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1417 
1418  // EXCEPTION CHECKS
1419 
1420  // check zero-size
1422  throw exception::ZeroSize("qpp::ptranspose()");
1423 
1424  // check that dims is a valid dimension vector
1425  if (!internal::check_dims(dims))
1426  throw exception::DimsInvalid("qpp::ptranspose()");
1427 
1428  // check that subsys are valid
1429  if (!internal::check_subsys_match_dims(subsys, dims))
1430  throw exception::SubsysMismatchDims("qpp::ptranspose()");
1431  // END EXCEPTION CHECKS
1432 
1433  idx D = static_cast<idx>(rA.rows());
1434  idx N = dims.size();
1435  idx Nsubsys = subsys.size();
1436  idx Cdims[maxn];
1437  idx Cmidxcol[maxn];
1438  idx Csubsys[maxn];
1439 
1440  // copy dims in Cdims and subsys in Csubsys
1441  for (idx i = 0; i < N; ++i)
1442  Cdims[i] = dims[i];
1443  for (idx i = 0; i < Nsubsys; ++i)
1444  Csubsys[i] = subsys[i];
1445 
1446  dyn_mat<typename Derived::Scalar> result(D, D);
1447 
1448  //************ ket ************//
1449  if (internal::check_cvector(rA)) // we have a ket
1450  {
1451  // check that dims match the dimension of A
1452  if (!internal::check_dims_match_cvect(dims, rA))
1453  throw exception::DimsMismatchCvector("qpp::ptranspose()");
1454 
1455  if (subsys.size() == dims.size())
1456  return (rA * adjoint(rA)).transpose();
1457 
1458  if (subsys.size() == 0)
1459  return rA * adjoint(rA);
1460 
1461  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1462  // use static allocation for speed!
1463  idx midxcoltmp[maxn];
1464  idx midxrow[maxn];
1465 
1466  for (idx k = 0; k < N; ++k)
1467  midxcoltmp[k] = Cmidxcol[k];
1468 
1469  /* compute the row multi-index */
1470  internal::n2multiidx(i, N, Cdims, midxrow);
1471 
1472  for (idx k = 0; k < Nsubsys; ++k)
1473  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1474 
1475  /* writes the result */
1476  return rA(internal::multiidx2n(midxrow, N, Cdims)) *
1477  std::conj(rA(internal::multiidx2n(midxcoltmp, N, Cdims)));
1478  }; /* end worker */
1479 
1480  for (idx j = 0; j < D; ++j) {
1481  // compute the column multi-index
1482  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1483 
1484 #ifdef WITH_OPENMP_
1485 #pragma omp parallel for
1486 #endif // WITH_OPENMP_
1487  for (idx i = 0; i < D; ++i)
1488  result(i, j) = worker(i);
1489  }
1490 
1491  return result;
1492  }
1493  //************ density matrix ************//
1494  else if (internal::check_square_mat(rA)) // we have a density operator
1495  {
1496  // check that dims match the dimension of A
1497  if (!internal::check_dims_match_mat(dims, rA))
1498  throw exception::DimsMismatchMatrix("qpp::ptranspose()");
1499 
1500  if (subsys.size() == dims.size())
1501  return rA.transpose();
1502 
1503  if (subsys.size() == 0)
1504  return rA;
1505 
1506  auto worker = [&](idx i) noexcept->typename Derived::Scalar {
1507  // use static allocation for speed!
1508  idx midxcoltmp[maxn];
1509  idx midxrow[maxn];
1510 
1511  for (idx k = 0; k < N; ++k)
1512  midxcoltmp[k] = Cmidxcol[k];
1513 
1514  /* compute the row multi-index */
1515  internal::n2multiidx(i, N, Cdims, midxrow);
1516 
1517  for (idx k = 0; k < Nsubsys; ++k)
1518  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1519 
1520  /* writes the result */
1521  return rA(internal::multiidx2n(midxrow, N, Cdims),
1522  internal::multiidx2n(midxcoltmp, N, Cdims));
1523  }; /* end worker */
1524 
1525  for (idx j = 0; j < D; ++j) {
1526  // compute the column multi-index
1527  internal::n2multiidx(j, N, Cdims, Cmidxcol);
1528 
1529 #ifdef WITH_OPENMP_
1530 #pragma omp parallel for
1531 #endif // WITH_OPENMP_
1532  for (idx i = 0; i < D; ++i)
1533  result(i, j) = worker(i);
1534  }
1535 
1536  return result;
1537  }
1538  //************ Exception: not ket nor density matrix ************//
1539  else
1540  throw exception::MatrixNotSquareNorCvector("qpp::ptranspose()");
1541 }
1542 
1556 template <typename Derived>
1558 ptranspose(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& subsys,
1559  idx d = 2) {
1560  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1561 
1562  // EXCEPTION CHECKS
1563 
1564  // check zero size
1566  throw exception::ZeroSize("qpp::ptranspose()");
1567 
1568  // check valid dims
1569  if (d < 2)
1570  throw exception::DimsInvalid("qpp::ptranspose()");
1571  // END EXCEPTION CHECKS
1572 
1573  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1574  std::vector<idx> dims(N, d); // local dimensions vector
1575 
1576  return ptranspose(rA, subsys, dims);
1577 }
1578 
1591 template <typename Derived>
1593 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1594  const std::vector<idx>& dims) {
1595  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1596 
1597  // EXCEPTION CHECKS
1598 
1599  // check zero-size
1601  throw exception::ZeroSize("qpp::syspermute()");
1602 
1603  // check that dims is a valid dimension vector
1604  if (!internal::check_dims(dims))
1605  throw exception::DimsInvalid("qpp::syspermute()");
1606 
1607  // check that we have a valid permutation
1608  if (!internal::check_perm(perm))
1609  throw exception::PermInvalid("qpp::syspermute()");
1610 
1611  // check that permutation match dimensions
1612  if (perm.size() != dims.size())
1613  throw exception::PermMismatchDims("qpp::syspermute()");
1614  // END EXCEPTION CHECKS
1615 
1616  idx D = static_cast<idx>(rA.rows());
1617  idx N = dims.size();
1618 
1620 
1621  //************ ket ************//
1622  if (internal::check_cvector(rA)) // we have a column vector
1623  {
1624  idx Cdims[maxn];
1625  idx Cperm[maxn];
1626 
1627  // check that dims match the dimension of rA
1628  if (!internal::check_dims_match_cvect(dims, rA))
1629  throw exception::DimsMismatchCvector("qpp::syspermute()");
1630 
1631  // copy dims in Cdims and perm in Cperm
1632  for (idx i = 0; i < N; ++i) {
1633  Cdims[i] = dims[i];
1634  Cperm[i] = perm[i];
1635  }
1636  result.resize(D, 1);
1637 
1638  auto worker = [&Cdims, &Cperm, N ](idx i) noexcept->idx {
1639  // use static allocation for speed,
1640  // double the size for matrices reshaped as vectors
1641  idx midx[maxn];
1642  idx midxtmp[maxn];
1643  idx permdims[maxn];
1644 
1645  /* compute the multi-index */
1646  internal::n2multiidx(i, N, Cdims, midx);
1647 
1648  for (idx k = 0; k < N; ++k) {
1649  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1650  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1651  }
1652 
1653  return internal::multiidx2n(midxtmp, N, permdims);
1654  }; /* end worker */
1655 
1656 #ifdef WITH_OPENMP_
1657 #pragma omp parallel for
1658 #endif // WITH_OPENMP_
1659  for (idx i = 0; i < D; ++i)
1660  result(worker(i)) = rA(i);
1661 
1662  return result;
1663  }
1664  //************ density matrix ************//
1665  else if (internal::check_square_mat(rA)) // we have a density operator
1666  {
1667  idx Cdims[2 * maxn];
1668  idx Cperm[2 * maxn];
1669 
1670  // check that dims match the dimension of rA
1671  if (!internal::check_dims_match_mat(dims, rA))
1672  throw exception::DimsMismatchMatrix("qpp::syspermute()");
1673 
1674  // copy dims in Cdims and perm in Cperm
1675  for (idx i = 0; i < N; ++i) {
1676  Cdims[i] = dims[i];
1677  Cdims[i + N] = dims[i];
1678  Cperm[i] = perm[i];
1679  Cperm[i + N] = perm[i] + N;
1680  }
1681  result.resize(D * D, 1);
1682  // map A to a column vector
1684  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1685  const_cast<typename Derived::Scalar*>(rA.data()), D * D, 1);
1686 
1687  auto worker = [&Cdims, &Cperm, N ](idx i) noexcept->idx {
1688  // use static allocation for speed,
1689  // double the size for matrices reshaped as vectors
1690  idx midx[2 * maxn];
1691  idx midxtmp[2 * maxn];
1692  idx permdims[2 * maxn];
1693 
1694  /* compute the multi-index */
1695  internal::n2multiidx(i, 2 * N, Cdims, midx);
1696 
1697  for (idx k = 0; k < 2 * N; ++k) {
1698  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1699  midxtmp[k] = midx[Cperm[k]]; // permuted multi-indexes
1700  }
1701 
1702  return internal::multiidx2n(midxtmp, 2 * N, permdims);
1703  }; /* end worker */
1704 
1705 #ifdef WITH_OPENMP_
1706 #pragma omp parallel for
1707 #endif // WITH_OPENMP_
1708  for (idx i = 0; i < D * D; ++i)
1709  result(worker(i)) = rA(i);
1710 
1711  return reshape(result, D, D);
1712  }
1713  //************ Exception: not ket nor density matrix ************//
1714  else
1715  throw exception::MatrixNotSquareNorCvector("qpp::syspermute()");
1716 }
1717 
1730 template <typename Derived>
1732 syspermute(const Eigen::MatrixBase<Derived>& A, const std::vector<idx>& perm,
1733  idx d = 2) {
1734  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1735 
1736  // EXCEPTION CHECKS
1737 
1738  // check zero size
1740  throw exception::ZeroSize("qpp::syspermute()");
1741 
1742  // check valid dims
1743  if (d < 2)
1744  throw exception::DimsInvalid("qpp::syspermute()");
1745  // END EXCEPTION CHECKS
1746 
1747  idx N = internal::get_num_subsys(static_cast<idx>(rA.rows()), d);
1748  std::vector<idx> dims(N, d); // local dimensions vector
1749 
1750  return syspermute(rA, perm, dims);
1751 }
1752 
1753 } /* namespace qpp */
1754 
1755 #endif /* OPERATIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
Dimensions not equal exception.
Definition: exception.h:282
Dimension(s) mismatch matrix size exception.
Definition: exception.h:298
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:73
bool check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:210
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:819
Eigen::MatrixXd dmat
Real (double precision) dynamic Eigen matrix.
Definition: types.h:69
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:66
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
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:1593
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
bool check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:150
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
Subsystems mismatch dimensions exception.
Definition: exception.h:363
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
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:440
Quantum++ main namespace.
Definition: codes.h:35
bool check_perm(const std::vector< idx > &perm)
Definition: util.h:259
idx get_num_subsys(idx sz, idx d)
Definition: util.h:366
Matrix is not square nor column vector exception.
Definition: exception.h:207
Invalid dimension(s) exception.
Definition: exception.h:267
cmat choi2super(const cmat &A)
Converts Choi matrix to superoperator matrix.
Definition: operations.h:852
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:886
bool check_dims_match_cvect(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:165
Not bi-partite exception.
Definition: exception.h:484
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:927
Invalid permutation exception.
Definition: exception.h:378
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1055
Matrix mismatch subsystems exception.
Definition: exception.h:252
Dimension(s) mismatch column vector size exception.
Definition: exception.h:314
Type mismatch exception.
Definition: exception.h:528
idx get_dim_subsys(idx sz, idx N)
Definition: util.h:378
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:764
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:1176
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
cmat kraus2super(const std::vector< cmat > &Ks)
Superoperator matrix.
Definition: operations.h:693
dyn_mat< typename Derived1::Scalar > applyCTRL(const Eigen::MatrixBase< Derived1 > &state, const Eigen::MatrixBase< Derived2 > &A, const std::vector< idx > &ctrl, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Applies the controlled-gate A to the part subsys of the multi-partite state vector or density matrix ...
Definition: operations.h:55
dyn_mat< typename Derived::Scalar > ptranspose(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &subsys, const std::vector< idx > &dims)
Partial transpose.
Definition: operations.h:1414
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
Permutation mismatch dimensions exception.
Definition: exception.h:394
std::size_t idx
Non-negative integer index.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:147
Object has zero size exception.
Definition: exception.h:132