Quantum++  v1.0
A modern C++11 quantum computing library
functions.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 FUNCTIONS_H_
33 #define FUNCTIONS_H_
34 
35 namespace qpp {
36 // Eigen function wrappers
44 template <typename Derived>
45 dyn_mat<typename Derived::Scalar>
46 transpose(const Eigen::MatrixBase<Derived>& A) {
47  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
48 
49  // check zero-size
51  throw exception::ZeroSize("qpp::transpose()");
52 
53  return rA.transpose();
54 }
55 
63 template <typename Derived>
65 conjugate(const Eigen::MatrixBase<Derived>& A) {
66  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
67 
68  // check zero-size
70  throw exception::ZeroSize("qpp::conjugate()");
71 
72  return rA.conjugate();
73 }
74 
82 template <typename Derived>
83 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A) {
84  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
85 
86  // EXCEPTION CHECKS
87 
88  // check zero-size
90  throw exception::ZeroSize("qpp::adjoint()");
91  // END EXCEPTION CHECKS
92 
93  return rA.adjoint();
94 }
95 
103 template <typename Derived>
104 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A) {
105  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
106 
107  // EXCEPTION CHECKS
108 
109  // check zero-size
111  throw exception::ZeroSize("qpp::inverse()");
112  // END EXCEPTION CHECKS
113 
114  return rA.inverse();
115 }
116 
123 template <typename Derived>
124 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A) {
125  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
126 
127  // EXCEPTION CHECKS
128 
129  // check zero-size
131  throw exception::ZeroSize("qpp::trace()");
132  // END EXCEPTION CHECKS
133 
134  return rA.trace();
135 }
136 
144 template <typename Derived>
145 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A) {
146  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
147 
148  // EXCEPTION CHECKS
149 
150  // check zero-size
152  throw exception::ZeroSize("qpp::det()");
153  // END EXCEPTION CHECKS
154 
155  return rA.determinant();
156 }
157 
167 template <typename Derived>
168 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A) {
169  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
170 
171  // EXCEPTION CHECKS
172 
173  // check zero-size
175  throw exception::ZeroSize("qpp::logdet()");
176 
177  // check square matrix
179  throw exception::MatrixNotSquare("qpp::logdet()");
180  // END EXCEPTION CHECKS
181 
182  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
184  lu.matrixLU().template triangularView<Eigen::Upper>();
185  typename Derived::Scalar result = std::log(U(0, 0));
186 
187  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
188  result += std::log(U(i, i));
189 
190  return result;
191 }
192 
200 template <typename Derived>
201 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A) {
202  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
203 
204  // EXCEPTION CHECKS
205 
206  // check zero-size
208  throw exception::ZeroSize("qpp::sum()");
209  // END EXCEPTION CHECKS
210 
211  return rA.sum();
212 }
213 
221 template <typename Derived>
222 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A) {
223  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
224 
225  // EXCEPTION CHECKS
226 
227  // check zero-size
229  throw exception::ZeroSize("qpp::prod()");
230  // END EXCEPTION CHECKS
231 
232  return rA.prod();
233 }
234 
241 template <typename Derived>
242 double norm(const Eigen::MatrixBase<Derived>& A) {
243  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
244 
245  // EXCEPTION CHECKS
246 
247  // check zero-size
249  throw exception::ZeroSize("qpp::norm()");
250  // END EXCEPTION CHECKS
251 
252  // convert matrix to complex then return its norm
253  return (rA.template cast<cplx>()).norm();
254 }
255 
264 template <typename Derived>
265 std::pair<dyn_col_vect<cplx>, cmat>
266 
267 eig(const Eigen::MatrixBase<Derived>& A) {
268  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
269 
270  // EXCEPTION CHECKS
271 
272  // check zero-size
274  throw exception::ZeroSize("qpp::eig()");
275 
276  // check square matrix
278  throw exception::MatrixNotSquare("qpp::eig()");
279  // END EXCEPTION CHECKS
280 
281  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
282 
283  return std::make_pair(es.eigenvalues(), es.eigenvectors());
284 }
285 
293 template <typename Derived>
294 dyn_col_vect<cplx> evals(const Eigen::MatrixBase<Derived>& A) {
295  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
296 
297  // EXCEPTION CHECKS
298 
299  // check zero-size
301  throw exception::ZeroSize("qpp::evals()");
302 
303  // check square matrix
305  throw exception::MatrixNotSquare("qpp::evals()");
306  // END EXCEPTION CHECKS
307 
308  return eig(rA).first;
309 }
310 
318 template <typename Derived>
319 cmat evects(const Eigen::MatrixBase<Derived>& A) {
320  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
321 
322  // EXCEPTION CHECKS
323 
324  // check zero-size
326  throw exception::ZeroSize("qpp::evects()");
327 
328  // check square matrix
330  throw exception::MatrixNotSquare("qpp::evects()");
331  // END EXCEPTION CHECKS
332 
333  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
334 
335  return eig(rA).second;
336 }
337 
346 template <typename Derived>
347 std::pair<dyn_col_vect<double>, cmat>
348 
349 heig(const Eigen::MatrixBase<Derived>& A) {
350  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
351 
352  // EXCEPTION CHECKS
353 
354  // check zero-size
356  throw exception::ZeroSize("qpp::heig()");
357 
358  // check square matrix
360  throw exception::MatrixNotSquare("qpp::heig()");
361  // END EXCEPTION CHECKS
362 
363  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
364 
365  return std::make_pair(es.eigenvalues(), es.eigenvectors());
366 }
367 
375 template <typename Derived>
376 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A) {
377  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
378 
379  // EXCEPTION CHECKS
380 
381  // check zero-size
383  throw exception::ZeroSize("qpp::hevals()");
384 
385  // check square matrix
387  throw exception::MatrixNotSquare("qpp::hevals()");
388  // END EXCEPTION CHECKS
389 
390  return heig(rA).first;
391 }
392 
400 template <typename Derived>
401 cmat hevects(const Eigen::MatrixBase<Derived>& A) {
402  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
403 
404  // EXCEPTION CHECKS
405 
406  // check zero-size
408  throw exception::ZeroSize("qpp::hevects()");
409 
410  // check square matrix
412  throw exception::MatrixNotSquare("qpp::hevects()");
413  // END EXCEPTION CHECKS
414 
415  return heig(rA).second;
416 }
417 
427 template <typename Derived>
428 std::tuple<cmat, dyn_col_vect<double>, cmat>
429 
430 svd(const Eigen::MatrixBase<Derived>& A) {
431  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
432 
433  // EXCEPTION CHECKS
434 
435  // check zero-size
437  throw exception::ZeroSize("qpp::svd()");
438  // END EXCEPTION CHECKS
439 
440  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
441  rA, Eigen::DecompositionOptions::ComputeFullU |
442  Eigen::DecompositionOptions::ComputeFullV);
443 
444  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
445 }
446 
454 template <typename Derived>
455 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A) {
456  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
457 
458  // EXCEPTION CHECKS
459 
460  // check zero-size
462  throw exception::ZeroSize("qpp::svals()");
463  // END EXCEPTION CHECKS
464 
465  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
466 
467  return sv.singularValues();
468 }
469 
477 template <typename Derived>
478 cmat svdU(const Eigen::MatrixBase<Derived>& A) {
479  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
480 
481  // EXCEPTION CHECKS
482 
483  // check zero-size
485  throw exception::ZeroSize("qpp::svdU()");
486  // END EXCEPTION CHECKS
487 
488  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
489  rA, Eigen::DecompositionOptions::ComputeFullU);
490 
491  return sv.matrixU();
492 }
493 
501 template <typename Derived>
502 cmat svdV(const Eigen::MatrixBase<Derived>& A) {
503  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
504 
505  // EXCEPTION CHECKS
506 
507  // check zero-size
509  throw exception::ZeroSize("qpp::svdV()");
510  // END EXCEPTION CHECKS
511 
512  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
513  rA, Eigen::DecompositionOptions::ComputeFullV);
514 
515  return sv.matrixV();
516 }
517 
518 // Matrix functional calculus
519 
527 template <typename Derived>
528 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (*f)(const cplx&)) {
529  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
530 
531  // EXCEPTION CHECKS
532 
533  // check zero-size
535  throw exception::ZeroSize("qpp::funm()");
536 
537  // check square matrix
539  throw exception::MatrixNotSquare("qpp::funm()");
540  // END EXCEPTION CHECKS
541 
542  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
543  cmat evects = es.eigenvectors();
544  cmat evals = es.eigenvalues();
545  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
546  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
547 
548  cmat evalsdiag = evals.asDiagonal();
549 
550  return evects * evalsdiag * evects.inverse();
551 }
552 
559 template <typename Derived>
560 cmat sqrtm(const Eigen::MatrixBase<Derived>& A) {
561  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
562 
563  // EXCEPTION CHECKS
564 
565  // check zero-size
567  throw exception::ZeroSize("qpp::sqrtm()");
568 
569  // check square matrix
571  throw exception::MatrixNotSquare("qpp::sqrtm()");
572  // END EXCEPTION CHECKS
573 
574  return funm(rA, &std::sqrt);
575 }
576 
583 template <typename Derived>
584 cmat absm(const Eigen::MatrixBase<Derived>& A) {
585  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
586 
587  // EXCEPTION CHECKS
588 
589  // check zero-size
591  throw exception::ZeroSize("qpp::absm()");
592 
593  // check square matrix
595  throw exception::MatrixNotSquare("qpp::absm()");
596  // END EXCEPTION CHECKS
597 
598  return sqrtm(adjoint(rA) * rA);
599 }
600 
607 template <typename Derived>
608 cmat expm(const Eigen::MatrixBase<Derived>& A) {
609  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
610 
611  // EXCEPTION CHECKS
612 
613  // check zero-size
615  throw exception::ZeroSize("qpp::expm()");
616 
617  // check square matrix
619  throw exception::MatrixNotSquare("qpp::expm()");
620  // END EXCEPTION CHECKS
621 
622  return funm(rA, &std::exp);
623 }
624 
631 template <typename Derived>
632 cmat logm(const Eigen::MatrixBase<Derived>& A) {
633  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
634 
635  // EXCEPTION CHECKS
636 
637  // check zero-size
639  throw exception::ZeroSize("qpp::logm()");
640 
641  // check square matrix
643  throw exception::MatrixNotSquare("qpp::logm()");
644  // END EXCEPTION CHECKS
645 
646  return funm(rA, &std::log);
647 }
648 
655 template <typename Derived>
656 cmat sinm(const Eigen::MatrixBase<Derived>& A) {
657  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
658 
659  // EXCEPTION CHECKS
660 
661  // check zero-size
663  throw exception::ZeroSize("qpp::sinm()");
664 
665  // check square matrix
667  throw exception::MatrixNotSquare("qpp::sinm()");
668  // END EXCEPTION CHECKS
669 
670  return funm(rA, &std::sin);
671 }
672 
679 template <typename Derived>
680 cmat cosm(const Eigen::MatrixBase<Derived>& A) {
681  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
682 
683  // EXCEPTION CHECKS
684 
685  // check zero-size
687  throw exception::ZeroSize("qpp::cosm()");
688 
689  // check square matrix
691  throw exception::MatrixNotSquare("qpp::cosm()");
692  // END EXCEPTION CHECKS
693 
694  return funm(rA, &std::cos);
695 }
696 
708 template <typename Derived>
709 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z) {
710  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
711 
712  // EXCEPTION CHECKS
713 
714  // check zero-size
716  throw exception::ZeroSize("qpp::spectralpowm()");
717 
718  // check square matrix
720  throw exception::MatrixNotSquare("qpp::spectralpowm()");
721  // END EXCEPTION CHECKS
722 
723  // Define A^0 = Id, for z IDENTICALLY zero
724  if (real(z) == 0 && imag(z) == 0)
725  return cmat::Identity(rA.rows(), rA.rows());
726 
727  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
728  cmat evects = es.eigenvectors();
729  cmat evals = es.eigenvalues();
730  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
731  evals(i) = std::pow(evals(i), z);
732 
733  cmat evalsdiag = evals.asDiagonal();
734 
735  return evects * evalsdiag * evects.inverse();
736 }
737 
750 template <typename Derived>
751 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
752  idx n) {
753  // EXCEPTION CHECKS
754 
755  // check zero-size
757  throw exception::ZeroSize("qpp::powm()");
758 
759  // check square matrix
761  throw exception::MatrixNotSquare("qpp::powm()");
762  // END EXCEPTION CHECKS
763 
764  // if n = 1, return the matrix unchanged
765  if (n == 1)
766  return A;
767 
770 
771  // if n = 0, return the identity (as just prepared in result)
772  if (n == 0)
773  return result;
774 
775  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
776 
777  // fast matrix power
778  for (; n > 0; n /= 2) {
779  if (n % 2)
780  result = (result * cA).eval();
781  cA = (cA * cA).eval();
782  }
783 
784  return result;
785 }
786 
795 template <typename Derived>
796 double schatten(const Eigen::MatrixBase<Derived>& A, double p) {
797  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
798 
799  // EXCEPTION CHECKS
800 
801  // check zero-size
803  throw exception::ZeroSize("qpp::schatten()");
804  if (p < 1)
805  throw exception::OutOfRange("qpp::schatten()");
806  // END EXCEPTION CHECKS
807 
808  if (p == infty) // infinity norm (largest singular value)
809  return svals(rA)(0);
810 
811  const dyn_col_vect<double> sv = svals(rA);
812  double result = 0;
813  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
814  result += std::pow(sv[i], p);
815 
816  return std::pow(result, 1. / p);
817 }
818 
819 // other functions
820 
829 template <typename OutputScalar, typename Derived>
831 cwise(const Eigen::MatrixBase<Derived>& A,
832  OutputScalar (*f)(const typename Derived::Scalar&)) {
833  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
834 
835  // EXCEPTION CHECKS
836 
837  // check zero-size
839  throw exception::ZeroSize("qpp::cwise()");
840  // END EXCEPTION CHECKS
841 
842  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
843 
844 #ifdef WITH_OPENMP_
845 #pragma omp parallel for collapse(2)
846 #endif // WITH_OPENMP_
847  // column major order for speed
848  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
849  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
850  result(i, j) = (*f)(rA(i, j));
851 
852  return result;
853 }
854 
855 // Kronecker product of multiple matrices, preserve return type
856 // variadic template
867 template <typename T>
869  return head;
870 }
871 
882 template <typename T, typename... Args>
883 dyn_mat<typename T::Scalar> kron(const T& head, const Args&... tail) {
884  return internal::kron2(head, kron(tail...));
885 }
886 
896 template <typename Derived>
897 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As) {
898  // EXCEPTION CHECKS
899 
900  if (As.size() == 0)
901  throw exception::ZeroSize("qpp::kron()");
902 
903  for (auto&& it : As)
905  throw exception::ZeroSize("qpp::kron()");
906  // END EXCEPTION CHECKS
907 
908  dyn_mat<typename Derived::Scalar> result = As[0].derived();
909  for (idx i = 1; i < As.size(); ++i) {
910  result = kron(result, As[i]);
911  }
912 
913  return result;
914 }
915 
916 // Kronecker product of a list of matrices, preserve return type
917 // deduce the template parameters from initializer_list
928 template <typename Derived>
930 kron(const std::initializer_list<Derived>& As) {
931  return kron(std::vector<Derived>(As));
932 }
933 
943 template <typename Derived>
944 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
945  idx n) {
946  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
947 
948  // EXCEPTION CHECKS
949 
950  // check zero-size
952  throw exception::ZeroSize("qpp::kronpow()");
953 
954  // check out of range
955  if (n == 0)
956  throw exception::OutOfRange("qpp::kronpow()");
957  // END EXCEPTION CHECKS
958 
959  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
960 
961  return kron(As);
962 }
963 
964 // Direct sum of multiple matrices, preserve return type
965 // variadic template
976 template <typename T>
978  return head;
979 }
980 
991 template <typename T, typename... Args>
992 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args&... tail) {
993  return internal::dirsum2(head, dirsum(tail...));
994 }
995 
1005 template <typename Derived>
1006 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As) {
1007  // EXCEPTION CHECKS
1008 
1009  if (As.size() == 0)
1010  throw exception::ZeroSize("qpp::dirsum()");
1011 
1012  for (auto&& it : As)
1014  throw exception::ZeroSize("qpp::dirsum()");
1015  // END EXCEPTION CHECKS
1016 
1017  idx total_rows = 0, total_cols = 0;
1018  for (idx i = 0; i < As.size(); ++i) {
1019  total_rows += static_cast<idx>(As[i].rows());
1020  total_cols += static_cast<idx>(As[i].cols());
1021  }
1023  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1024 
1025  idx cur_row = 0, cur_col = 0;
1026  for (idx i = 0; i < As.size(); ++i) {
1027  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1028  cur_row += static_cast<idx>(As[i].rows());
1029  cur_col += static_cast<idx>(As[i].cols());
1030  }
1031 
1032  return result;
1033 }
1034 
1035 // Direct sum of a list of matrices, preserve return type
1036 // deduce the template parameters from initializer_list
1047 template <typename Derived>
1049 dirsum(const std::initializer_list<Derived>& As) {
1050  return dirsum(std::vector<Derived>(As));
1051 }
1052 
1062 template <typename Derived>
1063 dyn_mat<typename Derived::Scalar> dirsumpow(const Eigen::MatrixBase<Derived>& A,
1064  idx n) {
1065  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1066 
1067  // EXCEPTION CHECKS
1068 
1069  // check zero-size
1071  throw exception::ZeroSize("qpp::dirsumpow()");
1072 
1073  // check out of range
1074  if (n == 0)
1075  throw exception::OutOfRange("qpp::dirsumpow()");
1076  // END EXCEPTION CHECKS
1077 
1078  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1079 
1080  return dirsum(As);
1081 }
1082 
1094 template <typename Derived>
1095 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1096  idx rows, idx cols) {
1097  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1098 
1099  idx Arows = static_cast<idx>(rA.rows());
1100  idx Acols = static_cast<idx>(rA.cols());
1101 
1102  // EXCEPTION CHECKS
1103 
1104  // check zero-size
1106  throw exception::ZeroSize("qpp::reshape()");
1107 
1108  if (Arows * Acols != rows * cols)
1109  throw exception::DimsMismatchMatrix("qpp::reshape()");
1110  // END EXCEPTION CHECKS
1111 
1112  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1113  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1114 }
1115 
1128 template <typename Derived1, typename Derived2>
1129 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1130  const Eigen::MatrixBase<Derived2>& B) {
1131  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1132  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1133 
1134  // EXCEPTION CHECKS
1135 
1136  // check types
1137  if (!std::is_same<typename Derived1::Scalar,
1138  typename Derived2::Scalar>::value)
1139  throw exception::TypeMismatch("qpp::comm()");
1140 
1141  // check zero-size
1143  throw exception::ZeroSize("qpp::comm()");
1144 
1145  // check square matrices
1147  throw exception::MatrixNotSquare("qpp::comm()");
1148 
1149  // check equal dimensions
1150  if (rA.rows() != rB.rows())
1151  throw exception::DimsNotEqual("qpp::comm()");
1152  // END EXCEPTION CHECKS
1153 
1154  return rA * rB - rB * rA;
1155 }
1156 
1169 template <typename Derived1, typename Derived2>
1171 anticomm(const Eigen::MatrixBase<Derived1>& A,
1172  const Eigen::MatrixBase<Derived2>& B) {
1173  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1174  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1175 
1176  // EXCEPTION CHECKS
1177 
1178  // check types
1179  if (!std::is_same<typename Derived1::Scalar,
1180  typename Derived2::Scalar>::value)
1181  throw exception::TypeMismatch("qpp::anticomm()");
1182 
1183  // check zero-size
1185  throw exception::ZeroSize("qpp::anticomm()");
1186 
1187  // check square matrices
1189  throw exception::MatrixNotSquare("qpp::anticomm()");
1190 
1191  // check equal dimensions
1192  if (rA.rows() != rB.rows())
1193  throw exception::DimsNotEqual("qpp::anticomm()");
1194  // END EXCEPTION CHECKS
1195 
1196  return rA * rB + rB * rA;
1197 }
1198 
1209 template <typename Derived>
1210 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A) {
1211  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1212 
1213  // EXCEPTION CHECKS
1214 
1215  // check zero-size
1217  throw exception::ZeroSize("qpp::prj()");
1218 
1219  // check column vector
1220  if (!internal::check_cvector(rA))
1221  throw exception::MatrixNotCvector("qpp::prj()");
1222  // END EXCEPTION CHECKS
1223 
1224  double normA = norm(rA);
1225  if (normA > eps)
1226  return rA * adjoint(rA) / (normA * normA);
1227  else
1228  return dyn_mat<typename Derived::Scalar>::Zero(rA.rows(), rA.rows());
1229 }
1230 
1238 template <typename Derived>
1239 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As) {
1240  // EXCEPTION CHECKS
1241 
1242  // check empty list
1244  throw exception::ZeroSize("qpp::grams()");
1245 
1246  for (auto&& it : As)
1248  throw exception::ZeroSize("qpp::grams()");
1249 
1250  // check that As[0] is a column vector
1251  if (!internal::check_cvector(As[0]))
1252  throw exception::MatrixNotCvector("qpp::grams()");
1253 
1254  // now check that all the rest match the size of the first vector
1255  for (auto&& it : As)
1256  if (it.rows() != As[0].rows() || it.cols() != 1)
1257  throw exception::DimsNotEqual("qpp::grams()");
1258  // END EXCEPTION CHECKS
1259 
1261  dyn_mat<typename Derived::Scalar>::Identity(As[0].rows(), As[0].rows());
1262 
1264  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1265 
1266  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1267  // find the first non-zero vector in the list
1268  idx pos = 0;
1269  for (pos = 0; pos < As.size(); ++pos) {
1270  if (norm(As[pos]) > eps) // add it as the first element
1271  {
1272  outvecs.push_back(As[pos]);
1273  break;
1274  }
1275  }
1276 
1277  // start the process
1278  for (idx i = pos + 1; i < As.size(); ++i) {
1279  cut -= prj(outvecs[i - 1 - pos]);
1280  vi = cut * As[i];
1281  outvecs.push_back(vi);
1282  }
1283 
1284  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1285 
1286  idx cnt = 0;
1287  for (auto&& it : outvecs) {
1288  double normA = norm(it);
1289  if (normA > eps) // we add only the non-zero vectors
1290  {
1291  result.col(cnt) = it / normA;
1292  cnt++;
1293  }
1294  }
1295 
1296  return result.block(0, 0, As[0].rows(), cnt);
1297 }
1298 
1299 // deduce the template parameters from initializer_list
1307 template <typename Derived>
1309 grams(const std::initializer_list<Derived>& As) {
1310  return grams(std::vector<Derived>(As));
1311 }
1312 
1320 template <typename Derived>
1321 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A) {
1322  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1323 
1324  // EXCEPTION CHECKS
1325 
1327  throw exception::ZeroSize("qpp::grams()");
1328  // END EXCEPTION CHECKS
1329 
1330  std::vector<dyn_mat<typename Derived::Scalar>> input;
1331 
1332  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1333  input.push_back(rA.col(i));
1334 
1335  return grams<dyn_mat<typename Derived::Scalar>>(input);
1336 }
1337 
1348 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims) {
1349  // EXCEPTION CHECKS
1350 
1351  if (!internal::check_dims(dims))
1352  throw exception::DimsInvalid("qpp::n2multiidx()");
1353 
1354  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1355  static_cast<idx>(1), std::multiplies<idx>()))
1356  throw exception::OutOfRange("qpp::n2multiidx()");
1357  // END EXCEPTION CHECKS
1358 
1359  // double the size for matrices reshaped as vectors
1360  idx result[2 * maxn];
1361  internal::n2multiidx(n, dims.size(), dims.data(), result);
1362 
1363  return std::vector<idx>(result, result + dims.size());
1364 }
1365 
1376 inline idx multiidx2n(const std::vector<idx>& midx,
1377  const std::vector<idx>& dims) {
1378  // EXCEPTION CHECKS
1379 
1380  if (!internal::check_dims(dims))
1381  throw exception::DimsInvalid("qpp::multiidx2n()");
1382 
1383  for (idx i = 0; i < dims.size(); ++i)
1384  if (midx[i] >= dims[i])
1385  throw exception::OutOfRange("qpp::multiidx2n()");
1386  // END EXCEPTION CHECKS
1387 
1388  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1389 }
1390 
1405 inline ket mket(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1406  idx N = mask.size();
1407 
1408  idx D = std::accumulate(std::begin(dims), std::end(dims),
1409  static_cast<idx>(1), std::multiplies<idx>());
1410 
1411  // EXCEPTION CHECKS
1412 
1413  // check zero size
1414  if (N == 0)
1415  throw exception::ZeroSize("qpp::mket()");
1416  // check valid dims
1417  if (!internal::check_dims(dims))
1418  throw exception::DimsInvalid("qpp::mket()");
1419  // check mask and dims have the same size
1420  if (mask.size() != dims.size())
1421  throw exception::SubsysMismatchDims("qpp::mket()");
1422  // check mask is a valid vector
1423  for (idx i = 0; i < N; ++i)
1424  if (mask[i] >= dims[i])
1425  throw exception::SubsysMismatchDims("qpp::mket()");
1426  // END EXCEPTION CHECKS
1427 
1428  ket result = ket::Zero(D);
1429  idx pos = multiidx2n(mask, dims);
1430  result(pos) = 1;
1431 
1432  return result;
1433 }
1434 
1448 inline ket mket(const std::vector<idx>& mask, idx d = 2) {
1449  idx N = mask.size();
1450  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1451 
1452  // EXCEPTION CHECKS
1453 
1454  // check zero size
1455  if (N == 0)
1456  throw exception::ZeroSize("qpp::mket()");
1457 
1458  // check valid dims
1459  if (d == 0)
1460  throw exception::DimsInvalid("qpp::mket()");
1461 
1462  // check mask is a valid vector
1463  for (idx i = 0; i < N; ++i)
1464  if (mask[i] >= d)
1465  throw exception::SubsysMismatchDims("qpp::mket()");
1466  // END EXCEPTION CHECKS
1467 
1468  ket result = ket::Zero(D);
1469  std::vector<idx> dims(N, d);
1470  idx pos = multiidx2n(mask, dims);
1471  result(pos) = 1;
1472 
1473  return result;
1474 }
1475 
1491 inline cmat mprj(const std::vector<idx>& mask, const std::vector<idx>& dims) {
1492  idx N = mask.size();
1493 
1494  idx D = std::accumulate(std::begin(dims), std::end(dims),
1495  static_cast<idx>(1), std::multiplies<idx>());
1496 
1497  // EXCEPTION CHECKS
1498 
1499  // check zero size
1500  if (N == 0)
1501  throw exception::ZeroSize("qpp::mprj()");
1502  // check valid dims
1503  if (!internal::check_dims(dims))
1504  throw exception::DimsInvalid("qpp::mprj()");
1505  // check mask and dims have the same size
1506  if (mask.size() != dims.size())
1507  throw exception::SubsysMismatchDims("qpp::mprj()");
1508  // check mask is a valid vector
1509  for (idx i = 0; i < N; ++i)
1510  if (mask[i] >= dims[i])
1511  throw exception::SubsysMismatchDims("qpp::mprj()");
1512  // END EXCEPTION CHECKS
1513 
1514  cmat result = cmat::Zero(D, D);
1515  idx pos = multiidx2n(mask, dims);
1516  result(pos, pos) = 1;
1517 
1518  return result;
1519 }
1520 
1536 inline cmat mprj(const std::vector<idx>& mask, idx d = 2) {
1537  idx N = mask.size();
1538  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1539 
1540  // EXCEPTION CHECKS
1541 
1542  // check zero size
1543  if (N == 0)
1544  throw exception::ZeroSize("qpp::mprj()");
1545 
1546  // check valid dims
1547  if (d == 0)
1548  throw exception::DimsInvalid("qpp::mprj()");
1549 
1550  // check mask is a valid vector
1551  for (idx i = 0; i < N; ++i)
1552  if (mask[i] >= d)
1553  throw exception::SubsysMismatchDims("qpp::mprj()");
1554  // END EXCEPTION CHECKS
1555 
1556  cmat result = cmat::Zero(D, D);
1557  std::vector<idx> dims(N, d);
1558  idx pos = multiidx2n(mask, dims);
1559  result(pos, pos) = 1;
1560 
1561  return result;
1562 }
1563 
1572 template <typename InputIterator>
1573 std::vector<double> abssq(InputIterator first, InputIterator last) {
1574  std::vector<double> weights(std::distance(first, last));
1575  std::transform(first, last, std::begin(weights),
1576  [](cplx z) -> double { return std::norm(z); });
1577 
1578  return weights;
1579 }
1580 
1587 template <typename Container>
1588 std::vector<double>
1589 abssq(const Container& c,
1590  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1591 // we need the std::enable_if to SFINAE out Eigen expressions
1592 // that will otherwise match, instead of matching
1593 // the overload below:
1594 
1595 // template<typename Derived>
1596 // abssq(const Eigen::MatrixBase<Derived>& A)
1597 {
1598  return abssq(std::begin(c), std::end(c));
1599 }
1600 
1607 template <typename Derived>
1608 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A) {
1609  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1610 
1611  // EXCEPTION CHECKS
1612 
1613  // check zero-size
1615  throw exception::ZeroSize("qpp::abssq()");
1616  // END EXCEPTION CHECKS
1617 
1618  return abssq(rA.data(), rA.data() + rA.size());
1619 }
1620 
1629 template <typename InputIterator>
1630 typename std::iterator_traits<InputIterator>::value_type
1631 sum(InputIterator first, InputIterator last) {
1632  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1633 
1634  return std::accumulate(first, last, static_cast<value_type>(0));
1635 }
1636 
1644 template <typename Container>
1645 typename Container::value_type
1646 sum(const Container& c,
1647  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1648  return sum(std::begin(c), std::end(c));
1649 }
1650 
1659 template <typename InputIterator>
1660 typename std::iterator_traits<InputIterator>::value_type
1661 prod(InputIterator first, InputIterator last) {
1662  using value_type = typename std::iterator_traits<InputIterator>::value_type;
1663 
1664  return std::accumulate(first, last, static_cast<value_type>(1),
1665  std::multiplies<value_type>());
1666 }
1667 
1675 template <typename Container>
1676 typename Container::value_type
1677 prod(const Container& c,
1678  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr) {
1679  return prod(std::begin(c), std::end(c));
1680 }
1681 
1694 template <typename Derived>
1696 rho2pure(const Eigen::MatrixBase<Derived>& A) {
1697  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1698 
1699  // EXCEPTION CHECKS
1700  // check zero-size
1702  throw exception::ZeroSize("qpp::rho2pure()");
1703  // check square matrix
1704  if (!internal::check_square_mat(rA))
1705  throw exception::MatrixNotSquare("qpp::rho2pure()");
1706  // END EXCEPTION CHECKS
1707 
1708  dyn_col_vect<double> tmp_evals = hevals(rA);
1709  cmat tmp_evects = hevects(rA);
1712  // find the non-zero eigenvector
1713  // there is only one, assuming the state is pure
1714  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k) {
1715  if (std::abs(tmp_evals(k)) > eps) {
1716  result = tmp_evects.col(k);
1717  break;
1718  }
1719  }
1720 
1721  return result;
1722 }
1723 
1732 template <typename T>
1733 std::vector<T> complement(std::vector<T> subsys, idx N) {
1734  // EXCEPTION CHECKS
1735 
1736  if (N < subsys.size())
1737  throw exception::OutOfRange("qpp::complement()");
1738  // END EXCEPTION CHECKS
1739 
1740  std::vector<T> all(N);
1741  std::vector<T> subsys_bar(N - subsys.size());
1742 
1743  std::iota(std::begin(all), std::end(all), 0);
1744  std::sort(std::begin(subsys), std::end(subsys));
1745  std::set_difference(std::begin(all), std::end(all), std::begin(subsys),
1746  std::end(subsys), std::begin(subsys_bar));
1747 
1748  return subsys_bar;
1749 }
1750 
1761 template <typename Derived>
1762 std::vector<double> rho2bloch(const Eigen::MatrixBase<Derived>& A) {
1763  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1764 
1765  // EXCEPTION CHECKS
1766 
1767  // check qubit matrix
1769  throw exception::NotQubitMatrix("qpp::rho2bloch()");
1770  // END EXCEPTION CHECKS
1771 
1772  std::vector<double> result(3);
1773  cmat X(2, 2), Y(2, 2), Z(2, 2);
1774  X << 0, 1, 1, 0;
1775  Y << 0, -1_i, 1_i, 0;
1776  Z << 1, 0, 0, -1;
1777  result[0] = std::real(trace(rA * X));
1778  result[1] = std::real(trace(rA * Y));
1779  result[2] = std::real(trace(rA * Z));
1780 
1781  return result;
1782 }
1783 
1792 inline cmat bloch2rho(const std::vector<double>& r) {
1793  // EXCEPTION CHECKS
1794 
1795  // check 3-dimensional vector
1796  if (r.size() != 3)
1797  throw exception::CustomException("qpp::bloch2rho",
1798  "r is not a 3-dimensional vector!");
1799  // END EXCEPTION CHECKS
1800 
1801  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1802  X << 0, 1, 1, 0;
1803  Y << 0, -1_i, 1_i, 0;
1804  Z << 1, 0, 0, -1;
1805  Id2 << 1, 0, 0, 1;
1806 
1807  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1808 }
1809 
1810 inline namespace literals {
1811 // Idea taken from http://techblog.altplus.co.jp/entry/2017/11/08/130921
1821 template <char... Bits>
1822 ket operator"" _ket() {
1823  constexpr idx n = sizeof...(Bits);
1824  constexpr char bits[n + 1] = {Bits..., '\0'};
1825  qpp::ket q = qpp::ket::Zero(std::pow(2, n));
1826  idx pos = 0;
1827 
1828  // EXCEPTION CHECKS
1829 
1830  // check valid multi-partite qubit state
1831  for (idx i = 0; i < n; ++i) {
1832  if (bits[i] != '0' && bits[i] != '1')
1833  throw exception::OutOfRange(R"xxx(qpp::operator "" _ket())xxx");
1834  }
1835  // END EXCEPTION CHECKS
1836 
1837  pos = std::stoi(bits, nullptr, 2);
1838  q(pos) = 1;
1839 
1840  return q;
1841 }
1842 
1852 template <char... Bits>
1853 bra operator"" _bra() {
1854  constexpr idx n = sizeof...(Bits);
1855  constexpr char bits[n + 1] = {Bits..., '\0'};
1856  qpp::bra q = qpp::ket::Zero(std::pow(2, n));
1857  idx pos = 0;
1858 
1859  // EXCEPTION CHECKS
1860 
1861  // check valid multi-partite qubit state
1862  for (idx i = 0; i < n; ++i) {
1863  if (bits[i] != '0' && bits[i] != '1')
1864  throw exception::OutOfRange(R"xxx(qpp::operator "" _bra())xxx");
1865  }
1866  // END EXCEPTION CHECKS
1867 
1868  pos = std::stoi(bits, nullptr, 2);
1869  q(pos) = 1;
1870 
1871  return q;
1872 }
1873 
1885 template <char... Bits>
1886 cmat operator"" _prj() {
1887  constexpr idx n = sizeof...(Bits);
1888  constexpr char bits[n + 1] = {Bits..., '\0'};
1889 
1890  // EXCEPTION CHECKS
1891 
1892  // check valid multi-partite qubit state
1893  for (idx i = 0; i < n; ++i) {
1894  if (bits[i] != '0' && bits[i] != '1')
1895  throw exception::OutOfRange(R"xxx(qpp::operator "" _prj())xxx");
1896  }
1897  // END EXCEPTION CHECKS
1898 
1899  return kron(operator""_ket<Bits...>(), operator""_bra<Bits...>());
1900 }
1901 } /* inline namespace literals */
1902 
1903 } /* namespace qpp */
1904 
1905 #endif /* FUNCTIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:349
Dimensions not equal exception.
Definition: exception.h:284
Dimension(s) mismatch matrix size exception.
Definition: exception.h:300
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:168
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:75
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolute value.
Definition: functions.h:584
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:273
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
Custom exception.
Definition: exception.h:571
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:608
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:222
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:68
Eigen::RowVectorXcd bra
Complex (double precision) dynamic Eigen row vector.
Definition: types.h:59
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:977
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:81
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
std::vector< double > rho2bloch(const Eigen::MatrixBase< Derived > &A)
Computes the 3-dimensional real Bloch vector corresponding to the qubit density matrix A...
Definition: functions.h:1762
Subsystems mismatch dimensions exception.
Definition: exception.h:365
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1171
Matrix is not a column vector exception.
Definition: exception.h:164
Quantum++ main namespace.
Definition: codes.h:35
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:267
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:242
cmat spectralpowm(const Eigen::MatrixBase< Derived > &A, const cplx z)
Matrix power.
Definition: functions.h:709
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:376
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:831
Invalid dimension(s) exception.
Definition: exception.h:269
Checks whether T is compatible with an STL-like iterable container.
Definition: traits.h:71
dyn_col_vect< typename Derived::Scalar > rho2pure(const Eigen::MatrixBase< Derived > &A)
Finds the pure state representation of a matrix proportional to a projector onto a pure state...
Definition: functions.h:1696
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1063
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1733
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:319
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:528
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:796
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:478
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1792
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1348
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:83
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1376
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:455
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:868
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1491
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:145
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:46
dyn_mat< typename Derived::Scalar > kronpow(const Eigen::MatrixBase< Derived > &A, idx n)
Kronecker power.
Definition: functions.h:944
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &A)
Projector.
Definition: functions.h:1210
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Fast matrix power based on the SQUARE-AND-MULTIPLY algorithm.
Definition: functions.h:751
Type mismatch exception.
Definition: exception.h:530
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:294
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:401
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:65
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:89
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:502
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:317
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:680
Parameter out of range exception.
Definition: exception.h:515
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:93
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:656
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of an STL-like range of complex numbers.
Definition: functions.h:1573
Matrix is not 2 x 2 exception.
Definition: exception.h:411
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &As)
Gram-Schmidt orthogonalization.
Definition: functions.h:1239
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1129
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1405
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:104
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:632
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:560
std::size_t idx
Non-negative integer index.
Definition: types.h:39
Matrix is not square exception.
Definition: exception.h:149
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:201
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1095
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:430
Object has zero size exception.
Definition: exception.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:236