Quantum++  v1.0.0-beta3
C++11 quantum computing library
functions.h
Go to the documentation of this file.
1 /*
2  * Quantum++
3  *
4  * Copyright (c) 2013 - 2017 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 FUNCTIONS_H_
28 #define FUNCTIONS_H_
29 
30 // Collection of generic quantum computing functions
31 
32 namespace qpp
33 {
34 
35 // Eigen function wrappers
43 template<typename Derived>
45  const Eigen::MatrixBase<Derived>& A)
46 {
47  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
48 
49  // check zero-size
51  throw Exception("qpp::transpose()", Exception::Type::ZERO_SIZE);
52 
53  return rA.transpose();
54 }
55 
63 template<typename Derived>
65  const Eigen::MatrixBase<Derived>& A)
66 {
67  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
68 
69  // check zero-size
71  throw Exception("qpp::conjugate()", Exception::Type::ZERO_SIZE);
72 
73  return rA.conjugate();
74 }
75 
83 template<typename Derived>
84 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A)
85 {
86  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
87 
88  // EXCEPTION CHECKS
89 
90  // check zero-size
92  throw Exception("qpp::adjoint()", Exception::Type::ZERO_SIZE);
93  // END EXCEPTION CHECKS
94 
95  return rA.adjoint();
96 }
97 
105 template<typename Derived>
106 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A)
107 {
108  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
109 
110  // EXCEPTION CHECKS
111 
112  // check zero-size
114  throw Exception("qpp::inverse()", Exception::Type::ZERO_SIZE);
115  // END EXCEPTION CHECKS
116 
117  return rA.inverse();
118 }
119 
126 template<typename Derived>
127 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A)
128 {
129  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
130 
131  // EXCEPTION CHECKS
132 
133  // check zero-size
135  throw Exception("qpp::trace()", Exception::Type::ZERO_SIZE);
136  // END EXCEPTION CHECKS
137 
138  return rA.trace();
139 }
140 
148 template<typename Derived>
149 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A)
150 {
151  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
152 
153  // EXCEPTION CHECKS
154 
155  // check zero-size
157  throw Exception("qpp::det()", Exception::Type::ZERO_SIZE);
158  // END EXCEPTION CHECKS
159 
160  return rA.determinant();
161 }
162 
172 template<typename Derived>
173 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A)
174 {
175  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
176 
177  // EXCEPTION CHECKS
178 
179  // check zero-size
181  throw Exception("qpp::logdet()", Exception::Type::ZERO_SIZE);
182 
183  // check square matrix
185  throw Exception("qpp::logdet()",
187  // END EXCEPTION CHECKS
188 
189  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
191  lu.matrixLU().template triangularView<Eigen::Upper>();
192  typename Derived::Scalar result = std::log(U(0, 0));
193 
194  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
195  result += std::log(U(i, i));
196 
197  return result;
198 
199 }
200 
208 template<typename Derived>
209 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A)
210 {
211  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
212 
213  // EXCEPTION CHECKS
214 
215  // check zero-size
217  throw Exception("qpp::sum()", Exception::Type::ZERO_SIZE);
218  // END EXCEPTION CHECKS
219 
220  return rA.sum();
221 }
222 
230 template<typename Derived>
231 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A)
232 {
233  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
234 
235  // EXCEPTION CHECKS
236 
237  // check zero-size
239  throw Exception("qpp::prod()", Exception::Type::ZERO_SIZE);
240  // END EXCEPTION CHECKS
241 
242  return rA.prod();
243 }
244 
251 template<typename Derived>
252 double norm(const Eigen::MatrixBase<Derived>& A)
253 {
254  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
255 
256  // EXCEPTION CHECKS
257 
258  // check zero-size
260  throw Exception("qpp::norm()", Exception::Type::ZERO_SIZE);
261  // END EXCEPTION CHECKS
262 
263  // convert matrix to complex then return its norm
264  return (rA.template cast<cplx>()).norm();
265 }
266 
275 template<typename Derived>
276 std::pair<dyn_col_vect<cplx>, cmat> eig(const Eigen::MatrixBase<Derived>& A)
277 {
278  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
279 
280  // EXCEPTION CHECKS
281 
282  // check zero-size
284  throw Exception("qpp::eig()", Exception::Type::ZERO_SIZE);
285 
286  // check square matrix
288  throw Exception("qpp::eig()", Exception::Type::MATRIX_NOT_SQUARE);
289  // END EXCEPTION CHECKS
290 
291  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
292 
293  return std::make_pair(es.eigenvalues(), es.eigenvectors());
294 }
295 
303 template<typename Derived>
304 dyn_col_vect <cplx> evals(const Eigen::MatrixBase<Derived>& A)
305 {
306  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
307 
308  // EXCEPTION CHECKS
309 
310  // check zero-size
312  throw Exception("qpp::evals()", Exception::Type::ZERO_SIZE);
313 
314  // check square matrix
316  throw Exception("qpp::evals()", Exception::Type::MATRIX_NOT_SQUARE);
317  // END EXCEPTION CHECKS
318 
319  return eig(rA).first;
320 }
321 
329 template<typename Derived>
330 cmat evects(const Eigen::MatrixBase<Derived>& A)
331 {
332  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
333 
334  // EXCEPTION CHECKS
335 
336  // check zero-size
338  throw Exception("qpp::evects()", Exception::Type::ZERO_SIZE);
339 
340  // check square matrix
342  throw Exception("qpp::evects()", Exception::Type::MATRIX_NOT_SQUARE);
343  // END EXCEPTION CHECKS
344 
345  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
346 
347  return eig(rA).second;
348 }
349 
358 template<typename Derived>
359 std::pair<dyn_col_vect<double>, cmat> heig(const Eigen::MatrixBase<Derived>& A)
360 {
361  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
362 
363  // EXCEPTION CHECKS
364 
365  // check zero-size
367  throw Exception("qpp::heig()", Exception::Type::ZERO_SIZE);
368 
369  // check square matrix
371  throw Exception("qpp::heig()", Exception::Type::MATRIX_NOT_SQUARE);
372  // END EXCEPTION CHECKS
373 
374  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
375 
376  return std::make_pair(es.eigenvalues(), es.eigenvectors());
377 }
378 
386 template<typename Derived>
387 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A)
388 {
389  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
390 
391  // EXCEPTION CHECKS
392 
393  // check zero-size
395  throw Exception("qpp::hevals()", Exception::Type::ZERO_SIZE);
396 
397  // check square matrix
399  throw Exception("qpp::hevals()", Exception::Type::MATRIX_NOT_SQUARE);
400  // END EXCEPTION CHECKS
401 
402  return heig(rA).first;
403 }
404 
412 template<typename Derived>
413 cmat hevects(const Eigen::MatrixBase<Derived>& A)
414 {
415  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
416 
417  // EXCEPTION CHECKS
418 
419  // check zero-size
421  throw Exception("qpp::hevects()", Exception::Type::ZERO_SIZE);
422 
423  // check square matrix
425  throw Exception("qpp::hevects()",
427  // END EXCEPTION CHECKS
428 
429  return heig(rA).second;
430 }
431 
441 template<typename Derived>
442 std::tuple<cmat, dyn_col_vect<double>, cmat>
443 svd(const Eigen::MatrixBase<Derived>& A)
444 {
445  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
446 
447  // EXCEPTION CHECKS
448 
449  // check zero-size
451  throw Exception("qpp::svd()", Exception::Type::ZERO_SIZE);
452  // END EXCEPTION CHECKS
453 
454  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
455  sv(rA,
456  Eigen::DecompositionOptions::ComputeFullU |
457  Eigen::DecompositionOptions::ComputeFullV);
458 
459  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
460 }
461 
469 template<typename Derived>
470 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A)
471 {
472  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
473 
474  // EXCEPTION CHECKS
475 
476  // check zero-size
478  throw Exception("qpp::svals()", Exception::Type::ZERO_SIZE);
479  // END EXCEPTION CHECKS
480 
481  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
482 
483  return sv.singularValues();
484 }
485 
493 template<typename Derived>
494 cmat svdU(const Eigen::MatrixBase<Derived>& A)
495 {
496  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
497 
498  // EXCEPTION CHECKS
499 
500  // check zero-size
502  throw Exception("qpp::svdU()", Exception::Type::ZERO_SIZE);
503  // END EXCEPTION CHECKS
504 
505  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
506  sv(rA, Eigen::DecompositionOptions::ComputeFullU);
507 
508  return sv.matrixU();
509 }
510 
518 template<typename Derived>
519 cmat svdV(const Eigen::MatrixBase<Derived>& A)
520 {
521  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
522 
523  // EXCEPTION CHECKS
524 
525  // check zero-size
527  throw Exception("qpp::svdV()", Exception::Type::ZERO_SIZE);
528  // END EXCEPTION CHECKS
529 
530  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
531  sv(rA, Eigen::DecompositionOptions::ComputeFullV);
532 
533  return sv.matrixV();
534 }
535 
536 // Matrix functional calculus
537 
545 template<typename Derived>
546 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (* f)(const cplx&))
547 {
548  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
549 
550  // EXCEPTION CHECKS
551 
552  // check zero-size
554  throw Exception("qpp::funm()", Exception::Type::ZERO_SIZE);
555 
556  // check square matrix
558  throw Exception("qpp::funm()", Exception::Type::MATRIX_NOT_SQUARE);
559  // END EXCEPTION CHECKS
560 
561  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
562  cmat evects = es.eigenvectors();
563  cmat evals = es.eigenvalues();
564  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
565  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
566 
567  cmat evalsdiag = evals.asDiagonal();
568 
569  return evects * evalsdiag * evects.inverse();
570 }
571 
578 template<typename Derived>
579 cmat sqrtm(const Eigen::MatrixBase<Derived>& A)
580 {
581  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
582 
583  // EXCEPTION CHECKS
584 
585  // check zero-size
587  throw Exception("qpp::sqrtm()", Exception::Type::ZERO_SIZE);
588 
589  // check square matrix
591  throw Exception("qpp::sqrtm()", Exception::Type::MATRIX_NOT_SQUARE);
592  // END EXCEPTION CHECKS
593 
594  return funm(rA, &std::sqrt);
595 }
596 
603 template<typename Derived>
604 cmat absm(const Eigen::MatrixBase<Derived>& A)
605 {
606  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
607 
608  // EXCEPTION CHECKS
609 
610  // check zero-size
612  throw Exception("qpp::absm()", Exception::Type::ZERO_SIZE);
613 
614  // check square matrix
616  throw Exception("qpp::absm()", Exception::Type::MATRIX_NOT_SQUARE);
617  // END EXCEPTION CHECKS
618 
619  return sqrtm(adjoint(rA) * rA);
620 }
621 
628 template<typename Derived>
629 cmat expm(const Eigen::MatrixBase<Derived>& A)
630 {
631  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
632 
633  // EXCEPTION CHECKS
634 
635  // check zero-size
637  throw Exception("qpp::expm()", Exception::Type::ZERO_SIZE);
638 
639  // check square matrix
641  throw Exception("qpp::expm()", Exception::Type::MATRIX_NOT_SQUARE);
642  // END EXCEPTION CHECKS
643 
644  return funm(rA, &std::exp);
645 }
646 
653 template<typename Derived>
654 cmat logm(const Eigen::MatrixBase<Derived>& A)
655 {
656  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
657 
658  // EXCEPTION CHECKS
659 
660  // check zero-size
662  throw Exception("qpp::logm()", Exception::Type::ZERO_SIZE);
663 
664  // check square matrix
666  throw Exception("qpp::logm()", Exception::Type::MATRIX_NOT_SQUARE);
667  // END EXCEPTION CHECKS
668 
669  return funm(rA, &std::log);
670 }
671 
678 template<typename Derived>
679 cmat sinm(const Eigen::MatrixBase<Derived>& A)
680 {
681  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
682 
683  // EXCEPTION CHECKS
684 
685  // check zero-size
687  throw Exception("qpp::sinm()", Exception::Type::ZERO_SIZE);
688 
689  // check square matrix
691  throw Exception("qpp::sinm()", Exception::Type::MATRIX_NOT_SQUARE);
692  // END EXCEPTION CHECKS
693 
694  return funm(rA, &std::sin);
695 }
696 
703 template<typename Derived>
704 cmat cosm(const Eigen::MatrixBase<Derived>& A)
705 {
706  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
707 
708  // EXCEPTION CHECKS
709 
710  // check zero-size
712  throw Exception("qpp::cosm()", Exception::Type::ZERO_SIZE);
713 
714  // check square matrix
716  throw Exception("qpp::cosm()", Exception::Type::MATRIX_NOT_SQUARE);
717  // END EXCEPTION CHECKS
718 
719  return funm(rA, &std::cos);
720 }
721 
733 template<typename Derived>
734 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z)
735 {
736  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
737 
738  // EXCEPTION CHECKS
739 
740  // check zero-size
742  throw Exception("qpp::spectralpowm()", Exception::Type::ZERO_SIZE);
743 
744  // check square matrix
746  throw Exception("qpp::spectralpowm()",
748  // END EXCEPTION CHECKS
749 
750  // Define A^0 = Id, for z IDENTICALLY zero
751  if (real(z) == 0 && imag(z) == 0)
752  return cmat::Identity(rA.rows(), rA.rows());
753 
754  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
755  cmat evects = es.eigenvectors();
756  cmat evals = es.eigenvalues();
757  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
758  evals(i) = std::pow(evals(i), z);
759 
760  cmat evalsdiag = evals.asDiagonal();
761 
762  return evects * evalsdiag * evects.inverse();
763 }
764 
777 template<typename Derived>
778 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
779  idx n)
780 {
781  // EXCEPTION CHECKS
782 
783  // check zero-size
785  throw Exception("qpp::powm()", Exception::Type::ZERO_SIZE);
786 
787  // check square matrix
789  throw Exception("qpp::powm()", Exception::Type::MATRIX_NOT_SQUARE);
790  // END EXCEPTION CHECKS
791 
792  // if n = 1, return the matrix unchanged
793  if (n == 1)
794  return A;
795 
798  A.rows(), A.rows());
799 
800  // if n = 0, return the identity (as just prepared in result)
801  if (n == 0)
802  return result;
803 
804  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
805 
806  // fast matrix power
807  for (; n > 0; n /= 2)
808  {
809  if (n % 2)
810  result = (result * cA).eval();
811  cA = (cA * cA).eval();
812  }
813 
814  return result;
815 }
816 
825 template<typename Derived>
826 double schatten(const Eigen::MatrixBase<Derived>& A, double p)
827 {
828  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
829 
830  // EXCEPTION CHECKS
831 
832  // check zero-size
834  throw Exception("qpp::schatten()", Exception::Type::ZERO_SIZE);
835  if (p < 1)
836  throw Exception("qpp::schatten()", Exception::Type::OUT_OF_RANGE);
837  // END EXCEPTION CHECKS
838 
839  if (p == infty) // infinity norm (largest singular value)
840  return svals(rA)(0);
841 
842  return std::pow(trace(powm(absm(rA), p)).real(), 1. / p);
843 }
844 
845 // other functions
846 
855 template<typename OutputScalar, typename Derived>
856 dyn_mat <OutputScalar> cwise(const Eigen::MatrixBase<Derived>& A,
857  OutputScalar (* f)(
858  const typename Derived::Scalar&))
859 {
860  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
861 
862  // EXCEPTION CHECKS
863 
864  // check zero-size
866  throw Exception("qpp::cwise()", Exception::Type::ZERO_SIZE);
867  // END EXCEPTION CHECKS
868 
869  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
870 
871 #ifdef WITH_OPENMP_
872 #pragma omp parallel for collapse(2)
873 #endif // WITH_OPENMP_
874  // column major order for speed
875  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
876  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
877  result(i, j) = (*f)(rA(i, j));
878 
879  return result;
880 }
881 
882 // Kronecker product of multiple matrices, preserve return type
883 // variadic template
894 template<typename T>
896 {
897  return head;
898 }
899 
910 template<typename T, typename ... Args>
911 dyn_mat<typename T::Scalar> kron(const T& head, const Args& ... tail)
912 {
913  return internal::kron2(head, kron(tail...));
914 }
915 
925 template<typename Derived>
926 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As)
927 {
928  // EXCEPTION CHECKS
929 
930  if (As.size() == 0)
931  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
932 
933  for (auto&& it : As)
935  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
936  // END EXCEPTION CHECKS
937 
938  dyn_mat<typename Derived::Scalar> result = As[0].derived();
939  for (idx i = 1; i < As.size(); ++i)
940  {
941  result = kron(result, As[i]);
942  }
943 
944  return result;
945 }
946 
947 // Kronecker product of a list of matrices, preserve return type
948 // deduce the template parameters from initializer_list
959 template<typename Derived>
961  const std::initializer_list<Derived>& As)
962 {
963  return kron(std::vector<Derived>(As));
964 }
965 
975 template<typename Derived>
976 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
977  idx n)
978 {
979  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
980 
981  // EXCEPTION CHECKS
982 
983  // check zero-size
985  throw Exception("qpp::kronpow()", Exception::Type::ZERO_SIZE);
986 
987  // check out of range
988  if (n == 0)
989  throw Exception("qpp::kronpow()", Exception::Type::OUT_OF_RANGE);
990  // END EXCEPTION CHECKS
991 
992  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
993 
994  return kron(As);
995 }
996 
997 // Direct sum of multiple matrices, preserve return type
998 // variadic template
1009 template<typename T>
1011 {
1012  return head;
1013 }
1014 
1025 template<typename T, typename ... Args>
1026 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args& ... tail)
1027 {
1028  return internal::dirsum2(head, dirsum(tail...));
1029 }
1030 
1040 template<typename Derived>
1041 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As)
1042 {
1043  // EXCEPTION CHECKS
1044 
1045  if (As.size() == 0)
1046  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1047 
1048  for (auto&& it : As)
1050  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1051  // END EXCEPTION CHECKS
1052 
1053  idx total_rows = 0, total_cols = 0;
1054  for (idx i = 0; i < As.size(); ++i)
1055  {
1056  total_rows += static_cast<idx>(As[i].rows());
1057  total_cols += static_cast<idx>(As[i].cols());
1058  }
1060  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1061 
1062  idx cur_row = 0, cur_col = 0;
1063  for (idx i = 0; i < As.size(); ++i)
1064  {
1065  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1066  cur_row += static_cast<idx>(As[i].rows());
1067  cur_col += static_cast<idx>(As[i].cols());
1068  }
1069 
1070  return result;
1071 }
1072 
1073 // Direct sum of a list of matrices, preserve return type
1074 // deduce the template parameters from initializer_list
1085 template<typename Derived>
1087  const std::initializer_list<Derived>& As)
1088 {
1089  return dirsum(std::vector<Derived>(As));
1090 }
1091 
1101 template<typename Derived>
1103  const Eigen::MatrixBase<Derived>& A, idx n)
1104 {
1105  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1106 
1107  // EXCEPTION CHECKS
1108 
1109  // check zero-size
1111  throw Exception("qpp::dirsumpow()", Exception::Type::ZERO_SIZE);
1112 
1113  // check out of range
1114  if (n == 0)
1115  throw Exception("qpp::dirsumpow()", Exception::Type::OUT_OF_RANGE);
1116  // END EXCEPTION CHECKS
1117 
1118  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1119 
1120  return dirsum(As);
1121 }
1122 
1134 template<typename Derived>
1135 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1136  idx rows, idx cols)
1137 {
1138  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1139 
1140  idx Arows = static_cast<idx>(rA.rows());
1141  idx Acols = static_cast<idx>(rA.cols());
1142 
1143  // EXCEPTION CHECKS
1144 
1145  // check zero-size
1147  throw Exception("qpp::reshape()", Exception::Type::ZERO_SIZE);
1148 
1149  if (Arows * Acols != rows * cols)
1150  throw Exception("qpp::reshape()",
1152  // END EXCEPTION CHECKS
1153 
1154  return Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1155  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1156 }
1157 
1170 template<typename Derived1, typename Derived2>
1171 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1172  const Eigen::MatrixBase<Derived2>& B)
1173 {
1174  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1175  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1176 
1177  // EXCEPTION CHECKS
1178 
1179  // check types
1180  if (!std::is_same<typename Derived1::Scalar,
1181  typename Derived2::Scalar>::value)
1182  throw Exception("qpp::comm()", Exception::Type::TYPE_MISMATCH);
1183 
1184  // check zero-size
1187  throw Exception("qpp::comm()", Exception::Type::ZERO_SIZE);
1188 
1189  // check square matrices
1191  throw Exception("qpp::comm()", Exception::Type::MATRIX_NOT_SQUARE);
1192 
1193  // check equal dimensions
1194  if (rA.rows() != rB.rows())
1195  throw Exception("qpp::comm()", Exception::Type::DIMS_NOT_EQUAL);
1196  // END EXCEPTION CHECKS
1197 
1198  return rA * rB - rB * rA;
1199 }
1200 
1213 template<typename Derived1, typename Derived2>
1215  const Eigen::MatrixBase<Derived1>& A,
1216  const Eigen::MatrixBase<Derived2>& B)
1217 {
1218  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1219  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1220 
1221  // EXCEPTION CHECKS
1222 
1223  // check types
1224  if (!std::is_same<typename Derived1::Scalar,
1225  typename Derived2::Scalar>::value)
1226  throw Exception("qpp::anticomm()", Exception::Type::TYPE_MISMATCH);
1227 
1228  // check zero-size
1231  throw Exception("qpp::anticomm()", Exception::Type::ZERO_SIZE);
1232 
1233  // check square matrices
1235  throw Exception("qpp::anticomm()",
1237 
1238  // check equal dimensions
1239  if (rA.rows() != rB.rows())
1240  throw Exception("qpp::anticomm()", Exception::Type::DIMS_NOT_EQUAL);
1241  // END EXCEPTION CHECKS
1242 
1243  return rA * rB + rB * rA;
1244 }
1245 
1256 template<typename Derived>
1257 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A)
1258 {
1259  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1260 
1261  // EXCEPTION CHECKS
1262 
1263  // check zero-size
1265  throw Exception("qpp::prj()", Exception::Type::ZERO_SIZE);
1266 
1267  // check column vector
1268  if (!internal::check_cvector(rA))
1269  throw Exception("qpp::prj()", Exception::Type::MATRIX_NOT_CVECTOR);
1270  // END EXCEPTION CHECKS
1271 
1272  double normA = norm(rA);
1273  if (normA > eps)
1274  return rA * adjoint(rA) / (normA * normA);
1275  else
1277  ::Zero(rA.rows(), rA.rows());
1278 
1279 }
1280 
1288 template<typename Derived>
1289 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As)
1290 {
1291  // EXCEPTION CHECKS
1292 
1293  // check empty list
1295  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1296 
1297  for (auto&& it : As)
1299  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1300 
1301  // check that As[0] is a column vector
1302  if (!internal::check_cvector(As[0]))
1303  throw Exception("qpp::grams()", Exception::Type::MATRIX_NOT_CVECTOR);
1304 
1305  // now check that all the rest match the size of the first vector
1306  for (auto&& it : As)
1307  if (it.rows() != As[0].rows() || it.cols() != 1)
1308  throw Exception("qpp::grams()", Exception::Type::DIMS_NOT_EQUAL);
1309  // END EXCEPTION CHECKS
1310 
1313  As[0].rows());
1314 
1316  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1317 
1318  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1319  // find the first non-zero vector in the list
1320  idx pos = 0;
1321  for (pos = 0; pos < As.size(); ++pos)
1322  {
1323  if (norm(As[pos]) > eps) // add it as the first element
1324  {
1325  outvecs.push_back(As[pos]);
1326  break;
1327  }
1328  }
1329 
1330  // start the process
1331  for (idx i = pos + 1; i < As.size(); ++i)
1332  {
1333  cut -= prj(outvecs[i - 1 - pos]);
1334  vi = cut * As[i];
1335  outvecs.push_back(vi);
1336  }
1337 
1338  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1339 
1340  idx cnt = 0;
1341  for (auto&& it : outvecs)
1342  {
1343  double normA = norm(it);
1344  if (normA > eps) // we add only the non-zero vectors
1345  {
1346  result.col(cnt) = it / normA;
1347  cnt++;
1348  }
1349  }
1350 
1351  return result.block(0, 0, As[0].rows(), cnt);
1352 }
1353 
1354 // deduce the template parameters from initializer_list
1362 template<typename Derived>
1364  const std::initializer_list<Derived>& As)
1365 {
1366  return grams(std::vector<Derived>(As));
1367 }
1368 
1376 template<typename Derived>
1377 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A)
1378 {
1379  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1380 
1381  // EXCEPTION CHECKS
1382 
1384  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1385  // END EXCEPTION CHECKS
1386 
1387  std::vector<dyn_mat<typename Derived::Scalar>> input;
1388 
1389  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1390  input.push_back(rA.col(i));
1391 
1392  return grams<dyn_mat<typename Derived::Scalar>>
1393  (input);
1394 }
1395 
1406 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims)
1407 {
1408  // EXCEPTION CHECKS
1409 
1410  if (!internal::check_dims(dims))
1411  throw Exception("qpp::n2multiidx()", Exception::Type::DIMS_INVALID);
1412 
1413  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1414  static_cast<idx>(1), std::multiplies<idx>()))
1415  throw Exception("qpp::n2multiidx()", Exception::Type::OUT_OF_RANGE);
1416  // END EXCEPTION CHECKS
1417 
1418  // double the size for matrices reshaped as vectors
1419  idx result[2 * maxn];
1420  internal::n2multiidx(n, dims.size(), dims.data(), result);
1421 
1422  return std::vector<idx>(result, result + dims.size());
1423 }
1424 
1435 inline idx multiidx2n(const std::vector<idx>& midx,
1436  const std::vector<idx>& dims)
1437 {
1438  // EXCEPTION CHECKS
1439 
1440  if (!internal::check_dims(dims))
1441  throw Exception("qpp::multiidx2n()", Exception::Type::DIMS_INVALID);
1442 
1443  for (idx i = 0; i < dims.size(); ++i)
1444  if (midx[i] >= dims[i])
1445  throw Exception("qpp::multiidx2n()",
1447  // END EXCEPTION CHECKS
1448 
1449  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1450 }
1451 
1464 inline ket mket(const std::vector<idx>& mask,
1465  const std::vector<idx>& dims)
1466 {
1467  idx N = mask.size();
1468 
1469  idx D = std::accumulate(std::begin(dims), std::end(dims),
1470  static_cast<idx>(1), std::multiplies<idx>());
1471 
1472  // EXCEPTION CHECKS
1473 
1474  // check zero size
1475  if (N == 0)
1476  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1477  // check valid dims
1478  if (!internal::check_dims(dims))
1479  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1480  // check mask and dims have the same size
1481  if (mask.size() != dims.size())
1482  throw Exception("qpp::mket()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1483  // check mask is a valid vector
1484  for (idx i = 0; i < N; ++i)
1485  if (mask[i] >= dims[i])
1486  throw Exception("qpp::mket()",
1488  // END EXCEPTION CHECKS
1489 
1490  ket result = ket::Zero(D);
1491  idx pos = multiidx2n(mask, dims);
1492  result(pos) = 1;
1493 
1494  return result;
1495 }
1496 
1509 inline ket mket(const std::vector<idx>& mask, idx d = 2)
1510 {
1511  idx N = mask.size();
1512  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1513 
1514  // EXCEPTION CHECKS
1515 
1516  // check zero size
1517  if (N == 0)
1518  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1519 
1520  // check valid dims
1521  if (d == 0)
1522  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1523 
1524  // check mask is a valid vector
1525  for (idx i = 0; i < N; ++i)
1526  if (mask[i] >= d)
1527  throw Exception("qpp::mket()",
1529  // END EXCEPTION CHECKS
1530 
1531  ket result = ket::Zero(D);
1532  std::vector<idx> dims(N, d);
1533  idx pos = multiidx2n(mask, dims);
1534  result(pos) = 1;
1535 
1536  return result;
1537 }
1538 
1553 inline cmat mprj(const std::vector<idx>& mask,
1554  const std::vector<idx>& dims)
1555 {
1556  idx N = mask.size();
1557 
1558  idx D = std::accumulate(std::begin(dims), std::end(dims),
1559  static_cast<idx>(1), std::multiplies<idx>());
1560 
1561  // EXCEPTION CHECKS
1562 
1563  // check zero size
1564  if (N == 0)
1565  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1566  // check valid dims
1567  if (!internal::check_dims(dims))
1568  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1569  // check mask and dims have the same size
1570  if (mask.size() != dims.size())
1571  throw Exception("qpp::mprj()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1572  // check mask is a valid vector
1573  for (idx i = 0; i < N; ++i)
1574  if (mask[i] >= dims[i])
1575  throw Exception("qpp::mprj()",
1577  // END EXCEPTION CHECKS
1578 
1579  cmat result = cmat::Zero(D, D);
1580  idx pos = multiidx2n(mask, dims);
1581  result(pos, pos) = 1;
1582 
1583  return result;
1584 }
1585 
1600 inline cmat mprj(const std::vector<idx>& mask, idx d = 2)
1601 {
1602  idx N = mask.size();
1603  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1604 
1605  // EXCEPTION CHECKS
1606 
1607  // check zero size
1608  if (N == 0)
1609  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1610 
1611  // check valid dims
1612  if (d == 0)
1613  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1614 
1615  // check mask is a valid vector
1616  for (idx i = 0; i < N; ++i)
1617  if (mask[i] >= d)
1618  throw Exception("qpp::mprj()",
1620  // END EXCEPTION CHECKS
1621 
1622  cmat result = cmat::Zero(D, D);
1623  std::vector<idx> dims(N, d);
1624  idx pos = multiidx2n(mask, dims);
1625  result(pos, pos) = 1;
1626 
1627  return result;
1628 }
1629 
1638 template<typename InputIterator>
1639 std::vector<double> abssq(InputIterator first, InputIterator last)
1640 {
1641  std::vector<double> weights(std::distance(first, last));
1642  std::transform(first, last, std::begin(weights),
1643  [](cplx z) -> double
1644  {
1645  return std::norm(z);
1646  });
1647 
1648  return weights;
1649 }
1650 
1657 template<typename Container>
1658 std::vector<double> abssq(const Container& c,
1659  typename std::enable_if<
1661  >::type* = nullptr)
1662 // we need the std::enable_if to SFINAE out Eigen expressions
1663 // that will otherwise match, instead of matching
1664 // the overload below:
1665 
1666 // template<typename Derived>
1667 // abssq(const Eigen::MatrixBase<Derived>& A)
1668 {
1669  return abssq(std::begin(c), std::end(c));
1670 }
1671 
1678 template<typename Derived>
1679 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A)
1680 {
1681  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1682 
1683  // EXCEPTION CHECKS
1684 
1685  // check zero-size
1687  throw Exception("qpp::abssq()", Exception::Type::ZERO_SIZE);
1688  // END EXCEPTION CHECKS
1689 
1690  return abssq(rA.data(), rA.data() + rA.size());
1691 }
1692 
1701 template<typename InputIterator>
1702 typename std::iterator_traits<InputIterator>::value_type
1703 sum(InputIterator first, InputIterator last)
1704 {
1705  using value_type =
1706  typename std::iterator_traits<InputIterator>::value_type;
1707 
1708  return std::accumulate(first, last, static_cast<value_type>(0));
1709 }
1710 
1718 template<typename Container>
1719 typename Container::value_type
1720 sum(const Container& c,
1721  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1722 {
1723  return sum(std::begin(c), std::end(c));
1724 }
1725 
1734 template<typename InputIterator>
1735 typename std::iterator_traits<InputIterator>::value_type
1736 prod(InputIterator first, InputIterator last)
1737 {
1738  using value_type =
1739  typename std::iterator_traits<InputIterator>::value_type;
1740 
1741  return std::accumulate(first, last, static_cast<value_type>(1),
1742  std::multiplies<value_type>());
1743 }
1744 
1752 template<typename Container>
1753 typename Container::value_type
1754 prod(const Container& c,
1755  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1756 {
1757  return prod(std::begin(c), std::end(c));
1758 }
1759 
1772 template<typename Derived>
1774  const Eigen::MatrixBase<Derived>& A)
1775 {
1776  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1777 
1778  // EXCEPTION CHECKS
1779  // check zero-size
1781  throw Exception("qpp::rho2pure()", Exception::Type::ZERO_SIZE);
1782  // check square matrix
1783  if (!internal::check_square_mat(rA))
1784  throw Exception("qpp::rho2pure()", Exception::Type::MATRIX_NOT_SQUARE);
1785  // END EXCEPTION CHECKS
1786 
1787  dyn_col_vect<double> tmp_evals = hevals(rA);
1788  cmat tmp_evects = hevects(rA);
1791  // find the non-zero eigenvector
1792  // there is only one, assuming the state is pure
1793  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1794  {
1795  if (std::abs(tmp_evals(k)) > eps)
1796  {
1797  result = tmp_evects.col(k);
1798  break;
1799  }
1800  }
1801 
1802  return result;
1803 }
1804 
1813 template<typename T>
1814 std::vector<T> complement(std::vector<T> subsys, idx N)
1815 {
1816  // EXCEPTION CHECKS
1817 
1818  if (N < subsys.size())
1819  throw Exception("qpp::complement()", Exception::Type::OUT_OF_RANGE);
1820  // END EXCEPTION CHECKS
1821 
1822  std::vector<T> all(N);
1823  std::vector<T> subsys_bar(N - subsys.size());
1824 
1825  std::iota(std::begin(all), std::end(all), 0);
1826  std::sort(std::begin(subsys), std::end(subsys));
1827  std::set_difference(std::begin(all), std::end(all),
1828  std::begin(subsys), std::end(subsys),
1829  std::begin(subsys_bar));
1830 
1831  return subsys_bar;
1832 }
1833 
1844 template<typename Derived>
1845 std::vector<double> rho2bloch(
1846  const Eigen::MatrixBase<Derived>& A)
1847 {
1848  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1849 
1850  // EXCEPTION CHECKS
1851 
1852  // check qubit matrix
1854  throw Exception("qpp::rho2bloch()", Exception::Type::NOT_QUBIT_MATRIX);
1855  // END EXCEPTION CHECKS
1856 
1857  std::vector<double> result(3);
1858  cmat X(2, 2), Y(2, 2), Z(2, 2);
1859  X << 0, 1, 1, 0;
1860  Y << 0, -1_i, 1_i, 0;
1861  Z << 1, 0, 0, -1;
1862  result[0] = std::real(trace(rA * X));
1863  result[1] = std::real(trace(rA * Y));
1864  result[2] = std::real(trace(rA * Z));
1865 
1866  return result;
1867 }
1868 
1877 inline cmat bloch2rho(const std::vector<double>& r)
1878 {
1879  // EXCEPTION CHECKS
1880 
1881  // check 3-dimensional vector
1882  if (r.size() != 3)
1883  throw Exception("qpp::bloch2rho", "r is not a 3-dimensional vector!");
1884  // END EXCEPTION CHECKS
1885 
1886  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1887  X << 0, 1, 1, 0;
1888  Y << 0, -1_i, 1_i, 0;
1889  Z << 1, 0, 0, -1;
1890  Id2 << 1, 0, 0, 1;
1891 
1892  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1893 }
1894 
1895 } /* namespace qpp */
1896 
1897 #endif /* FUNCTIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:127
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:359
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:173
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:74
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)( const typename Derived::Scalar &))
Functor.
Definition: functions.h:856
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolute value.
Definition: functions.h:604
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:272
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:120
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:629
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:231
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:67
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:1010
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:78
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:48
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:1845
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:51
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1214
Quantum++ main namespace.
Definition: codes.h:30
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:276
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:252
cmat spectralpowm(const Eigen::MatrixBase< Derived > &A, const cplx z)
Matrix power.
Definition: functions.h:734
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:387
Checks whether T is compatible with an STL-like iterable container.
Definition: traits.h:52
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:1773
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1102
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1814
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:330
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:546
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:826
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:494
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1877
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1406
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:84
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1435
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:470
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:895
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:46
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1553
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:149
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:44
dyn_mat< typename Derived::Scalar > kronpow(const Eigen::MatrixBase< Derived > &A, idx n)
Kronecker power.
Definition: functions.h:976
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &A)
Projector.
Definition: functions.h:1257
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:778
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:304
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:413
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:64
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:88
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:127
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:519
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:315
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:704
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:90
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:679
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of an STL-like range of complex numbers.
Definition: functions.h:1639
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &As)
Gram-Schmidt orthogonalization.
Definition: functions.h:1289
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1171
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:61
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1464
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:106
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:654
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:140
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:579
std::size_t idx
Non-negative integer index.
Definition: types.h:36
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:1135
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:443
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:228