Quantum++  v0.8.8.2
C++11 quantum computing library
operations.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2016 Vlad Gheorghiu (vgheorgh@gmail.com)
5  *
6  * This file is part of Quantum++.
7  *
8  * Quantum++ is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * Quantum++ is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Quantum++. If not, see <http://www.gnu.org/licenses/>.
20  */
21 
27 #ifndef OPERATIONS_H_
28 #define OPERATIONS_H_
29 
30 // silence g++4.8 bogus warning -Wunused-but-set-variable in lambda functions
31 #if (__GNUC__ && !__clang__)
32 #pragma GCC diagnostic ignored "-Wunused-but-set-variable"
33 #endif
34 
35 namespace qpp
36 {
53 template<typename Derived1, typename Derived2>
55  const Eigen::MatrixBase<Derived1>& state,
56  const Eigen::MatrixBase<Derived2>& A,
57  const std::vector<idx>& ctrl,
58  const std::vector<idx>& subsys,
59  const std::vector<idx>& dims)
60 {
61  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
62  = state.derived();
63  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
64 
65  // EXCEPTION CHECKS
66 
67  // check types
68  if ( !std::is_same<typename Derived1::Scalar,
69  typename Derived2::Scalar>::value )
70  throw Exception("qpp::applyCTRL()", Exception::Type::TYPE_MISMATCH);
71 
72  // check zero sizes
74  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
75 
76  // check zero sizes
77  if ( !internal::_check_nonzero_size(rstate))
78  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
79 
80  // check square matrix for the gate
82  throw Exception("qpp::applyCTRL()",
84 
85  // check that all control subsystems have the same dimension
86  idx d = ctrl.size() > 0 ? dims[ctrl[0]] : 1;
87  for ( idx i = 1; i < ctrl.size(); ++i )
88  if ( dims[ctrl[i]] != d )
89  throw Exception("qpp::applyCTRL()",
91 
92  // check that dimension is valid
93  if ( !internal::_check_dims(dims))
94  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
95 
96  // check subsys is valid w.r.t. dims
97  if ( !internal::_check_subsys_match_dims(subsys, dims))
98  throw Exception("qpp::applyCTRL()",
100 
101  // check that gate matches the dimensions of the subsys
102  std::vector<idx> subsys_dims(subsys.size());
103  for ( idx i = 0; i < subsys.size(); ++i )
104  subsys_dims[i] = dims[subsys[i]];
105  if ( !internal::_check_dims_match_mat(subsys_dims, rA))
106  throw Exception("qpp::applyCTRL()",
108 
109  std::vector<idx> ctrlgate = ctrl; // ctrl + gate subsystem vector
110  ctrlgate.insert(std::end(ctrlgate), std::begin(subsys), std::end(subsys));
111  std::sort(std::begin(ctrlgate), std::end(ctrlgate));
112 
113  // check that ctrl + gate subsystem is valid
114  // with respect to local dimensions
115  if ( !internal::_check_subsys_match_dims(ctrlgate, dims))
116  throw Exception("qpp::applyCTRL()",
118  // END EXCEPTION CHECKS
119 
120  // construct the table of A^i and (A^dagger)^i
121  std::vector<dyn_mat<typename Derived1::Scalar>> Ai;
122  std::vector<dyn_mat<typename Derived1::Scalar>> Aidagger;
123  for ( idx i = 0; i < std::max(d, static_cast<idx>(2)); ++i )
124  {
125  Ai.push_back(powm(rA, i));
126  Aidagger.push_back(powm(adjoint(rA), i));
127  }
128 
129  idx D = static_cast<idx>(rstate.rows()); // total dimension
130  idx n = dims.size(); // total number of subsystems
131  idx ctrlsize = ctrl.size(); // number of ctrl subsystem
132  idx ctrlgatesize = ctrlgate.size(); // number of ctrl+gate subsystems
133  idx subsyssize = subsys.size(); // number of subsystems of the target
134  // dimension of ctrl subsystem
135  idx Dctrl = static_cast<idx>(std::llround(std::pow(d, ctrlsize)));
136  idx DA = static_cast<idx>(rA.rows()); // dimension of gate subsystem
137 
138  idx Cdims[maxn]; // local dimensions
139  idx CdimsA[maxn]; // local dimensions
140  idx CdimsCTRL[maxn]; // local dimensions
141  idx CdimsCTRLAbar[maxn]; // local dimensions
142 
143  // compute the complementary subsystem of ctrlgate w.r.t. dims
144  std::vector<idx> ctrlgatebar = complement(ctrlgate, n);
145  // number of subsystems that are complementary to the ctrl+gate
146  idx ctrlgatebarsize = ctrlgatebar.size();
147 
148  idx DCTRLAbar = 1; // dimension of the rest
149  for ( idx i = 0; i < ctrlgatebarsize; ++i )
150  DCTRLAbar *= dims[ctrlgatebar[i]];
151 
152  for ( idx k = 0; k < n; ++k )
153  Cdims[k] = dims[k];
154  for ( idx k = 0; k < subsyssize; ++k )
155  CdimsA[k] = dims[subsys[k]];
156  for ( idx k = 0; k < ctrlsize; ++k )
157  CdimsCTRL[k] = d;
158  for ( idx k = 0; k < ctrlgatebarsize; ++k )
159  CdimsCTRLAbar[k] = dims[ctrlgatebar[k]];
160 
161  // worker, computes the coefficient and the index for the ket case
162  // used in #pragma omp parallel for collapse
163  auto coeff_idx_ket = [&](idx _i, idx _m, idx _r) noexcept
164  -> std::pair<typename Derived1::Scalar, idx>
165  {
166  idx indx = 0;
167  typename Derived1::Scalar coeff = 0;
168 
169  idx Cmidx[maxn]; // the total multi-index
170  idx CmidxA[maxn]; // the gate part multi-index
171  idx CmidxCTRLAbar[maxn]; // the rest multi-index
172 
173  // compute the index
174 
175  // set the CTRL part
176  for ( idx k = 0; k < ctrlsize; ++k )
177  {
178  Cmidx[ctrl[k]] = _i;
179  }
180 
181  // set the rest
182  internal::_n2multiidx(_r, n - ctrlgatesize,
183  CdimsCTRLAbar, CmidxCTRLAbar);
184  for ( idx k = 0; k < n - ctrlgatesize; ++k )
185  {
186  Cmidx[ctrlgatebar[k]] = CmidxCTRLAbar[k];
187  }
188 
189  // set the A part
190  internal::_n2multiidx(_m, subsyssize, CdimsA, CmidxA);
191  for ( idx k = 0; k < subsyssize; ++k )
192  {
193  Cmidx[subsys[k]] = CmidxA[k];
194  }
195 
196  // we now got the total index
197  indx = internal::_multiidx2n(Cmidx, n, Cdims);
198 
199  // compute the coefficient
200  for ( idx _n = 0; _n < DA; ++_n )
201  {
202  internal::_n2multiidx(_n, subsyssize, CdimsA, CmidxA);
203  for ( idx k = 0; k < subsyssize; ++k )
204  {
205  Cmidx[subsys[k]] = CmidxA[k];
206  }
207  coeff += Ai[_i](_m, _n) *
208  rstate(internal::_multiidx2n(Cmidx, n, Cdims));
209  }
210 
211  return std::make_pair(coeff, indx);
212  }; /* end coeff_idx_ket */
213 
214  // worker, computes the coefficient and the index
215  // for the density matrix case
216  // used in #pragma omp parallel for collapse
217  auto coeff_idx_rho = [&](idx _i1, idx _m1,
218  idx _r1, idx _i2, idx _m2, idx _r2) noexcept
219  -> std::tuple<typename Derived1::Scalar, idx, idx>
220  {
221  idx idxrow = 0;
222  idx idxcol = 0;
223  typename Derived1::Scalar coeff = 0, lhs = 1, rhs = 1;
224 
225  idx Cmidxrow[maxn]; // the total row multi-index
226  idx Cmidxcol[maxn]; // the total col multi-index
227  idx CmidxArow[maxn]; // the gate part row multi-index
228  idx CmidxAcol[maxn]; // the gate part col multi-index
229  idx CmidxCTRLrow[maxn]; // the control row multi-index
230  idx CmidxCTRLcol[maxn]; // the control col multi-index
231  idx CmidxCTRLAbarrow[maxn]; // the rest row multi-index
232  idx CmidxCTRLAbarcol[maxn]; // the rest col multi-index
233 
234  // compute the ket/bra indexes
235 
236  // set the CTRL part
237  internal::_n2multiidx(_i1, ctrlsize,
238  CdimsCTRL, CmidxCTRLrow);
239  internal::_n2multiidx(_i2, ctrlsize,
240  CdimsCTRL, CmidxCTRLcol);
241 
242  for ( idx k = 0; k < ctrlsize; ++k )
243  {
244  Cmidxrow[ctrl[k]] = CmidxCTRLrow[k];
245  Cmidxcol[ctrl[k]] = CmidxCTRLcol[k];
246  }
247 
248  // set the rest
249  internal::_n2multiidx(_r1, n - ctrlgatesize,
250  CdimsCTRLAbar, CmidxCTRLAbarrow);
251  internal::_n2multiidx(_r2, n - ctrlgatesize,
252  CdimsCTRLAbar, CmidxCTRLAbarcol);
253  for ( idx k = 0; k < n - ctrlgatesize; ++k )
254  {
255  Cmidxrow[ctrlgatebar[k]] = CmidxCTRLAbarrow[k];
256  Cmidxcol[ctrlgatebar[k]] = CmidxCTRLAbarcol[k];
257  }
258 
259  // set the A part
260  internal::_n2multiidx(_m1, subsyssize, CdimsA, CmidxArow);
261  internal::_n2multiidx(_m2, subsyssize, CdimsA, CmidxAcol);
262  for ( idx k = 0; k < subsys.size(); ++k )
263  {
264  Cmidxrow[subsys[k]] = CmidxArow[k];
265  Cmidxcol[subsys[k]] = CmidxAcol[k];
266  }
267 
268  // we now got the total row/col indexes
269  idxrow = internal::_multiidx2n(Cmidxrow, n, Cdims);
270  idxcol = internal::_multiidx2n(Cmidxcol, n, Cdims);
271 
272  // check whether all CTRL row and col multi indexes are equal
273  bool all_ctrl_rows_equal = true;
274  bool all_ctrl_cols_equal = true;
275 
276  idx first_ctrl_row, first_ctrl_col;
277  if ( ctrlsize > 0 )
278  {
279  first_ctrl_row = CmidxCTRLrow[0];
280  first_ctrl_col = CmidxCTRLcol[0];
281  }
282  else
283  {
284  first_ctrl_row = first_ctrl_col = 1;
285  }
286 
287  for ( idx k = 1; k < ctrlsize; ++k )
288  {
289  if ( CmidxCTRLrow[k] != first_ctrl_row )
290  {
291  all_ctrl_rows_equal = false;
292  break;
293  }
294  }
295  for ( idx k = 1; k < ctrlsize; ++k )
296  {
297  if ( CmidxCTRLcol[k] != first_ctrl_col )
298  {
299  all_ctrl_cols_equal = false;
300  break;
301  }
302  }
303 
304  // at least one control activated, compute the coefficient
305  for ( idx _n1 = 0; _n1 < DA; ++_n1 )
306  {
307  internal::_n2multiidx(_n1, subsyssize, CdimsA, CmidxArow);
308  for ( idx k = 0; k < subsyssize; ++k )
309  {
310  Cmidxrow[subsys[k]] = CmidxArow[k];
311  }
312  idx idxrowtmp = internal::_multiidx2n(Cmidxrow, n, Cdims);
313 
314  if ( all_ctrl_rows_equal )
315  {
316  lhs = Ai[first_ctrl_row](_m1, _n1);
317  }
318  else
319  {
320  lhs = (_m1 == _n1) ? 1 : 0; // identity matrix
321  }
322 
323  for ( idx _n2 = 0; _n2 < DA; ++_n2 )
324  {
325  internal::_n2multiidx(_n2, subsyssize, CdimsA, CmidxAcol);
326  for ( idx k = 0; k < subsyssize; ++k )
327  {
328  Cmidxcol[subsys[k]] = CmidxAcol[k];
329  }
330 
331  if ( all_ctrl_cols_equal )
332  {
333  rhs = Aidagger[first_ctrl_col](_n2, _m2);
334  }
335  else
336  {
337  rhs = (_n2 == _m2) ? 1 : 0; // identity matrix
338  }
339 
340  idx idxcoltmp = internal::_multiidx2n(Cmidxcol, n, Cdims);
341 
342  coeff += lhs * rstate(idxrowtmp, idxcoltmp) * rhs;
343  }
344  }
345 
346  return std::make_tuple(coeff, idxrow, idxcol);
347  }; /* end coeff_idx_rho */
348 
349  //************ ket ************//
350  if ( internal::_check_cvector(rstate)) // we have a ket
351  {
352  // check that dims match state vector
353  if ( !internal::_check_dims_match_cvect(dims, rstate))
354  throw Exception("qpp::applyCTRL()",
356  if ( D == 1 )
357  return rstate;
358 
359  dyn_mat<typename Derived1::Scalar> result = rstate;
360 
361 #ifdef _WITH_OPENMP_
362 #pragma omp parallel for collapse(2)
363 #endif
364  for ( idx m = 0; m < DA; ++m )
365  for ( idx r = 0; r < DCTRLAbar; ++r )
366  {
367  if ( ctrlsize == 0 ) // no control
368  {
369  result(coeff_idx_ket(1, m, r).second) =
370  coeff_idx_ket(1, m, r).first;
371  }
372  else
373  for ( idx i = 0; i < d; ++i )
374  {
375  result(coeff_idx_ket(i, m, r).second) =
376  coeff_idx_ket(i, m, r).first;
377  }
378  }
379 
380  return result;
381  }
382  //************ density matrix ************//
383  else if ( internal::_check_square_mat(rstate)) // we have a density operator
384  {
385  // check that dims match state matrix
386  if ( !internal::_check_dims_match_mat(dims, rstate))
387  throw Exception("qpp::applyCTRL()",
389 
390  if ( D == 1 )
391  return rstate;
392 
393  dyn_mat<typename Derived1::Scalar> result = rstate;
394 
395 #ifdef _WITH_OPENMP_
396 #pragma omp parallel for collapse(4)
397 #endif
398  for ( idx m1 = 0; m1 < DA; ++m1 )
399  for ( idx r1 = 0; r1 < DCTRLAbar; ++r1 )
400  for ( idx m2 = 0; m2 < DA; ++m2 )
401  for ( idx r2 = 0; r2 < DCTRLAbar; ++r2 )
402  if ( ctrlsize == 0 ) // no control
403  {
404  auto coeff_idxes = coeff_idx_rho(1, m1, r1,
405  1, m2, r2);
406  result(std::get<1>(coeff_idxes),
407  std::get<2>(coeff_idxes)) =
408  std::get<0>(coeff_idxes);
409  }
410  else
411  {
412  for ( idx i1 = 0; i1 < Dctrl; ++i1 )
413  for ( idx i2 = 0; i2 < Dctrl; ++i2 )
414  {
415  auto coeff_idxes = coeff_idx_rho(
416  i1, m1, r1,
417  i2, m2, r2);
418  result(std::get<1>(coeff_idxes),
419  std::get<2>(coeff_idxes)) =
420  std::get<0>(coeff_idxes);
421  }
422  }
423 
424  return result;
425  }
426  //************ Exception: not ket nor density matrix ************//
427  else
428  throw Exception("qpp::applyCTRL()",
430 }
431 
447 template<typename Derived1, typename Derived2>
449  const Eigen::MatrixBase<Derived1>& state,
450  const Eigen::MatrixBase<Derived2>& A,
451  const std::vector<idx>& ctrl,
452  const std::vector<idx>& subsys,
453  idx d = 2)
454 {
455  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
456  = state.derived();
457  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
458 
459  // EXCEPTION CHECKS
460 
461  // check zero size
462  if ( !internal::_check_nonzero_size(rstate))
463  throw Exception("qpp::applyCTRL()", Exception::Type::ZERO_SIZE);
464 
465  // check valid dims
466  if ( d == 0 )
467  throw Exception("qpp::applyCTRL()", Exception::Type::DIMS_INVALID);
468  // END EXCEPTION CHECKS
469 
470  idx n =
471  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
472  std::log2(d)));
473  std::vector<idx> dims(n, d); // local dimensions vector
474 
475  return applyCTRL(rstate, rA, ctrl, subsys, dims);
476 }
477 
491 template<typename Derived1, typename Derived2>
493  const Eigen::MatrixBase<Derived1>& state,
494  const Eigen::MatrixBase<Derived2>& A,
495  const std::vector<idx>& subsys,
496  const std::vector<idx>& dims)
497 {
498  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
499  = state.derived();
500  const dyn_mat<typename Derived2::Scalar>& rA = A.derived();
501 
502  // EXCEPTION CHECKS
503 
504  // check types
505  if ( !std::is_same<typename Derived1::Scalar,
506  typename Derived2::Scalar>::value )
507  throw Exception("qpp::apply()", Exception::Type::TYPE_MISMATCH);
508 
509  // check zero sizes
511  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
512 
513  // check zero sizes
514  if ( !internal::_check_nonzero_size(rstate))
515  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
516 
517  // check square matrix for the gate
518  if ( !internal::_check_square_mat(rA))
519  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
520 
521  // check that dimension is valid
522  if ( !internal::_check_dims(dims))
523  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
524 
525  // check subsys is valid w.r.t. dims
526  if ( !internal::_check_subsys_match_dims(subsys, dims))
527  throw Exception("qpp::apply()", Exception::Type::SUBSYS_MISMATCH_DIMS);
528 
529  // check that gate matches the dimensions of the subsys
530  std::vector<idx> subsys_dims(subsys.size());
531  for ( idx i = 0; i < subsys.size(); ++i )
532  subsys_dims[i] = dims[subsys[i]];
533  if ( !internal::_check_dims_match_mat(subsys_dims, rA))
534  throw Exception("qpp::apply()",
536  // END EXCEPTION CHECKS
537 
538  //************ ket ************//
539  if ( internal::_check_cvector(rstate)) // we have a ket
540  {
541  // check that dims match state vector
542  if ( !internal::_check_dims_match_cvect(dims, rstate))
543  throw Exception("qpp::apply()",
545 
546  return applyCTRL(rstate, rA, {}, subsys, dims);
547  }
548  //************ density matrix ************//
549  else if ( internal::_check_square_mat(rstate)) // we have a density operator
550  {
551 
552  // check that dims match state matrix
553  if ( !internal::_check_dims_match_mat(dims, rstate))
554  throw Exception("qpp::apply()",
556 
557  return applyCTRL(rstate, rA, {}, subsys, dims);
558  }
559  //************ Exception: not ket nor density matrix ************//
560  else
561  throw Exception("qpp::apply()",
563 }
564 
578 template<typename Derived1, typename Derived2>
580  const Eigen::MatrixBase<Derived1>& state,
581  const Eigen::MatrixBase<Derived2>& A,
582  const std::vector<idx>& subsys,
583  idx d = 2)
584 {
585  const typename Eigen::MatrixBase<Derived1>::EvalReturnType& rstate
586  = state.derived();
587  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
588 
589  // EXCEPTION CHECKS
590 
591  // check zero size
592  if ( !internal::_check_nonzero_size(rstate))
593  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
594 
595  // check valid dims
596  if ( d == 0 )
597  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
598  // END EXCEPTION CHECKS
599 
600  idx n =
601  static_cast<idx>(std::llround(std::log2(rstate.rows()) /
602  std::log2(d)));
603  std::vector<idx> dims(n, d); // local dimensions vector
604 
605  return apply(rstate, rA, subsys, dims);
606 }
607 
616 template<typename Derived>
617 cmat apply(const Eigen::MatrixBase<Derived>& rho,
618  const std::vector<cmat>& Ks)
619 {
620  const cmat& rrho = rho.derived();
621 
622  // EXCEPTION CHECKS
623 
624  if ( !internal::_check_nonzero_size(rrho))
625  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
626  if ( !internal::_check_square_mat(rrho))
627  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
628  if ( Ks.size() == 0 )
629  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
630  if ( !internal::_check_square_mat(Ks[0]))
631  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
632  if ( Ks[0].rows() != rrho.rows())
633  throw Exception("qpp::apply()",
635  for ( auto&& it : Ks )
636  if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
637  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
638  // END EXCEPTION CHECKS
639 
640  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
641 
642 #ifdef _WITH_OPENMP_
643 #pragma omp parallel for
644 #endif
645  for ( idx i = 0; i < Ks.size(); ++i )
646  {
647 #ifdef _WITH_OPENMP_
648 #pragma omp critical
649 #endif
650  {
651  result += Ks[i] * rrho * adjoint(Ks[i]);
652  }
653  }
654 
655  return result;
656 }
657 
668 template<typename Derived>
669 cmat apply(const Eigen::MatrixBase<Derived>& rho,
670  const std::vector<cmat>& Ks,
671  const std::vector<idx>& subsys,
672  const std::vector<idx>& dims)
673 {
674  const cmat& rrho = rho.derived();
675 
676  // EXCEPTION CHECKS
677 
678  // check zero sizes
679  if ( !internal::_check_nonzero_size(rrho))
680  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
681 
682  // check square matrix for the rho
683  if ( !internal::_check_square_mat(rrho))
684  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
685 
686  // check that dimension is valid
687  if ( !internal::_check_dims(dims))
688  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
689 
690  // check that dims match rho matrix
691  if ( !internal::_check_dims_match_mat(dims, rrho))
692  throw Exception("qpp::apply()",
694 
695  // check subsys is valid w.r.t. dims
696  if ( !internal::_check_subsys_match_dims(subsys, dims))
697  throw Exception("qpp::apply()",
699 
700  std::vector<idx> subsys_dims(subsys.size());
701  for ( idx i = 0; i < subsys.size(); ++i )
702  subsys_dims[i] = dims[subsys[i]];
703 
704  // check the Kraus operators
705  if ( Ks.size() == 0 )
706  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
707  if ( !internal::_check_square_mat(Ks[0]))
708  throw Exception("qpp::apply()", Exception::Type::MATRIX_NOT_SQUARE);
709  if ( !internal::_check_dims_match_mat(subsys_dims, Ks[0]))
710  throw Exception("qpp::apply()",
712  for ( auto&& it : Ks )
713  if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
714  throw Exception("qpp::apply()", Exception::Type::DIMS_NOT_EQUAL);
715  // END EXCEPTION CHECKS
716 
717  cmat result = cmat::Zero(rrho.rows(), rrho.rows());
718 
719  for ( idx i = 0; i < Ks.size(); ++i )
720  result += apply(rrho, Ks[i], subsys, dims);
721 
722  return result;
723 }
724 
735 template<typename Derived>
736 cmat apply(const Eigen::MatrixBase<Derived>& rho,
737  const std::vector<cmat>& Ks,
738  const std::vector<idx>& subsys,
739  idx d = 2)
740 {
741  const cmat& rrho = rho.derived();
742 
743  // EXCEPTION CHECKS
744 
745  // check zero sizes
746  if ( !internal::_check_nonzero_size(rrho))
747  throw Exception("qpp::apply()", Exception::Type::ZERO_SIZE);
748 
749  // check valid dims
750  if ( d == 0 )
751  throw Exception("qpp::apply()", Exception::Type::DIMS_INVALID);
752  // END EXCEPTION CHECKS
753 
754  idx n =
755  static_cast<idx>(std::llround(std::log2(rrho.rows()) /
756  std::log2(d)));
757  std::vector<idx> dims(n, d); // local dimensions vector
758 
759  return apply(rrho, Ks, subsys, dims);
760 }
761 
773 inline cmat kraus2super(const std::vector<cmat>& Ks)
774 {
775  // EXCEPTION CHECKS
776 
777  if ( Ks.size() == 0 )
778  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
779  if ( !internal::_check_nonzero_size(Ks[0]))
780  throw Exception("qpp::kraus2super()", Exception::Type::ZERO_SIZE);
781  if ( !internal::_check_square_mat(Ks[0]))
782  throw Exception("qpp::kraus2super()",
784  for ( auto&& it : Ks )
785  if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
786  throw Exception("qpp::kraus2super()",
788  // END EXCEPTION CHECKS
789 
790  idx D = static_cast<idx>(Ks[0].rows());
791 
792  cmat result(D * D, D * D);
793  cmat MN = cmat::Zero(D, D);
794  bra A = bra::Zero(D);
795  ket B = ket::Zero(D);
796  cmat EMN = cmat::Zero(D, D);
797 
798 #ifdef _WITH_OPENMP_
799 #pragma omp parallel for collapse(2)
800 #endif
801  for ( idx m = 0; m < D; ++m )
802  {
803  for ( idx n = 0; n < D; ++n )
804  {
805 #ifdef _WITH_OPENMP_
806 #pragma omp critical
807 #endif
808  {
809  // compute E(|m><n|)
810  MN(m, n) = 1;
811  for ( idx i = 0; i < Ks.size(); ++i )
812  EMN += Ks[i] * MN * adjoint(Ks[i]);
813  MN(m, n) = 0;
814 
815  for ( idx a = 0; a < D; ++a )
816  {
817  A(a) = 1;
818  for ( idx b = 0; b < D; ++b )
819  {
820  // compute result(ab,mn)=<a|E(|m><n)|b>
821  B(b) = 1;
822  result(a * D + b, m * D + n) =
823  static_cast<cmat>(A * EMN * B).value();
824  B(b) = 0;
825  }
826  A(a) = 0;
827  }
828  EMN = cmat::Zero(D, D);
829  }
830  }
831  }
832 
833  return result;
834 }
835 
851 inline cmat kraus2choi(const std::vector<cmat>& Ks)
852 {
853  // EXCEPTION CHECKS
854 
855  if ( Ks.size() == 0 )
856  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
857  if ( !internal::_check_nonzero_size(Ks[0]))
858  throw Exception("qpp::kraus2choi()", Exception::Type::ZERO_SIZE);
859  if ( !internal::_check_square_mat(Ks[0]))
860  throw Exception("qpp::kraus2choi()",
862  for ( auto&& it : Ks )
863  if ( it.rows() != Ks[0].rows() || it.cols() != Ks[0].rows())
864  throw Exception("qpp::kraus2choi()",
866  // END EXCEPTION CHECKS
867 
868  idx D = static_cast<idx>(Ks[0].rows());
869 
870  // construct the D x D \sum |jj> vector
871  // (un-normalized maximally entangled state)
872  cmat MES = cmat::Zero(D * D, 1);
873  for ( idx a = 0; a < D; ++a )
874  MES(a * D + a) = 1;
875 
876  cmat Omega = MES * adjoint(MES);
877 
878  cmat result = cmat::Zero(D * D, D * D);
879 
880 #ifdef _WITH_OPENMP_
881 #pragma omp parallel for
882 #endif
883  for ( idx i = 0; i < Ks.size(); ++i )
884  {
885 #ifdef _WITH_OPENMP_
886 #pragma omp critical
887 #endif
888  {
889  result += kron(cmat::Identity(D, D), Ks[i]) * Omega
890  * adjoint(kron(cmat::Identity(D, D), Ks[i]));
891  }
892  }
893 
894  return result;
895 }
896 
910 inline std::vector<cmat> choi2kraus(const cmat& A)
911 {
912  // EXCEPTION CHECKS
913 
915  throw Exception("qpp::choi2kraus()", Exception::Type::ZERO_SIZE);
917  throw Exception("qpp::choi2kraus()",
919  idx D = static_cast<idx>(std::llround(
920  std::sqrt(static_cast<double>(A.rows()))));
921  if ( D * D != static_cast<idx>(A.rows()))
922  throw Exception("qpp::choi2kraus()", Exception::Type::DIMS_INVALID);
923  // END EXCEPTION CHECKS
924 
925  dmat ev = hevals(A);
926  cmat evec = hevects(A);
927  std::vector<cmat> result;
928 
929  for ( idx i = 0; i < D * D; ++i )
930  {
931  if ( std::abs(ev(i)) > eps )
932  result.push_back(
933  std::sqrt(std::abs(ev(i))) * reshape(evec.col(i), D, D));
934  }
935 
936  return result;
937 }
938 
946 inline cmat choi2super(const cmat& A)
947 {
948  // EXCEPTION CHECKS
949 
951  throw Exception("qpp::choi2super()", Exception::Type::ZERO_SIZE);
953  throw Exception("qpp::choi2super()",
955  idx D = static_cast<idx>(std::llround(
956  std::sqrt(static_cast<double>(A.rows()))));
957  if ( D * D != static_cast<idx>(A.rows()))
958  throw Exception("qpp::choi2super()", Exception::Type::DIMS_INVALID);
959  // END EXCEPTION CHECKS
960 
961  cmat result(D * D, D * D);
962 
963 #ifdef _WITH_OPENMP_
964 #pragma omp parallel for collapse(4)
965 #endif
966  for ( idx a = 0; a < D; ++a )
967  for ( idx b = 0; b < D; ++b )
968  for ( idx m = 0; m < D; ++m )
969  for ( idx n = 0; n < D; ++n )
970  result(a * D + b, m * D + n) = A(m * D + a, n * D + b);
971 
972  return result;
973 }
974 
982 inline cmat super2choi(const cmat& A)
983 {
984  // EXCEPTION CHECKS
985 
987  throw Exception("qpp::super2choi()", Exception::Type::ZERO_SIZE);
989  throw Exception("qpp::super2choi()",
991  idx D = static_cast<idx>(std::llround(
992  std::sqrt(static_cast<double>(A.rows()))));
993  if ( D * D != static_cast<idx>(A.rows()))
994  throw Exception("qpp::super2choi()", Exception::Type::DIMS_INVALID);
995  // END EXCEPTION CHECKS
996 
997  cmat result(D * D, D * D);
998 
999 #ifdef _WITH_OPENMP_
1000 #pragma omp parallel for collapse(4)
1001 #endif
1002  for ( idx a = 0; a < D; ++a )
1003  for ( idx b = 0; b < D; ++b )
1004  for ( idx m = 0; m < D; ++m )
1005  for ( idx n = 0; n < D; ++n )
1006  result(m * D + a, n * D + b) = A(a * D + b, m * D + n);
1007 
1008  return result;
1009 }
1010 
1025 template<typename Derived>
1026 dyn_mat<typename Derived::Scalar> ptrace1(const Eigen::MatrixBase<Derived>& A,
1027  const std::vector<idx>& dims)
1028 {
1029  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1030  = A.derived();
1031 
1032  // EXCEPTION CHECKS
1033 
1034  // check zero-size
1036  throw Exception("qpp::ptrace1()", Exception::Type::ZERO_SIZE);
1037 
1038  // check that dims is a valid dimension vector
1039  if ( !internal::_check_dims(dims))
1040  throw Exception("qpp::ptrace1()", Exception::Type::DIMS_INVALID);
1041 
1042  // check dims has only 2 elements
1043  if ( dims.size() != 2 )
1044  throw Exception("qpp::ptrace1()", Exception::Type::NOT_BIPARTITE);
1045  // END EXCEPTION CHECKS
1046 
1047  idx DA = dims[0];
1048  idx DB = dims[1];
1049 
1052 
1053  //************ ket ************//
1054  if ( internal::_check_cvector(rA)) // we have a ket
1055  {
1056  // check that dims match the dimension of A
1057  if ( !internal::_check_dims_match_cvect(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) * std::conj(rA(m * DB + j));
1067 
1068  return sum;
1069  }; /* end worker */
1070 
1071 #ifdef _WITH_OPENMP_
1072 #pragma omp parallel for collapse(2)
1073 #endif
1074  for ( idx j = 0; j < DB; ++j ) // column major order for speed
1075  for ( idx i = 0; i < DB; ++i )
1076  result(i, j) = worker(i, j);
1077 
1078  return result;
1079  }
1080  //************ density matrix ************//
1081  else if ( internal::_check_square_mat(rA)) // we have a density operator
1082  {
1083  // check that dims match the dimension of A
1084  if ( !internal::_check_dims_match_mat(dims, rA))
1085  throw Exception("qpp::ptrace1()",
1087 
1088  auto worker = [=](idx i, idx j) noexcept
1089  -> typename Derived::Scalar
1090  {
1091  typename Derived::Scalar sum = 0;
1092  for ( idx m = 0; m < DA; ++m )
1093  sum += rA(m * DB + i, m * DB + j);
1094 
1095  return sum;
1096  }; /* end worker */
1097 
1098 #ifdef _WITH_OPENMP_
1099 #pragma omp parallel for collapse(2)
1100 #endif
1101  for ( idx j = 0; j < DB; ++j ) // column major order for speed
1102  for ( idx i = 0; i < DB; ++i )
1103  result(i, j) = worker(i, j);
1104 
1105  return result;
1106  }
1107  //************ Exception: not ket nor density matrix ************//
1108  else
1109  throw Exception("qpp::ptrace1()",
1111 }
1112 
1127 template<typename Derived>
1128 dyn_mat<typename Derived::Scalar> ptrace2(const Eigen::MatrixBase<Derived>& A,
1129  const std::vector<idx>& dims)
1130 {
1131  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1132  = A.derived();
1133 
1134  // EXCEPTION CHECKS
1135 
1136  // check zero-size
1138  throw Exception("qpp::ptrace2()", Exception::Type::ZERO_SIZE);
1139 
1140  // check that dims is a valid dimension vector
1141  if ( !internal::_check_dims(dims))
1142  throw Exception("qpp::ptrace2()", Exception::Type::DIMS_INVALID);
1143 
1144  // check dims has only 2 elements
1145  if ( dims.size() != 2 )
1146  throw Exception("qpp::ptrace2()", Exception::Type::NOT_BIPARTITE);
1147  // END EXCEPTION CHECKS
1148 
1149  idx DA = dims[0];
1150  idx DB = dims[1];
1151 
1154 
1155  //************ ket ************//
1156  if ( internal::_check_cvector(rA)) // we have a ket
1157  {
1158  // check that dims match the dimension of A
1159  if ( !internal::_check_dims_match_cvect(dims, rA))
1160  throw Exception("qpp::ptrace2()",
1162 
1163  auto worker = [=](idx i, idx j) noexcept
1164  -> typename Derived::Scalar
1165  {
1166  typename Derived::Scalar sum = 0;
1167  for ( idx m = 0; m < DB; ++m )
1168  sum += rA(i * DB + m) * std::conj(rA(j * DB + m));
1169 
1170  return sum;
1171  }; /* end worker */
1172 
1173 #ifdef _WITH_OPENMP_
1174 #pragma omp parallel for collapse(2)
1175 #endif
1176  for ( idx j = 0; j < DA; ++j ) // column major order for speed
1177  for ( idx i = 0; i < DA; ++i )
1178  result(i, j) = worker(i, j);
1179 
1180  return result;
1181  }
1182  //************ density matrix ************//
1183  else if ( internal::_check_square_mat(rA)) // we have a density operator
1184  {
1185  // check that dims match the dimension of A
1186  if ( !internal::_check_dims_match_mat(dims, rA))
1187  throw Exception("qpp::ptrace2()",
1189 
1190 #ifdef _WITH_OPENMP_
1191 #pragma omp parallel for collapse(2)
1192 #endif
1193  for ( idx j = 0; j < DA; ++j ) // column major order for speed
1194  for ( idx i = 0; i < DA; ++i )
1195  result(i, j) = trace(rA.block(i * DB, j * DB, DB, DB));
1196 
1197  return result;
1198  }
1199  //************ Exception: not ket nor density matrix ************//
1200  else
1201  throw Exception("qpp::ptrace1()",
1203 }
1204 
1219 template<typename Derived>
1220 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1221  const std::vector<idx>& subsys,
1222  const std::vector<idx>& dims)
1223 {
1224  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA
1225  = A.derived();
1226 
1227  // EXCEPTION CHECKS
1228 
1229  // check zero-size
1231  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1232 
1233  // check that dims is a valid dimension vector
1234  if ( !internal::_check_dims(dims))
1235  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1236 
1237  // check that subsys are valid
1238  if ( !internal::_check_subsys_match_dims(subsys, dims))
1239  throw Exception("qpp::ptrace()",
1241  // END EXCEPTION CHECKS
1242 
1243  idx D = static_cast<idx>(rA.rows());
1244  idx n = dims.size();
1245  idx nsubsys = subsys.size();
1246  idx nsubsysbar = n - nsubsys;
1247  idx Dsubsys = 1;
1248  for ( idx i = 0; i < nsubsys; ++i )
1249  Dsubsys *= dims[subsys[i]];
1250  idx Dsubsysbar = D / Dsubsys;
1251 
1252  idx Cdims[maxn];
1253  idx Csubsys[maxn];
1254  idx Cdimssubsys[maxn];
1255  idx Csubsysbar[maxn];
1256  idx Cdimssubsysbar[maxn];
1257 
1258  idx Cmidxcolsubsysbar[maxn];
1259 
1260  std::vector<idx> subsys_bar = complement(subsys, n);
1261  std::copy(std::begin(subsys_bar), std::end(subsys_bar),
1262  std::begin(Csubsysbar));
1263 
1264  for ( idx i = 0; i < n; ++i )
1265  {
1266  Cdims[i] = dims[i];
1267  }
1268  for ( idx i = 0; i < nsubsys; ++i )
1269  {
1270  Csubsys[i] = subsys[i];
1271  Cdimssubsys[i] = dims[subsys[i]];
1272  }
1273  for ( idx i = 0; i < nsubsysbar; ++i )
1274  {
1275  Cdimssubsysbar[i] = dims[subsys_bar[i]];
1276  }
1277 
1279  dyn_mat<typename Derived::Scalar>(Dsubsysbar, Dsubsysbar);
1280 
1281  //************ ket ************//
1282  if ( internal::_check_cvector(rA)) // we have a ket
1283  {
1284  // check that dims match the dimension of A
1285  if ( !internal::_check_dims_match_cvect(dims, rA))
1286  throw Exception("qpp::ptrace()",
1288 
1289  if ( subsys.size() == dims.size())
1290  {
1291  result(0, 0) = (adjoint(rA) * rA).value();
1292  return result;
1293  }
1294 
1295  if ( subsys.size() == 0 )
1296  return rA * adjoint(rA);
1297 
1298  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1299  -> typename Derived::Scalar
1300  {
1301  // use static allocation for speed!
1302 
1303  idx Cmidxrow[maxn];
1304  idx Cmidxcol[maxn];
1305  idx Cmidxrowsubsysbar[maxn];
1306  idx Cmidxsubsys[maxn];
1307 
1308  /* get the row multi-indexes of the complement */
1309  internal::_n2multiidx(i, nsubsysbar,
1310  Cdimssubsysbar, Cmidxrowsubsysbar);
1311  /* write them in the global row/col multi-indexes */
1312  for ( idx k = 0; k < nsubsysbar; ++k )
1313  {
1314  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1315  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1316  }
1317  typename Derived::Scalar sm = 0;
1318  for ( idx a = 0; a < Dsubsys; ++a )
1319  {
1320  // get the multi-index over which we do the summation
1321  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1322  // write it into the global row/col multi-indexes
1323  for ( idx k = 0; k < nsubsys; ++k )
1324  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1325  = Cmidxsubsys[k];
1326 
1327  // now do the sum
1328  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims)) *
1329  std::conj(rA(internal::_multiidx2n(Cmidxcol, n,
1330  Cdims)));
1331  }
1332 
1333  return sm;
1334  }; /* end worker */
1335 
1336  for ( idx j = 0; j < Dsubsysbar; ++j ) // column major order for speed
1337  {
1338  // compute the column multi-indexes of the complement
1339  internal::_n2multiidx(j, nsubsysbar,
1340  Cdimssubsysbar, Cmidxcolsubsysbar);
1341 #ifdef _WITH_OPENMP_
1342 #pragma omp parallel for
1343 #endif
1344  for ( idx i = 0; i < Dsubsysbar; ++i )
1345  {
1346  result(i, j) = worker(i);
1347  }
1348  }
1349 
1350  return result;
1351  }
1352  //************ density matrix ************//
1353  else if ( internal::_check_square_mat(rA)) // we have a density operator
1354  {
1355  // check that dims match the dimension of A
1356  if ( !internal::_check_dims_match_mat(dims, rA))
1357  throw Exception("qpp::ptrace()",
1359 
1360  if ( subsys.size() == dims.size())
1361  {
1362  result(0, 0) = rA.trace();
1363  return result;
1364  }
1365 
1366  if ( subsys.size() == 0 )
1367  return rA;
1368 
1369  auto worker = [=, &Cmidxcolsubsysbar](idx i) noexcept
1370  -> typename Derived::Scalar
1371  {
1372  // use static allocation for speed!
1373 
1374  idx Cmidxrow[maxn];
1375  idx Cmidxcol[maxn];
1376  idx Cmidxrowsubsysbar[maxn];
1377  idx Cmidxsubsys[maxn];
1378 
1379  /* get the row/col multi-indexes of the complement */
1380  internal::_n2multiidx(i, nsubsysbar,
1381  Cdimssubsysbar, Cmidxrowsubsysbar);
1382  /* write them in the global row/col multi-indexes */
1383  for ( idx k = 0; k < nsubsysbar; ++k )
1384  {
1385  Cmidxrow[Csubsysbar[k]] = Cmidxrowsubsysbar[k];
1386  Cmidxcol[Csubsysbar[k]] = Cmidxcolsubsysbar[k];
1387  }
1388  typename Derived::Scalar sm = 0;
1389  for ( idx a = 0; a < Dsubsys; ++a )
1390  {
1391  // get the multi-index over which we do the summation
1392  internal::_n2multiidx(a, nsubsys, Cdimssubsys, Cmidxsubsys);
1393  // write it into the global row/col multi-indexes
1394  for ( idx k = 0; k < nsubsys; ++k )
1395  Cmidxrow[Csubsys[k]] = Cmidxcol[Csubsys[k]]
1396  = Cmidxsubsys[k];
1397 
1398  // now do the sum
1399  sm += rA(internal::_multiidx2n(Cmidxrow, n, Cdims),
1400  internal::_multiidx2n(Cmidxcol, n, Cdims));
1401  }
1402 
1403  return sm;
1404  }; /* end worker */
1405 
1406  for ( idx j = 0; j < Dsubsysbar; ++j ) // column major order for speed
1407  {
1408  // compute the column multi-indexes of the complement
1409  internal::_n2multiidx(j, nsubsysbar,
1410  Cdimssubsysbar, Cmidxcolsubsysbar);
1411 #ifdef _WITH_OPENMP_
1412 #pragma omp parallel for
1413 #endif
1414  for ( idx i = 0; i < Dsubsysbar; ++i )
1415  {
1416  result(i, j) = worker(i);
1417  }
1418  }
1419 
1420  return result;
1421  }
1422  //************ Exception: not ket nor density matrix ************//
1423  else
1424  throw Exception("qpp::ptrace()",
1426 }
1427 
1442 template<typename Derived>
1443 dyn_mat<typename Derived::Scalar> ptrace(const Eigen::MatrixBase<Derived>& A,
1444  const std::vector<idx>& subsys,
1445  idx d = 2)
1446 {
1447  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1448 
1449  // EXCEPTION CHECKS
1450 
1451  // check zero size
1453  throw Exception("qpp::ptrace()", Exception::Type::ZERO_SIZE);
1454 
1455  // check valid dims
1456  if ( d == 0 )
1457  throw Exception("qpp::ptrace()", Exception::Type::DIMS_INVALID);
1458  // END EXCEPTION CHECKS
1459 
1460  idx n =
1461  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1462  std::log2(d)));
1463  std::vector<idx> dims(n, d); // local dimensions vector
1464 
1465  return ptrace(rA, subsys, dims);
1466 }
1467 
1481 template<typename Derived>
1483  const Eigen::MatrixBase<Derived>& A,
1484  const std::vector<idx>& subsys,
1485  const std::vector<idx>& dims)
1486 {
1487  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1488 
1489  // EXCEPTION CHECKS
1490 
1491  // check zero-size
1493  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1494 
1495  // check that dims is a valid dimension vector
1496  if ( !internal::_check_dims(dims))
1497  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1498 
1499  // check that subsys are valid
1500  if ( !internal::_check_subsys_match_dims(subsys, dims))
1501  throw Exception("qpp::ptranspose()",
1503  // END EXCEPTION CHECKS
1504 
1505  idx D = static_cast<idx>(rA.rows());
1506  idx numdims = dims.size();
1507  idx numsubsys = subsys.size();
1508  idx Cdims[maxn];
1509  idx Cmidxcol[maxn];
1510  idx Csubsys[maxn];
1511 
1512  // copy dims in Cdims and subsys in Csubsys
1513  for ( idx i = 0; i < numdims; ++i )
1514  Cdims[i] = dims[i];
1515  for ( idx i = 0; i < numsubsys; ++i )
1516  Csubsys[i] = subsys[i];
1517 
1518  dyn_mat<typename Derived::Scalar> result(D, D);
1519 
1520  //************ ket ************//
1521  if ( internal::_check_cvector(rA)) // we have a ket
1522  {
1523  // check that dims match the dimension of A
1524  if ( !internal::_check_dims_match_cvect(dims, rA))
1525  throw Exception("qpp::ptranspose()",
1527 
1528  if ( subsys.size() == dims.size())
1529  return (rA * adjoint(rA)).transpose();
1530 
1531  if ( subsys.size() == 0 )
1532  return rA * adjoint(rA);
1533 
1534  auto worker = [=, &Cmidxcol](idx i) noexcept
1535  -> typename Derived::Scalar
1536  {
1537  // use static allocation for speed!
1538  idx midxcoltmp[maxn];
1539  idx midxrow[maxn];
1540 
1541  for ( idx k = 0; k < numdims; ++k )
1542  midxcoltmp[k] = Cmidxcol[k];
1543 
1544  /* compute the row multi-index */
1545  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1546 
1547  for ( idx k = 0; k < numsubsys; ++k )
1548  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1549 
1550  /* writes the result */
1551  return rA(internal::_multiidx2n(midxrow, numdims, Cdims)) *
1552  std::conj(rA(internal::_multiidx2n(midxcoltmp, numdims,
1553  Cdims)));
1554  }; /* end worker */
1555 
1556  for ( idx j = 0; j < D; ++j )
1557  {
1558  // compute the column multi-index
1559  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1560 
1561 #ifdef _WITH_OPENMP_
1562 #pragma omp parallel for
1563 #endif
1564  for ( idx i = 0; i < D; ++i )
1565  result(i, j) = worker(i);
1566  }
1567 
1568  return result;
1569  }
1570  //************ density matrix ************//
1571  else if ( internal::_check_square_mat(rA)) // we have a density operator
1572  {
1573  // check that dims match the dimension of A
1574  if ( !internal::_check_dims_match_mat(dims, rA))
1575  throw Exception("qpp::ptranspose()",
1577 
1578  if ( subsys.size() == dims.size())
1579  return rA.transpose();
1580 
1581  if ( subsys.size() == 0 )
1582  return rA;
1583 
1584  auto worker = [=, &Cmidxcol](idx i) noexcept
1585  -> typename Derived::Scalar
1586  {
1587  // use static allocation for speed!
1588  idx midxcoltmp[maxn];
1589  idx midxrow[maxn];
1590 
1591  for ( idx k = 0; k < numdims; ++k )
1592  midxcoltmp[k] = Cmidxcol[k];
1593 
1594  /* compute the row multi-index */
1595  internal::_n2multiidx(i, numdims, Cdims, midxrow);
1596 
1597  for ( idx k = 0; k < numsubsys; ++k )
1598  std::swap(midxcoltmp[Csubsys[k]], midxrow[Csubsys[k]]);
1599 
1600  /* writes the result */
1601  return rA(internal::_multiidx2n(midxrow, numdims, Cdims),
1602  internal::_multiidx2n(midxcoltmp, numdims, Cdims));
1603  }; /* end worker */
1604 
1605  for ( idx j = 0; j < D; ++j )
1606  {
1607  // compute the column multi-index
1608  internal::_n2multiidx(j, numdims, Cdims, Cmidxcol);
1609 
1610 #ifdef _WITH_OPENMP_
1611 #pragma omp parallel for
1612 #endif
1613  for ( idx i = 0; i < D; ++i )
1614  result(i, j) = worker(i);
1615  }
1616 
1617  return result;
1618  }
1619  //************ Exception: not ket nor density matrix ************//
1620  else
1621  throw Exception("qpp::ptranspose()",
1623 }
1624 
1638 template<typename Derived>
1640  const Eigen::MatrixBase<Derived>& A,
1641  const std::vector<idx>& subsys,
1642  idx d = 2)
1643 {
1644  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1645 
1646  // EXCEPTION CHECKS
1647 
1648  // check zero size
1650  throw Exception("qpp::ptranspose()", Exception::Type::ZERO_SIZE);
1651 
1652  // check valid dims
1653  if ( d == 0 )
1654  throw Exception("qpp::ptranspose()", Exception::Type::DIMS_INVALID);
1655  // END EXCEPTION CHECKS
1656 
1657  idx n =
1658  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1659  std::log2(d)));
1660  std::vector<idx> dims(n, d); // local dimensions vector
1661 
1662  return ptranspose(rA, subsys, dims);
1663 }
1664 
1677 template<typename Derived>
1679  const Eigen::MatrixBase<Derived>& A,
1680  const std::vector<idx>& perm,
1681  const std::vector<idx>& dims)
1682 {
1683  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1684 
1685  // EXCEPTION CHECKS
1686 
1687  // check zero-size
1689  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1690 
1691  // check that dims is a valid dimension vector
1692  if ( !internal::_check_dims(dims))
1693  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1694 
1695  // check that we have a valid permutation
1696  if ( !internal::_check_perm(perm))
1697  throw Exception("qpp::syspermute()", Exception::Type::PERM_INVALID);
1698 
1699  // check that permutation match dimensions
1700  if ( perm.size() != dims.size())
1701  throw Exception("qpp::syspermute()",
1703  // END EXCEPTION CHECKS
1704 
1705  idx D = static_cast<idx>(rA.rows());
1706  idx numdims = dims.size();
1707 
1709 
1710  //************ ket ************//
1711  if ( internal::_check_cvector(rA)) // we have a column vector
1712  {
1713  idx Cdims[maxn];
1714  idx Cperm[maxn];
1715 
1716  // check that dims match the dimension of rA
1717  if ( !internal::_check_dims_match_cvect(dims, rA))
1718  throw Exception("qpp::syspermute()",
1720 
1721  // copy dims in Cdims and perm in Cperm
1722  for ( idx i = 0; i < numdims; ++i )
1723  {
1724  Cdims[i] = dims[i];
1725  Cperm[i] = perm[i];
1726  }
1727  result.resize(D, 1);
1728 
1729  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1730  -> idx
1731  {
1732  // use static allocation for speed,
1733  // double the size for matrices reshaped as vectors
1734  idx midx[maxn];
1735  idx midxtmp[maxn];
1736  idx permdims[maxn];
1737 
1738  /* compute the multi-index */
1739  internal::_n2multiidx(i, numdims, Cdims, midx);
1740 
1741  for ( idx k = 0; k < numdims; ++k )
1742  {
1743  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1744  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1745  }
1746 
1747  return internal::_multiidx2n(midxtmp, numdims, permdims);
1748  }; /* end worker */
1749 
1750 #ifdef _WITH_OPENMP_
1751 #pragma omp parallel for
1752 #endif
1753  for ( idx i = 0; i < D; ++i )
1754  result(worker(i)) = rA(i);
1755 
1756  return result;
1757  }
1758  //************ density matrix ************//
1759  else if ( internal::_check_square_mat(rA)) // we have a density operator
1760  {
1761  idx Cdims[2 * maxn];
1762  idx Cperm[2 * maxn];
1763 
1764  // check that dims match the dimension of rA
1765  if ( !internal::_check_dims_match_mat(dims, rA))
1766  throw Exception("qpp::syspermute()",
1768 
1769  // copy dims in Cdims and perm in Cperm
1770  for ( idx i = 0; i < numdims; ++i )
1771  {
1772  Cdims[i] = dims[i];
1773  Cdims[i + numdims] = dims[i];
1774  Cperm[i] = perm[i];
1775  Cperm[i + numdims] = perm[i] + numdims;
1776  }
1777  result.resize(D * D, 1);
1778  // map A to a column vector
1780  Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1781  const_cast<typename Derived::Scalar*>(rA.data()), D * D,
1782  1);
1783 
1784  auto worker = [&Cdims, &Cperm, numdims](idx i) noexcept
1785  -> idx
1786  {
1787  // use static allocation for speed,
1788  // double the size for matrices reshaped as vectors
1789  idx midx[2 * maxn];
1790  idx midxtmp[2 * maxn];
1791  idx permdims[2 * maxn];
1792 
1793  /* compute the multi-index */
1794  internal::_n2multiidx(i, 2 * numdims, Cdims, midx);
1795 
1796  for ( idx k = 0; k < 2 * numdims; ++k )
1797  {
1798  permdims[k] = Cdims[Cperm[k]]; // permuted dimensions
1799  midxtmp[k] = midx[Cperm[k]];// permuted multi-indexes
1800  }
1801 
1802  return internal::_multiidx2n(midxtmp, 2 * numdims, permdims);
1803  }; /* end worker */
1804 
1805 #ifdef _WITH_OPENMP_
1806 #pragma omp parallel for
1807 #endif
1808  for ( idx i = 0; i < D * D; ++i )
1809  result(worker(i)) = rA(i);
1810 
1811  return reshape(result, D, D);
1812  }
1813  //************ Exception: not ket nor density matrix ************//
1814  else
1815  throw Exception("qpp::syspermute()",
1817 }
1818 
1831 template<typename Derived>
1833  const Eigen::MatrixBase<Derived>& A,
1834  const std::vector<idx>& perm,
1835  idx d = 2)
1836 {
1837  const typename Eigen::MatrixBase<Derived>::EvalReturnType& rA = A.derived();
1838 
1839  // EXCEPTION CHECKS
1840 
1841  // check zero size
1843  throw Exception("qpp::syspermute()", Exception::Type::ZERO_SIZE);
1844 
1845  // check valid dims
1846  if ( d == 0 )
1847  throw Exception("qpp::syspermute()", Exception::Type::DIMS_INVALID);
1848  // END EXCEPTION CHECKS
1849 
1850  idx n =
1851  static_cast<idx>(std::llround(std::log2(rA.rows()) /
1852  std::log2(d)));
1853  std::vector<idx> dims(n, d); // local dimensions vector
1854 
1855  return syspermute(rA, perm, dims);
1856 }
1857 
1858 } /* namespace qpp */
1859 
1860 #endif /* OPERATIONS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:74
std::vector< cmat > choi2kraus(const cmat &A)
Orthogonal Kraus operators from Choi matrix.
Definition: operations.h:910
bool _check_subsys_match_dims(const std::vector< idx > &subsys, const std::vector< idx > &dims)
Definition: util.h:183
bool _check_dims_match_mat(const std::vector< idx > &dims, const Eigen::MatrixBase< Derived > &A)
Definition: util.h:141
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:152
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:67
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h: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:1678
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:492
Quantum++ main namespace.
Definition: codes.h:30
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:84
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:946
cmat super2choi(const cmat &A)
Converts superoperator matrix to Choi matrix.
Definition: operations.h:982
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1818
dyn_mat< typename Derived::Scalar > ptrace1(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1026
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:900
dyn_mat< typename Derived::Scalar > ptrace2(const Eigen::MatrixBase< Derived > &A, const std::vector< idx > &dims)
Partial trace.
Definition: operations.h:1128
bool _check_perm(const std::vector< idx > &perm)
Definition: util.h:242
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:61
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:112
cmat kraus2choi(const std::vector< cmat > &Ks)
Choi matrix.
Definition: operations.h:851
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:1220
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:773
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:1482
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:48
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:125
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:1140