Quantum++  v0.8.6
C++11 quantum computing library
functions.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 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 {
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 {
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 {
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;
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;
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;
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;
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;
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;
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;
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>
277 
278 eig(const Eigen::MatrixBase<Derived>& A)
279 {
280  const dyn_mat<typename Derived::Scalar>& rA = A;
281 
282  // EXCEPTION CHECKS
283 
284  // check zero-size
286  throw Exception("qpp::eig()", Exception::Type::ZERO_SIZE);
287 
288  // check square matrix
290  throw Exception("qpp::eig()", Exception::Type::MATRIX_NOT_SQUARE);
291  // END EXCEPTION CHECKS
292 
293  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
294 
295  return std::make_pair(es.eigenvalues(), es.eigenvectors());
296 }
297 
305 template<typename Derived>
306 dyn_col_vect <cplx> evals(const Eigen::MatrixBase<Derived>& A)
307 {
308  const dyn_mat<typename Derived::Scalar>& rA = A;
309 
310  // EXCEPTION CHECKS
311 
312  // check zero-size
314  throw Exception("qpp::evals()", Exception::Type::ZERO_SIZE);
315 
316  // check square matrix
318  throw Exception("qpp::evals()", Exception::Type::MATRIX_NOT_SQUARE);
319  // END EXCEPTION CHECKS
320 
321  return eig(rA).first;
322 }
323 
331 template<typename Derived>
332 cmat evects(const Eigen::MatrixBase<Derived>& A)
333 {
334  const dyn_mat<typename Derived::Scalar>& rA = A;
335 
336  // EXCEPTION CHECKS
337 
338  // check zero-size
340  throw Exception("qpp::evects()", Exception::Type::ZERO_SIZE);
341 
342  // check square matrix
344  throw Exception("qpp::evects()", Exception::Type::MATRIX_NOT_SQUARE);
345  // END EXCEPTION CHECKS
346 
347  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
348 
349  return eig(rA).second;
350 }
351 
360 template<typename Derived>
361 std::pair<dyn_col_vect < double>, cmat>
362 
363 heig(const Eigen::MatrixBase<Derived>& A)
364 {
365  const dyn_mat<typename Derived::Scalar>& rA = A;
366 
367  // EXCEPTION CHECKS
368 
369  // check zero-size
371  throw Exception("qpp::heig()", Exception::Type::ZERO_SIZE);
372 
373  // check square matrix
375  throw Exception("qpp::heig()", Exception::Type::MATRIX_NOT_SQUARE);
376  // END EXCEPTION CHECKS
377 
378  Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
379 
380  return std::make_pair(es.eigenvalues(), es.eigenvectors());
381 }
382 
390 template<typename Derived>
391 dyn_col_vect<double> hevals(const Eigen::MatrixBase<Derived>& A)
392 {
393  const dyn_mat<typename Derived::Scalar>& rA = A;
394 
395  // EXCEPTION CHECKS
396 
397  // check zero-size
399  throw Exception("qpp::hevals()", Exception::Type::ZERO_SIZE);
400 
401  // check square matrix
403  throw Exception("qpp::hevals()", Exception::Type::MATRIX_NOT_SQUARE);
404  // END EXCEPTION CHECKS
405 
406  return heig(rA).first;
407 }
408 
416 template<typename Derived>
417 cmat hevects(const Eigen::MatrixBase<Derived>& A)
418 {
419  const dyn_mat<typename Derived::Scalar>& rA = A;
420 
421  // EXCEPTION CHECKS
422 
423  // check zero-size
425  throw Exception("qpp::hevects()", Exception::Type::ZERO_SIZE);
426 
427  // check square matrix
429  throw Exception("qpp::hevects()",
431  // END EXCEPTION CHECKS
432 
433  return heig(rA).second;
434 }
435 
445 template<typename Derived>
446 std::tuple<cmat, dyn_col_vect < double>, cmat>
447 
448 svd(const Eigen::MatrixBase<Derived>& A)
449 {
450  const dyn_mat<typename Derived::Scalar>& rA = A;
451 
452  // EXCEPTION CHECKS
453 
454  // check zero-size
456  throw Exception("qpp::svd()", Exception::Type::ZERO_SIZE);
457  // END EXCEPTION CHECKS
458 
459  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
460  sv(rA,
461  Eigen::DecompositionOptions::ComputeFullU |
462  Eigen::DecompositionOptions::ComputeFullV);
463 
464  return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
465 }
466 
474 template<typename Derived>
475 dyn_col_vect<double> svals(const Eigen::MatrixBase<Derived>& A)
476 {
477  const dyn_mat<typename Derived::Scalar>& rA = A;
478 
479  // EXCEPTION CHECKS
480 
481  // check zero-size
483  throw Exception("qpp::svals()", Exception::Type::ZERO_SIZE);
484  // END EXCEPTION CHECKS
485 
486  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
487 
488  return sv.singularValues();
489 }
490 
498 template<typename Derived>
499 cmat svdU(const Eigen::MatrixBase<Derived>& A)
500 {
501  const dyn_mat<typename Derived::Scalar>& rA = A;
502 
503  // EXCEPTION CHECKS
504 
505  // check zero-size
507  throw Exception("qpp::svdU()", Exception::Type::ZERO_SIZE);
508  // END EXCEPTION CHECKS
509 
510  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
511  sv(rA, Eigen::DecompositionOptions::ComputeFullU);
512 
513  return sv.matrixU();
514 }
515 
523 template<typename Derived>
524 cmat svdV(const Eigen::MatrixBase<Derived>& A)
525 {
526  const dyn_mat<typename Derived::Scalar>& rA = A;
527 
528  // EXCEPTION CHECKS
529 
530  // check zero-size
532  throw Exception("qpp::svdV()", Exception::Type::ZERO_SIZE);
533  // END EXCEPTION CHECKS
534 
535  Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
536  sv(rA, Eigen::DecompositionOptions::ComputeFullV);
537 
538  return sv.matrixV();
539 }
540 
541 // Matrix functional calculus
542 
550 template<typename Derived>
551 cmat funm(const Eigen::MatrixBase<Derived>& A, cplx (* f)(const cplx&))
552 {
553  const dyn_mat<typename Derived::Scalar>& rA = A;
554 
555  // EXCEPTION CHECKS
556 
557  // check zero-size
559  throw Exception("qpp::funm()", Exception::Type::ZERO_SIZE);
560 
561  // check square matrix
563  throw Exception("qpp::funm()", Exception::Type::MATRIX_NOT_SQUARE);
564  // END EXCEPTION CHECKS
565 
566  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
567  cmat evects = es.eigenvectors();
568  cmat evals = es.eigenvalues();
569  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
570  evals(i) = (*f)(evals(i)); // apply f(x) to each eigenvalue
571 
572  cmat evalsdiag = evals.asDiagonal();
573 
574  return evects * evalsdiag * evects.inverse();
575 }
576 
583 template<typename Derived>
584 cmat sqrtm(const Eigen::MatrixBase<Derived>& A)
585 {
586  const dyn_mat<typename Derived::Scalar>& rA = A;
587 
588  // EXCEPTION CHECKS
589 
590  // check zero-size
592  throw Exception("qpp::sqrtm()", Exception::Type::ZERO_SIZE);
593 
594  // check square matrix
596  throw Exception("qpp::sqrtm()", Exception::Type::MATRIX_NOT_SQUARE);
597  // END EXCEPTION CHECKS
598 
599  return funm(rA, &std::sqrt);
600 }
601 
608 template<typename Derived>
609 cmat absm(const Eigen::MatrixBase<Derived>& A)
610 {
611  const dyn_mat<typename Derived::Scalar>& rA = A;
612 
613  // EXCEPTION CHECKS
614 
615  // check zero-size
617  throw Exception("qpp::absm()", Exception::Type::ZERO_SIZE);
618 
619  // check square matrix
621  throw Exception("qpp::absm()", Exception::Type::MATRIX_NOT_SQUARE);
622  // END EXCEPTION CHECKS
623 
624  return sqrtm(adjoint(rA) * rA);
625 }
626 
633 template<typename Derived>
634 cmat expm(const Eigen::MatrixBase<Derived>& A)
635 {
636  const dyn_mat<typename Derived::Scalar>& rA = A;
637 
638  // EXCEPTION CHECKS
639 
640  // check zero-size
642  throw Exception("qpp::expm()", Exception::Type::ZERO_SIZE);
643 
644  // check square matrix
646  throw Exception("qpp::expm()", Exception::Type::MATRIX_NOT_SQUARE);
647  // END EXCEPTION CHECKS
648 
649  return funm(rA, &std::exp);
650 }
651 
658 template<typename Derived>
659 cmat logm(const Eigen::MatrixBase<Derived>& A)
660 {
661  const dyn_mat<typename Derived::Scalar>& rA = A;
662 
663  // EXCEPTION CHECKS
664 
665  // check zero-size
667  throw Exception("qpp::logm()", Exception::Type::ZERO_SIZE);
668 
669  // check square matrix
671  throw Exception("qpp::logm()", Exception::Type::MATRIX_NOT_SQUARE);
672  // END EXCEPTION CHECKS
673 
674  return funm(rA, &std::log);
675 }
676 
683 template<typename Derived>
684 cmat sinm(const Eigen::MatrixBase<Derived>& A)
685 {
686  const dyn_mat<typename Derived::Scalar>& rA = A;
687 
688  // EXCEPTION CHECKS
689 
690  // check zero-size
692  throw Exception("qpp::sinm()", Exception::Type::ZERO_SIZE);
693 
694  // check square matrix
696  throw Exception("qpp::sinm()", Exception::Type::MATRIX_NOT_SQUARE);
697  // END EXCEPTION CHECKS
698 
699  return funm(rA, &std::sin);
700 }
701 
708 template<typename Derived>
709 cmat cosm(const Eigen::MatrixBase<Derived>& A)
710 {
711  const dyn_mat<typename Derived::Scalar>& rA = A;
712 
713  // EXCEPTION CHECKS
714 
715  // check zero-size
717  throw Exception("qpp::cosm()", Exception::Type::ZERO_SIZE);
718 
719  // check square matrix
721  throw Exception("qpp::cosm()", Exception::Type::MATRIX_NOT_SQUARE);
722  // END EXCEPTION CHECKS
723 
724  return funm(rA, &std::cos);
725 }
726 
738 template<typename Derived>
739 cmat spectralpowm(const Eigen::MatrixBase<Derived>& A, const cplx z)
740 {
741  const dyn_mat<typename Derived::Scalar>& rA = A;
742 
743  // EXCEPTION CHECKS
744 
745  // check zero-size
747  throw Exception("qpp::spectralpowm()", Exception::Type::ZERO_SIZE);
748 
749  // check square matrix
751  throw Exception("qpp::spectralpowm()",
753  // END EXCEPTION CHECKS
754 
755  // Define A^0 = Id, for z IDENTICALLY zero
756  if (real(z) == 0 && imag(z) == 0)
757  return cmat::Identity(rA.rows(), rA.rows());
758 
759  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
760  cmat evects = es.eigenvectors();
761  cmat evals = es.eigenvalues();
762  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
763  evals(i) = std::pow(evals(i), z);
764 
765  cmat evalsdiag = evals.asDiagonal();
766 
767  return evects * evalsdiag * evects.inverse();
768 }
769 
782 template<typename Derived>
783 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
784  idx n)
785 {
787 
788  // EXCEPTION CHECKS
789 
790  // check zero-size
792  throw Exception("qpp::powm()", Exception::Type::ZERO_SIZE);
793 
794  // check square matrix
796  throw Exception("qpp::powm()", Exception::Type::MATRIX_NOT_SQUARE);
797  // END EXCEPTION CHECKS
798 
801  cA.rows(), cA.rows());
802 
803  // fast matrix power
804  for (; n > 0; n /= 2)
805  {
806  if (n % 2)
807  result = (result * cA).eval();
808  cA = (cA * cA).eval();
809  }
810 
811  return result;
812 }
813 
822 template<typename Derived>
823 double schatten(const Eigen::MatrixBase<Derived>& A, double p)
824 {
825  const dyn_mat<typename Derived::Scalar>& rA = A;
826 
827  // EXCEPTION CHECKS
828 
829  // check zero-size
831  throw Exception("qpp::schatten()", Exception::Type::ZERO_SIZE);
832  if (p < 1)
833  throw Exception("qpp::schatten()", Exception::Type::OUT_OF_RANGE);
834  // END EXCEPTION CHECKS
835 
836  if (p == infty) // infinity norm (largest singular value)
837  return svals(rA)(0);
838 
839  return std::pow(trace(powm(absm(rA), p)).real(), 1. / p);
840 }
841 
842 // other functions
843 
852 template<typename OutputScalar, typename Derived>
853 dyn_mat <OutputScalar> cwise(const Eigen::MatrixBase<Derived>& A,
854  OutputScalar (* f)(
855  const typename Derived::Scalar&))
856 {
857  const dyn_mat<typename Derived::Scalar>& rA = A;
858 
859  // EXCEPTION CHECKS
860 
861  // check zero-size
863  throw Exception("qpp::cwise()", Exception::Type::ZERO_SIZE);
864  // END EXCEPTION CHECKS
865 
866  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
867 
868 #pragma omp parallel for collapse(2)
869  // column major order for speed
870  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
871  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
872  result(i, j) = (*f)(rA(i, j));
873 
874  return result;
875 }
876 
877 // Kronecker product of multiple matrices, preserve return type
878 // variadic template
889 template<typename T>
891 {
892  return head;
893 }
894 
905 template<typename T, typename ... Args>
906 dyn_mat<typename T::Scalar> kron(const T& head, const Args& ... tail)
907 {
908  return internal::_kron2(head, kron(tail...));
909 }
910 
920 template<typename Derived>
921 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As)
922 {
923  // EXCEPTION CHECKS
924 
925  if (As.size() == 0)
926  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
927 
928  for (auto&& it : As)
930  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
931  // END EXCEPTION CHECKS
932 
933  dyn_mat<typename Derived::Scalar> result = As[0];
934  for (idx i = 1; i < As.size(); ++i)
935  {
936  result = kron(result, As[i]);
937  }
938 
939  return result;
940 }
941 
942 // Kronecker product of a list of matrices, preserve return type
943 // deduce the template parameters from initializer_list
954 template<typename Derived>
956  const std::initializer_list<Derived>& As)
957 {
958  return kron(std::vector<Derived>(As));
959 }
960 
970 template<typename Derived>
971 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
972  idx n)
973 {
974  const dyn_mat<typename Derived::Scalar>& rA = A;
975 
976  // EXCEPTION CHECKS
977 
978  // check zero-size
980  throw Exception("qpp::kronpow()", Exception::Type::ZERO_SIZE);
981 
982  // check out of range
983  if (n == 0)
984  throw Exception("qpp::kronpow()", Exception::Type::OUT_OF_RANGE);
985  // END EXCEPTION CHECKS
986 
987  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
988 
989  return kron(As);
990 }
991 
992 // Direct sum of multiple matrices, preserve return type
993 // variadic template
1004 template<typename T>
1006 {
1007  return head;
1008 }
1009 
1020 template<typename T, typename ... Args>
1021 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args& ... tail)
1022 {
1023  return internal::_dirsum2(head, dirsum(tail...));
1024 }
1025 
1035 template<typename Derived>
1036 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As)
1037 {
1038  // EXCEPTION CHECKS
1039 
1040  if (As.size() == 0)
1041  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1042 
1043  for (auto&& it : As)
1045  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1046  // END EXCEPTION CHECKS
1047 
1048  idx total_rows = 0, total_cols = 0;
1049  for (idx i = 0; i < As.size(); ++i)
1050  {
1051  total_rows += static_cast<idx>(As[i].rows());
1052  total_cols += static_cast<idx>(As[i].cols());
1053  }
1055  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1056 
1057  idx cur_row = 0, cur_col = 0;
1058  for (idx i = 0; i < As.size(); ++i)
1059  {
1060  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1061  cur_row += static_cast<idx>(As[i].rows());
1062  cur_col += static_cast<idx>(As[i].cols());
1063  }
1064 
1065  return result;
1066 }
1067 
1068 // Direct sum of a list of matrices, preserve return type
1069 // deduce the template parameters from initializer_list
1080 template<typename Derived>
1082  const std::initializer_list<Derived>& As)
1083 {
1084  return dirsum(std::vector<Derived>(As));
1085 }
1086 
1096 template<typename Derived>
1098  const Eigen::MatrixBase<Derived>& A, idx n)
1099 {
1100  const dyn_mat<typename Derived::Scalar>& rA = A;
1101 
1102  // EXCEPTION CHECKS
1103 
1104  // check zero-size
1106  throw Exception("qpp::dirsumpow()", Exception::Type::ZERO_SIZE);
1107 
1108  // check out of range
1109  if (n == 0)
1110  throw Exception("qpp::dirsumpow()", Exception::Type::OUT_OF_RANGE);
1111  // END EXCEPTION CHECKS
1112 
1113  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1114 
1115  return dirsum(As);
1116 }
1117 
1129 template<typename Derived>
1130 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1131  idx rows, idx cols)
1132 {
1133  const dyn_mat<typename Derived::Scalar>& rA = A;
1134 
1135  idx Arows = static_cast<idx>(rA.rows());
1136  idx Acols = static_cast<idx>(rA.cols());
1137 
1138  // EXCEPTION CHECKS
1139 
1140  // check zero-size
1142  throw Exception("qpp::reshape()", Exception::Type::ZERO_SIZE);
1143 
1144  if (Arows * Acols != rows * cols)
1145  throw Exception("qpp::reshape()",
1147  // END EXCEPTION CHECKS
1148 
1149  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1150  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1151 }
1152 
1165 template<typename Derived1, typename Derived2>
1166 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1167  const Eigen::MatrixBase<Derived2>& B)
1168 {
1169  const dyn_mat<typename Derived1::Scalar>& rA = A;
1170  const dyn_mat<typename Derived2::Scalar>& rB = B;
1171 
1172  // EXCEPTION CHECKS
1173 
1174  // check types
1175  if (!std::is_same<typename Derived1::Scalar,
1176  typename Derived2::Scalar>::value)
1177  throw Exception("qpp::comm()", Exception::Type::TYPE_MISMATCH);
1178 
1179  // check zero-size
1182  throw Exception("qpp::comm()", Exception::Type::ZERO_SIZE);
1183 
1184  // check square matrices
1186  throw Exception("qpp::comm()", Exception::Type::MATRIX_NOT_SQUARE);
1187 
1188  // check equal dimensions
1189  if (rA.rows() != rB.rows())
1190  throw Exception("qpp::comm()", Exception::Type::DIMS_NOT_EQUAL);
1191  // END EXCEPTION CHECKS
1192 
1193  return rA * rB - rB * rA;
1194 }
1195 
1208 template<typename Derived1, typename Derived2>
1210  const Eigen::MatrixBase<Derived1>& A,
1211  const Eigen::MatrixBase<Derived2>& B)
1212 {
1213  const dyn_mat<typename Derived1::Scalar>& rA = A;
1214  const dyn_mat<typename Derived2::Scalar>& rB = B;
1215 
1216  // EXCEPTION CHECKS
1217 
1218  // check types
1219  if (!std::is_same<typename Derived1::Scalar,
1220  typename Derived2::Scalar>::value)
1221  throw Exception("qpp::anticomm()", Exception::Type::TYPE_MISMATCH);
1222 
1223  // check zero-size
1226  throw Exception("qpp::anticomm()", Exception::Type::ZERO_SIZE);
1227 
1228  // check square matrices
1230  throw Exception("qpp::anticomm()",
1232 
1233  // check equal dimensions
1234  if (rA.rows() != rB.rows())
1235  throw Exception("qpp::anticomm()", Exception::Type::DIMS_NOT_EQUAL);
1236  // END EXCEPTION CHECKS
1237 
1238  return rA * rB + rB * rA;
1239 }
1240 
1251 template<typename Derived>
1252 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& V)
1253 {
1254  const dyn_mat<typename Derived::Scalar>& rV = V;
1255 
1256  // EXCEPTION CHECKS
1257 
1258  // check zero-size
1260  throw Exception("qpp::prj()", Exception::Type::ZERO_SIZE);
1261 
1262  // check column vector
1263  if (!internal::_check_cvector(rV))
1264  throw Exception("qpp::prj()", Exception::Type::MATRIX_NOT_CVECTOR);
1265  // END EXCEPTION CHECKS
1266 
1267  double normV = norm(rV);
1268  if (normV > eps)
1269  return rV * adjoint(rV) / (normV * normV);
1270  else
1272  ::Zero(rV.rows(), rV.rows());
1273 
1274 }
1275 
1283 template<typename Derived>
1284 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& Vs)
1285 {
1286  // EXCEPTION CHECKS
1287 
1288  // check empty list
1290  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1291 
1292  for (auto&& it : Vs)
1294  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1295 
1296  // check that Vs[0] is a column vector
1297  if (!internal::_check_cvector(Vs[0]))
1298  throw Exception("qpp::grams()", Exception::Type::MATRIX_NOT_CVECTOR);
1299 
1300  // now check that all the rest match the size of the first vector
1301  for (auto&& it : Vs)
1302  if (it.rows() != Vs[0].rows() || it.cols() != 1)
1303  throw Exception("qpp::grams()", Exception::Type::DIMS_NOT_EQUAL);
1304  // END EXCEPTION CHECKS
1305 
1308  Vs[0].rows());
1309 
1311  dyn_mat<typename Derived::Scalar>::Zero(Vs[0].rows(), 1);
1312 
1313  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1314  // find the first non-zero vector in the list
1315  idx pos = 0;
1316  for (pos = 0; pos < Vs.size(); ++pos)
1317  {
1318  if (norm(Vs[pos]) > eps) // add it as the first element
1319  {
1320  outvecs.push_back(Vs[pos]);
1321  break;
1322  }
1323  }
1324 
1325  // start the process
1326  for (idx i = pos + 1; i < Vs.size(); ++i)
1327  {
1328  cut -= prj(outvecs[i - 1 - pos]);
1329  vi = cut * Vs[i];
1330  outvecs.push_back(vi);
1331  }
1332 
1333  dyn_mat<typename Derived::Scalar> result(Vs[0].rows(), outvecs.size());
1334 
1335  idx cnt = 0;
1336  for (auto&& it : outvecs)
1337  {
1338  double normV = norm(it);
1339  if (normV > eps) // we add only the non-zero vectors
1340  {
1341  result.col(cnt) = it / normV;
1342  cnt++;
1343  }
1344  }
1345 
1346  return result.block(0, 0, Vs[0].rows(), cnt);
1347 }
1348 
1349 // deduce the template parameters from initializer_list
1357 template<typename Derived>
1359  const std::initializer_list<Derived>& Vs)
1360 {
1361  return grams(std::vector<Derived>(Vs));
1362 }
1363 
1371 template<typename Derived>
1372 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A)
1373 {
1374  const dyn_mat<typename Derived::Scalar>& rA = A;
1375 
1376  // EXCEPTION CHECKS
1377 
1379  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1380  // END EXCEPTION CHECKS
1381 
1382  std::vector<dyn_mat<typename Derived::Scalar>> input;
1383 
1384  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1385  input.push_back(rA.col(i));
1386 
1387  return grams<dyn_mat<typename Derived::Scalar>>(input);
1388 }
1389 
1400 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims)
1401 {
1402  // EXCEPTION CHECKS
1403 
1404  if (!internal::_check_dims(dims))
1405  throw Exception("qpp::n2multiidx()", Exception::Type::DIMS_INVALID);
1406 
1407  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1408  static_cast<idx>(1), std::multiplies<idx>()))
1409  throw Exception("qpp::n2multiidx()", Exception::Type::OUT_OF_RANGE);
1410  // END EXCEPTION CHECKS
1411 
1412  // double the size for matrices reshaped as vectors
1413  idx result[2 * maxn];
1414  internal::_n2multiidx(n, dims.size(), dims.data(), result);
1415 
1416  return std::vector<idx>(result, result + dims.size());
1417 }
1418 
1429 inline idx multiidx2n(const std::vector<idx>& midx,
1430  const std::vector<idx>& dims)
1431 {
1432  // EXCEPTION CHECKS
1433 
1434  if (!internal::_check_dims(dims))
1435  throw Exception("qpp::multiidx2n()", Exception::Type::DIMS_INVALID);
1436 
1437  for (idx i = 0; i < dims.size(); ++i)
1438  if (midx[i] >= dims[i])
1439  throw Exception("qpp::multiidx2n()",
1441  // END EXCEPTION CHECKS
1442 
1443  return internal::_multiidx2n(midx.data(), dims.size(), dims.data());
1444 }
1445 
1458 inline ket mket(const std::vector<idx>& mask,
1459  const std::vector<idx>& dims)
1460 {
1461  idx n = mask.size();
1462 
1463  idx D = std::accumulate(std::begin(dims), std::end(dims),
1464  static_cast<idx>(1), std::multiplies<idx>());
1465 
1466  // EXCEPTION CHECKS
1467 
1468  // check zero size
1469  if (n == 0)
1470  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1471  // check valid dims
1472  if (!internal::_check_dims(dims))
1473  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1474  // check mask and dims have the same size
1475  if (mask.size() != dims.size())
1476  throw Exception("qpp::mket()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1477  // check mask is a valid vector
1478  for (idx i = 0; i < n; ++i)
1479  if (mask[i] >= dims[i])
1480  throw Exception("qpp::mket()",
1482  // END EXCEPTION CHECKS
1483 
1484  ket result = ket::Zero(D);
1485  idx pos = multiidx2n(mask, dims);
1486  result(pos) = 1;
1487 
1488  return result;
1489 }
1490 
1503 inline ket mket(const std::vector<idx>& mask, idx d = 2)
1504 {
1505  idx n = mask.size();
1506  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1507 
1508  // EXCEPTION CHECKS
1509 
1510  // check zero size
1511  if (n == 0)
1512  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1513 
1514  // check valid dims
1515  if (d == 0)
1516  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1517 
1518  // check mask is a valid vector
1519  for (idx i = 0; i < n; ++i)
1520  if (mask[i] >= d)
1521  throw Exception("qpp::mket()",
1523  // END EXCEPTION CHECKS
1524 
1525  ket result = ket::Zero(D);
1526  std::vector<idx> dims(n, d);
1527  idx pos = multiidx2n(mask, dims);
1528  result(pos) = 1;
1529 
1530  return result;
1531 }
1532 
1547 inline cmat mprj(const std::vector<idx>& mask,
1548  const std::vector<idx>& dims)
1549 {
1550  idx n = mask.size();
1551 
1552  idx D = std::accumulate(std::begin(dims), std::end(dims),
1553  static_cast<idx>(1), std::multiplies<idx>());
1554 
1555  // EXCEPTION CHECKS
1556 
1557  // check zero size
1558  if (n == 0)
1559  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1560  // check valid dims
1561  if (!internal::_check_dims(dims))
1562  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1563  // check mask and dims have the same size
1564  if (mask.size() != dims.size())
1565  throw Exception("qpp::mprj()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1566  // check mask is a valid vector
1567  for (idx i = 0; i < n; ++i)
1568  if (mask[i] >= dims[i])
1569  throw Exception("qpp::mprj()",
1571  // END EXCEPTION CHECKS
1572 
1573  cmat result = cmat::Zero(D, D);
1574  idx pos = multiidx2n(mask, dims);
1575  result(pos, pos) = 1;
1576 
1577  return result;
1578 }
1579 
1594 inline cmat mprj(const std::vector<idx>& mask, idx d = 2)
1595 {
1596  idx n = mask.size();
1597  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1598 
1599  // EXCEPTION CHECKS
1600 
1601  // check zero size
1602  if (n == 0)
1603  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1604 
1605  // check valid dims
1606  if (d == 0)
1607  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1608 
1609  // check mask is a valid vector
1610  for (idx i = 0; i < n; ++i)
1611  if (mask[i] >= d)
1612  throw Exception("qpp::mprj()",
1614  // END EXCEPTION CHECKS
1615 
1616  cmat result = cmat::Zero(D, D);
1617  std::vector<idx> dims(n, d);
1618  idx pos = multiidx2n(mask, dims);
1619  result(pos, pos) = 1;
1620 
1621  return result;
1622 }
1623 
1632 template<typename InputIterator>
1633 std::vector<double> abssq(InputIterator first, InputIterator last)
1634 {
1635  std::vector<double> weights(std::distance(first, last));
1636  std::transform(first, last, std::begin(weights),
1637  [](cplx z) -> double
1638  {
1639  return std::norm(z);
1640  });
1641 
1642  return weights;
1643 }
1644 
1651 template<typename Container>
1652 std::vector<double> abssq(const Container& c,
1653  // need the enable_if to SFINAE out Eigen expressions
1654  // that will otherwise match, instead of matching
1655  // the overload below:
1656 
1657  // template<typename Derived>
1658  // abssq(const Eigen::MatrixBase<Derived>& V)
1659  typename std::enable_if<
1661  >::type* = nullptr)
1662 {
1663  return abssq(std::begin(c), std::end(c));
1664 }
1665 
1672 template<typename Derived>
1673 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A)
1674 {
1675  const dyn_mat<typename Derived::Scalar>& rA = A;
1676 
1677  // EXCEPTION CHECKS
1678 
1679  // check zero-size
1681  throw Exception("qpp::abssq()", Exception::Type::ZERO_SIZE);
1682  // END EXCEPTION CHECKS
1683 
1684  return abssq(rA.data(), rA.data() + rA.size());
1685 }
1686 
1695 template<typename InputIterator>
1696 typename std::iterator_traits<InputIterator>::value_type
1697 sum(InputIterator first, InputIterator last)
1698 {
1699  using value_type =
1700  typename std::iterator_traits<InputIterator>::value_type;
1701 
1702  return std::accumulate(first, last, static_cast<value_type>(0));
1703 }
1704 
1712 template<typename Container>
1713 typename Container::value_type
1714 sum(const Container& c)
1715 {
1716  return sum(std::begin(c), std::end(c));
1717 }
1718 
1727 template<typename InputIterator>
1728 typename std::iterator_traits<InputIterator>::value_type
1729 prod(InputIterator first, InputIterator last)
1730 {
1731  using value_type =
1732  typename std::iterator_traits<InputIterator>::value_type;
1733 
1734  return std::accumulate(first, last, static_cast<value_type>(1),
1735  std::multiplies<value_type>());
1736 }
1737 
1745 template<typename Container>
1746 typename Container::value_type
1747 prod(const Container& c)
1748 {
1749  return prod(std::begin(c), std::end(c));
1750 }
1751 
1764 template<typename Derived>
1766  const Eigen::MatrixBase<Derived>& A)
1767 {
1768  const dyn_mat<typename Derived::Scalar>& rA = A;
1769 
1770  // EXCEPTION CHECKS
1771  // check zero-size
1773  throw Exception("qpp::rho2pure()", Exception::Type::ZERO_SIZE);
1774  // check square matrix
1775  if (!internal::_check_square_mat(rA))
1776  throw Exception("qpp::rho2pure()", Exception::Type::MATRIX_NOT_SQUARE);
1777  // END EXCEPTION CHECKS
1778 
1779  dyn_col_vect<double> tmp_evals = hevals(rA);
1780  cmat tmp_evects = hevects(rA);
1783  // find the non-zero eigenvector
1784  // there is only one, assuming the state is pure
1785  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1786  {
1787  if (std::abs(tmp_evals(k)) > eps)
1788  {
1789  result = tmp_evects.col(k);
1790  break;
1791  }
1792  }
1793 
1794  return result;
1795 }
1796 
1805 template<typename T>
1806 std::vector<T> complement(std::vector<T> subsys, idx N)
1807 {
1808  // EXCEPTION CHECKS
1809 
1810  if (N < subsys.size())
1811  throw Exception("qpp::complement()", Exception::Type::OUT_OF_RANGE);
1812  // END EXCEPTION CHECKS
1813 
1814  std::vector<T> all(N);
1815  std::vector<T> subsys_bar(N - subsys.size());
1816 
1817  std::iota(std::begin(all), std::end(all), 0);
1818  std::sort(std::begin(subsys), std::end(subsys));
1819  std::set_difference(std::begin(all), std::end(all),
1820  std::begin(subsys), std::end(subsys),
1821  std::begin(subsys_bar));
1822 
1823  return subsys_bar;
1824 }
1825 
1836 template<typename Derived>
1837 std::vector<double> rho2bloch(
1838  const Eigen::MatrixBase<Derived>& A)
1839 {
1840  const dyn_mat<typename Derived::Scalar>& rA = A;
1841 
1842  // EXCEPTION CHECKS
1843 
1844  // check qubit matrix
1846  throw Exception("qpp::rho2bloch()", Exception::Type::NOT_QUBIT_MATRIX);
1847  // END EXCEPTION CHECKS
1848 
1849  std::vector<double> result(3);
1850  cmat X(2, 2), Y(2, 2), Z(2, 2);
1851  X << 0, 1, 1, 0;
1852  Y << 0, -1_i, 1_i, 0;
1853  Z << 1, 0, 0, -1;
1854  result[0] = std::real(trace(rA * X));
1855  result[1] = std::real(trace(rA * Y));
1856  result[2] = std::real(trace(rA * Z));
1857 
1858  return result;
1859 }
1860 
1869 inline cmat bloch2rho(const std::vector<double>& r)
1870 {
1871  // EXCEPTION CHECKS
1872 
1873  // check 3-dimensional vector
1874  if (r.size() != 3)
1875  throw Exception("qpp::bloch2rho", "r is not a 3-dimensional vector!");
1876  // END EXCEPTION CHECKS
1877 
1878  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1879  X << 0, 1, 1, 0;
1880  Y << 0, -1_i, 1_i, 0;
1881  Z << 1, 0, 0, -1;
1882  Id2 << 1, 0, 0, 1;
1883 
1884  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1885 }
1886 
1887 } /* namespace qpp */
1888 
1889 #endif /* FUNCTIONS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:173
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:80
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolut value.
Definition: functions.h:609
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:272
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:634
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:73
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:1005
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:83
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:1837
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:56
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1209
Quantum++ main namespace.
Definition: codes.h:30
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:739
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:391
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:853
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:1765
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1097
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1806
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:332
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:551
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:823
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:499
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1869
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:278
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1400
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:1429
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:475
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:448
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:890
bool _check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:220
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:51
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:54
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1547
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:971
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Matrix power.
Definition: functions.h:783
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &V)
Projector.
Definition: functions.h:1252
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:306
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:417
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:64
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:94
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:313
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:524
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:709
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:95
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:684
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of an STL-like range of complex numbers.
Definition: functions.h:1633
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1166
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:66
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1458
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:41
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:659
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:584
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:126
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:209
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1130
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &Vs)
Gram-Schmidt orthogonalization.
Definition: functions.h:1284
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:363