Quantum++  v0.8.8.2
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 {
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
184  if ( !internal::_check_square_mat(rA))
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>
277 
278 eig(const Eigen::MatrixBase<Derived>& A)
279 {
280  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
281 
282  // EXCEPTION CHECKS
283 
284  // check zero-size
286  throw Exception("qpp::eig()", Exception::Type::ZERO_SIZE);
287 
288  // check square matrix
289  if ( !internal::_check_square_mat(rA))
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.derived();
309 
310  // EXCEPTION CHECKS
311 
312  // check zero-size
314  throw Exception("qpp::evals()", Exception::Type::ZERO_SIZE);
315 
316  // check square matrix
317  if ( !internal::_check_square_mat(rA))
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.derived();
335 
336  // EXCEPTION CHECKS
337 
338  // check zero-size
340  throw Exception("qpp::evects()", Exception::Type::ZERO_SIZE);
341 
342  // check square matrix
343  if ( !internal::_check_square_mat(rA))
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.derived();
366 
367  // EXCEPTION CHECKS
368 
369  // check zero-size
371  throw Exception("qpp::heig()", Exception::Type::ZERO_SIZE);
372 
373  // check square matrix
374  if ( !internal::_check_square_mat(rA))
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.derived();
394 
395  // EXCEPTION CHECKS
396 
397  // check zero-size
399  throw Exception("qpp::hevals()", Exception::Type::ZERO_SIZE);
400 
401  // check square matrix
402  if ( !internal::_check_square_mat(rA))
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.derived();
420 
421  // EXCEPTION CHECKS
422 
423  // check zero-size
425  throw Exception("qpp::hevects()", Exception::Type::ZERO_SIZE);
426 
427  // check square matrix
428  if ( !internal::_check_square_mat(rA))
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.derived();
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.derived();
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.derived();
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.derived();
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.derived();
554 
555  // EXCEPTION CHECKS
556 
557  // check zero-size
559  throw Exception("qpp::funm()", Exception::Type::ZERO_SIZE);
560 
561  // check square matrix
562  if ( !internal::_check_square_mat(rA))
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.derived();
587 
588  // EXCEPTION CHECKS
589 
590  // check zero-size
592  throw Exception("qpp::sqrtm()", Exception::Type::ZERO_SIZE);
593 
594  // check square matrix
595  if ( !internal::_check_square_mat(rA))
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.derived();
612 
613  // EXCEPTION CHECKS
614 
615  // check zero-size
617  throw Exception("qpp::absm()", Exception::Type::ZERO_SIZE);
618 
619  // check square matrix
620  if ( !internal::_check_square_mat(rA))
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.derived();
637 
638  // EXCEPTION CHECKS
639 
640  // check zero-size
642  throw Exception("qpp::expm()", Exception::Type::ZERO_SIZE);
643 
644  // check square matrix
645  if ( !internal::_check_square_mat(rA))
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.derived();
662 
663  // EXCEPTION CHECKS
664 
665  // check zero-size
667  throw Exception("qpp::logm()", Exception::Type::ZERO_SIZE);
668 
669  // check square matrix
670  if ( !internal::_check_square_mat(rA))
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.derived();
687 
688  // EXCEPTION CHECKS
689 
690  // check zero-size
692  throw Exception("qpp::sinm()", Exception::Type::ZERO_SIZE);
693 
694  // check square matrix
695  if ( !internal::_check_square_mat(rA))
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.derived();
712 
713  // EXCEPTION CHECKS
714 
715  // check zero-size
717  throw Exception("qpp::cosm()", Exception::Type::ZERO_SIZE);
718 
719  // check square matrix
720  if ( !internal::_check_square_mat(rA))
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.derived();
742 
743  // EXCEPTION CHECKS
744 
745  // check zero-size
747  throw Exception("qpp::spectralpowm()", Exception::Type::ZERO_SIZE);
748 
749  // check square matrix
750  if ( !internal::_check_square_mat(rA))
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 {
786  // EXCEPTION CHECKS
787 
788  // check zero-size
790  throw Exception("qpp::powm()", Exception::Type::ZERO_SIZE);
791 
792  // check square matrix
794  throw Exception("qpp::powm()", Exception::Type::MATRIX_NOT_SQUARE);
795  // END EXCEPTION CHECKS
796 
797  // if n = 1, return the matrix unchanged
798  if ( n == 1 )
799  return A;
800 
803  A.rows(), A.rows());
804 
805  // if n = 0, return the identity (as just prepared in result)
806  if ( n == 0 )
807  return result;
808 
809  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
810 
811  // fast matrix power
812  for ( ; n > 0; n /= 2 )
813  {
814  if ( n % 2 )
815  result = (result * cA).eval();
816  cA = (cA * cA).eval();
817  }
818 
819  return result;
820 }
821 
830 template<typename Derived>
831 double schatten(const Eigen::MatrixBase<Derived>& A, double p)
832 {
833  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
834 
835  // EXCEPTION CHECKS
836 
837  // check zero-size
839  throw Exception("qpp::schatten()", Exception::Type::ZERO_SIZE);
840  if ( p < 1 )
841  throw Exception("qpp::schatten()", Exception::Type::OUT_OF_RANGE);
842  // END EXCEPTION CHECKS
843 
844  if ( p == infty ) // infinity norm (largest singular value)
845  return svals(rA)(0);
846 
847  return std::pow(trace(powm(absm(rA), p)).real(), 1. / p);
848 }
849 
850 // other functions
851 
860 template<typename OutputScalar, typename Derived>
861 dyn_mat <OutputScalar> cwise(const Eigen::MatrixBase<Derived>& A,
862  OutputScalar (* f)(
863  const typename Derived::Scalar&))
864 {
865  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
866 
867  // EXCEPTION CHECKS
868 
869  // check zero-size
871  throw Exception("qpp::cwise()", Exception::Type::ZERO_SIZE);
872  // END EXCEPTION CHECKS
873 
874  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
875 
876 #ifdef _WITH_OPENMP_
877 #pragma omp parallel for collapse(2)
878 #endif
879  // column major order for speed
880  for ( idx j = 0; j < static_cast<idx>(rA.cols()); ++j )
881  for ( idx i = 0; i < static_cast<idx>(rA.rows()); ++i )
882  result(i, j) = (*f)(rA(i, j));
883 
884  return result;
885 }
886 
887 // Kronecker product of multiple matrices, preserve return type
888 // variadic template
899 template<typename T>
901 {
902  return head;
903 }
904 
915 template<typename T, typename ... Args>
916 dyn_mat<typename T::Scalar> kron(const T& head, const Args& ... tail)
917 {
918  return internal::_kron2(head, kron(tail...));
919 }
920 
930 template<typename Derived>
931 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As)
932 {
933  // EXCEPTION CHECKS
934 
935  if ( As.size() == 0 )
936  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
937 
938  for ( auto&& it : As )
940  throw Exception("qpp::kron()", Exception::Type::ZERO_SIZE);
941  // END EXCEPTION CHECKS
942 
943  dyn_mat<typename Derived::Scalar> result = As[0].derived();
944  for ( idx i = 1; i < As.size(); ++i )
945  {
946  result = kron(result, As[i]);
947  }
948 
949  return result;
950 }
951 
952 // Kronecker product of a list of matrices, preserve return type
953 // deduce the template parameters from initializer_list
964 template<typename Derived>
966  const std::initializer_list<Derived>& As)
967 {
968  return kron(std::vector<Derived>(As));
969 }
970 
980 template<typename Derived>
981 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
982  idx n)
983 {
984  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
985 
986  // EXCEPTION CHECKS
987 
988  // check zero-size
990  throw Exception("qpp::kronpow()", Exception::Type::ZERO_SIZE);
991 
992  // check out of range
993  if ( n == 0 )
994  throw Exception("qpp::kronpow()", Exception::Type::OUT_OF_RANGE);
995  // END EXCEPTION CHECKS
996 
997  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
998 
999  return kron(As);
1000 }
1001 
1002 // Direct sum of multiple matrices, preserve return type
1003 // variadic template
1014 template<typename T>
1016 {
1017  return head;
1018 }
1019 
1030 template<typename T, typename ... Args>
1031 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args& ... tail)
1032 {
1033  return internal::_dirsum2(head, dirsum(tail...));
1034 }
1035 
1045 template<typename Derived>
1046 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As)
1047 {
1048  // EXCEPTION CHECKS
1049 
1050  if ( As.size() == 0 )
1051  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1052 
1053  for ( auto&& it : As )
1055  throw Exception("qpp::dirsum()", Exception::Type::ZERO_SIZE);
1056  // END EXCEPTION CHECKS
1057 
1058  idx total_rows = 0, total_cols = 0;
1059  for ( idx i = 0; i < As.size(); ++i )
1060  {
1061  total_rows += static_cast<idx>(As[i].rows());
1062  total_cols += static_cast<idx>(As[i].cols());
1063  }
1065  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1066 
1067  idx cur_row = 0, cur_col = 0;
1068  for ( idx i = 0; i < As.size(); ++i )
1069  {
1070  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1071  cur_row += static_cast<idx>(As[i].rows());
1072  cur_col += static_cast<idx>(As[i].cols());
1073  }
1074 
1075  return result;
1076 }
1077 
1078 // Direct sum of a list of matrices, preserve return type
1079 // deduce the template parameters from initializer_list
1090 template<typename Derived>
1092  const std::initializer_list<Derived>& As)
1093 {
1094  return dirsum(std::vector<Derived>(As));
1095 }
1096 
1106 template<typename Derived>
1108  const Eigen::MatrixBase<Derived>& A, idx n)
1109 {
1110  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1111 
1112  // EXCEPTION CHECKS
1113 
1114  // check zero-size
1116  throw Exception("qpp::dirsumpow()", Exception::Type::ZERO_SIZE);
1117 
1118  // check out of range
1119  if ( n == 0 )
1120  throw Exception("qpp::dirsumpow()", Exception::Type::OUT_OF_RANGE);
1121  // END EXCEPTION CHECKS
1122 
1123  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1124 
1125  return dirsum(As);
1126 }
1127 
1139 template<typename Derived>
1140 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1141  idx rows, idx cols)
1142 {
1143  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1144 
1145  idx Arows = static_cast<idx>(rA.rows());
1146  idx Acols = static_cast<idx>(rA.cols());
1147 
1148  // EXCEPTION CHECKS
1149 
1150  // check zero-size
1152  throw Exception("qpp::reshape()", Exception::Type::ZERO_SIZE);
1153 
1154  if ( Arows * Acols != rows * cols )
1155  throw Exception("qpp::reshape()",
1157  // END EXCEPTION CHECKS
1158 
1159  return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1160  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1161 }
1162 
1175 template<typename Derived1, typename Derived2>
1176 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1177  const Eigen::MatrixBase<Derived2>& B)
1178 {
1179  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1180  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1181 
1182  // EXCEPTION CHECKS
1183 
1184  // check types
1185  if ( !std::is_same<typename Derived1::Scalar,
1186  typename Derived2::Scalar>::value )
1187  throw Exception("qpp::comm()", Exception::Type::TYPE_MISMATCH);
1188 
1189  // check zero-size
1192  throw Exception("qpp::comm()", Exception::Type::ZERO_SIZE);
1193 
1194  // check square matrices
1196  throw Exception("qpp::comm()", Exception::Type::MATRIX_NOT_SQUARE);
1197 
1198  // check equal dimensions
1199  if ( rA.rows() != rB.rows())
1200  throw Exception("qpp::comm()", Exception::Type::DIMS_NOT_EQUAL);
1201  // END EXCEPTION CHECKS
1202 
1203  return rA * rB - rB * rA;
1204 }
1205 
1218 template<typename Derived1, typename Derived2>
1220  const Eigen::MatrixBase<Derived1>& A,
1221  const Eigen::MatrixBase<Derived2>& B)
1222 {
1223  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1224  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1225 
1226  // EXCEPTION CHECKS
1227 
1228  // check types
1229  if ( !std::is_same<typename Derived1::Scalar,
1230  typename Derived2::Scalar>::value )
1231  throw Exception("qpp::anticomm()", Exception::Type::TYPE_MISMATCH);
1232 
1233  // check zero-size
1236  throw Exception("qpp::anticomm()", Exception::Type::ZERO_SIZE);
1237 
1238  // check square matrices
1240  throw Exception("qpp::anticomm()",
1242 
1243  // check equal dimensions
1244  if ( rA.rows() != rB.rows())
1245  throw Exception("qpp::anticomm()", Exception::Type::DIMS_NOT_EQUAL);
1246  // END EXCEPTION CHECKS
1247 
1248  return rA * rB + rB * rA;
1249 }
1250 
1261 template<typename Derived>
1262 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& V)
1263 {
1264  const dyn_mat<typename Derived::Scalar>& rV = V.derived();
1265 
1266  // EXCEPTION CHECKS
1267 
1268  // check zero-size
1270  throw Exception("qpp::prj()", Exception::Type::ZERO_SIZE);
1271 
1272  // check column vector
1273  if ( !internal::_check_cvector(rV))
1274  throw Exception("qpp::prj()", Exception::Type::MATRIX_NOT_CVECTOR);
1275  // END EXCEPTION CHECKS
1276 
1277  double normV = norm(rV);
1278  if ( normV > eps )
1279  return rV * adjoint(rV) / (normV * normV);
1280  else
1282  ::Zero(rV.rows(), rV.rows());
1283 
1284 }
1285 
1293 template<typename Derived>
1294 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& Vs)
1295 {
1296  // EXCEPTION CHECKS
1297 
1298  // check empty list
1300  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1301 
1302  for ( auto&& it : Vs )
1304  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1305 
1306  // check that Vs[0] is a column vector
1307  if ( !internal::_check_cvector(Vs[0]))
1308  throw Exception("qpp::grams()", Exception::Type::MATRIX_NOT_CVECTOR);
1309 
1310  // now check that all the rest match the size of the first vector
1311  for ( auto&& it : Vs )
1312  if ( it.rows() != Vs[0].rows() || it.cols() != 1 )
1313  throw Exception("qpp::grams()", Exception::Type::DIMS_NOT_EQUAL);
1314  // END EXCEPTION CHECKS
1315 
1318  Vs[0].rows());
1319 
1321  dyn_mat<typename Derived::Scalar>::Zero(Vs[0].rows(), 1);
1322 
1323  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1324  // find the first non-zero vector in the list
1325  idx pos = 0;
1326  for ( pos = 0; pos < Vs.size(); ++pos )
1327  {
1328  if ( norm(Vs[pos]) > eps ) // add it as the first element
1329  {
1330  outvecs.push_back(Vs[pos]);
1331  break;
1332  }
1333  }
1334 
1335  // start the process
1336  for ( idx i = pos + 1; i < Vs.size(); ++i )
1337  {
1338  cut -= prj(outvecs[i - 1 - pos]);
1339  vi = cut * Vs[i];
1340  outvecs.push_back(vi);
1341  }
1342 
1343  dyn_mat<typename Derived::Scalar> result(Vs[0].rows(), outvecs.size());
1344 
1345  idx cnt = 0;
1346  for ( auto&& it : outvecs )
1347  {
1348  double normV = norm(it);
1349  if ( normV > eps ) // we add only the non-zero vectors
1350  {
1351  result.col(cnt) = it / normV;
1352  cnt++;
1353  }
1354  }
1355 
1356  return result.block(0, 0, Vs[0].rows(), cnt);
1357 }
1358 
1359 // deduce the template parameters from initializer_list
1367 template<typename Derived>
1369  const std::initializer_list<Derived>& Vs)
1370 {
1371  return grams(std::vector<Derived>(Vs));
1372 }
1373 
1381 template<typename Derived>
1382 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A)
1383 {
1384  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1385 
1386  // EXCEPTION CHECKS
1387 
1389  throw Exception("qpp::grams()", Exception::Type::ZERO_SIZE);
1390  // END EXCEPTION CHECKS
1391 
1392  std::vector<dyn_mat<typename Derived::Scalar>> input;
1393 
1394  for ( idx i = 0; i < static_cast<idx>(rA.cols()); ++i )
1395  input.push_back(rA.col(i));
1396 
1397  return grams<dyn_mat<typename Derived::Scalar>>(input);
1398 }
1399 
1410 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims)
1411 {
1412  // EXCEPTION CHECKS
1413 
1414  if ( !internal::_check_dims(dims))
1415  throw Exception("qpp::n2multiidx()", Exception::Type::DIMS_INVALID);
1416 
1417  if ( n >= std::accumulate(std::begin(dims), std::end(dims),
1418  static_cast<idx>(1), std::multiplies<idx>()))
1419  throw Exception("qpp::n2multiidx()", Exception::Type::OUT_OF_RANGE);
1420  // END EXCEPTION CHECKS
1421 
1422  // double the size for matrices reshaped as vectors
1423  idx result[2 * maxn];
1424  internal::_n2multiidx(n, dims.size(), dims.data(), result);
1425 
1426  return std::vector<idx>(result, result + dims.size());
1427 }
1428 
1439 inline idx multiidx2n(const std::vector<idx>& midx,
1440  const std::vector<idx>& dims)
1441 {
1442  // EXCEPTION CHECKS
1443 
1444  if ( !internal::_check_dims(dims))
1445  throw Exception("qpp::multiidx2n()", Exception::Type::DIMS_INVALID);
1446 
1447  for ( idx i = 0; i < dims.size(); ++i )
1448  if ( midx[i] >= dims[i] )
1449  throw Exception("qpp::multiidx2n()",
1451  // END EXCEPTION CHECKS
1452 
1453  return internal::_multiidx2n(midx.data(), dims.size(), dims.data());
1454 }
1455 
1468 inline ket mket(const std::vector<idx>& mask,
1469  const std::vector<idx>& dims)
1470 {
1471  idx n = mask.size();
1472 
1473  idx D = std::accumulate(std::begin(dims), std::end(dims),
1474  static_cast<idx>(1), std::multiplies<idx>());
1475 
1476  // EXCEPTION CHECKS
1477 
1478  // check zero size
1479  if ( n == 0 )
1480  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1481  // check valid dims
1482  if ( !internal::_check_dims(dims))
1483  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1484  // check mask and dims have the same size
1485  if ( mask.size() != dims.size())
1486  throw Exception("qpp::mket()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1487  // check mask is a valid vector
1488  for ( idx i = 0; i < n; ++i )
1489  if ( mask[i] >= dims[i] )
1490  throw Exception("qpp::mket()",
1492  // END EXCEPTION CHECKS
1493 
1494  ket result = ket::Zero(D);
1495  idx pos = multiidx2n(mask, dims);
1496  result(pos) = 1;
1497 
1498  return result;
1499 }
1500 
1513 inline ket mket(const std::vector<idx>& mask, idx d = 2)
1514 {
1515  idx n = mask.size();
1516  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1517 
1518  // EXCEPTION CHECKS
1519 
1520  // check zero size
1521  if ( n == 0 )
1522  throw Exception("qpp::mket()", Exception::Type::ZERO_SIZE);
1523 
1524  // check valid dims
1525  if ( d == 0 )
1526  throw Exception("qpp::mket()", Exception::Type::DIMS_INVALID);
1527 
1528  // check mask is a valid vector
1529  for ( idx i = 0; i < n; ++i )
1530  if ( mask[i] >= d )
1531  throw Exception("qpp::mket()",
1533  // END EXCEPTION CHECKS
1534 
1535  ket result = ket::Zero(D);
1536  std::vector<idx> dims(n, d);
1537  idx pos = multiidx2n(mask, dims);
1538  result(pos) = 1;
1539 
1540  return result;
1541 }
1542 
1557 inline cmat mprj(const std::vector<idx>& mask,
1558  const std::vector<idx>& dims)
1559 {
1560  idx n = mask.size();
1561 
1562  idx D = std::accumulate(std::begin(dims), std::end(dims),
1563  static_cast<idx>(1), std::multiplies<idx>());
1564 
1565  // EXCEPTION CHECKS
1566 
1567  // check zero size
1568  if ( n == 0 )
1569  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1570  // check valid dims
1571  if ( !internal::_check_dims(dims))
1572  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1573  // check mask and dims have the same size
1574  if ( mask.size() != dims.size())
1575  throw Exception("qpp::mprj()", Exception::Type::SUBSYS_MISMATCH_DIMS);
1576  // check mask is a valid vector
1577  for ( idx i = 0; i < n; ++i )
1578  if ( mask[i] >= dims[i] )
1579  throw Exception("qpp::mprj()",
1581  // END EXCEPTION CHECKS
1582 
1583  cmat result = cmat::Zero(D, D);
1584  idx pos = multiidx2n(mask, dims);
1585  result(pos, pos) = 1;
1586 
1587  return result;
1588 }
1589 
1604 inline cmat mprj(const std::vector<idx>& mask, idx d = 2)
1605 {
1606  idx n = mask.size();
1607  idx D = static_cast<idx>(std::llround(std::pow(d, n)));
1608 
1609  // EXCEPTION CHECKS
1610 
1611  // check zero size
1612  if ( n == 0 )
1613  throw Exception("qpp::mprj()", Exception::Type::ZERO_SIZE);
1614 
1615  // check valid dims
1616  if ( d == 0 )
1617  throw Exception("qpp::mprj()", Exception::Type::DIMS_INVALID);
1618 
1619  // check mask is a valid vector
1620  for ( idx i = 0; i < n; ++i )
1621  if ( mask[i] >= d )
1622  throw Exception("qpp::mprj()",
1624  // END EXCEPTION CHECKS
1625 
1626  cmat result = cmat::Zero(D, D);
1627  std::vector<idx> dims(n, d);
1628  idx pos = multiidx2n(mask, dims);
1629  result(pos, pos) = 1;
1630 
1631  return result;
1632 }
1633 
1642 template<typename InputIterator>
1643 std::vector<double> abssq(InputIterator first, InputIterator last)
1644 {
1645  std::vector<double> weights(std::distance(first, last));
1646  std::transform(first, last, std::begin(weights),
1647  [](cplx z) -> double
1648  {
1649  return std::norm(z);
1650  });
1651 
1652  return weights;
1653 }
1654 
1661 template<typename Container>
1662 std::vector<double> abssq(const Container& c,
1663  typename std::enable_if<
1665  >::type* = nullptr)
1666 // we need the std::enable_if to SFINAE out Eigen expressions
1667 // that will otherwise match, instead of matching
1668 // the overload below:
1669 
1670 // template<typename Derived>
1671 // abssq(const Eigen::MatrixBase<Derived>& V)
1672 {
1673  return abssq(std::begin(c), std::end(c));
1674 }
1675 
1682 template<typename Derived>
1683 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A)
1684 {
1685  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1686 
1687  // EXCEPTION CHECKS
1688 
1689  // check zero-size
1691  throw Exception("qpp::abssq()", Exception::Type::ZERO_SIZE);
1692  // END EXCEPTION CHECKS
1693 
1694  return abssq(rA.data(), rA.data() + rA.size());
1695 }
1696 
1705 template<typename InputIterator>
1706 typename std::iterator_traits<InputIterator>::value_type
1707 sum(InputIterator first, InputIterator last)
1708 {
1709  using value_type =
1710  typename std::iterator_traits<InputIterator>::value_type;
1711 
1712  return std::accumulate(first, last, static_cast<value_type>(0));
1713 }
1714 
1722 template<typename Container>
1723 typename Container::value_type
1724 sum(const Container& c,
1725  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1726 {
1727  return sum(std::begin(c), std::end(c));
1728 }
1729 
1738 template<typename InputIterator>
1739 typename std::iterator_traits<InputIterator>::value_type
1740 prod(InputIterator first, InputIterator last)
1741 {
1742  using value_type =
1743  typename std::iterator_traits<InputIterator>::value_type;
1744 
1745  return std::accumulate(first, last, static_cast<value_type>(1),
1746  std::multiplies<value_type>());
1747 }
1748 
1756 template<typename Container>
1757 typename Container::value_type
1758 prod(const Container& c,
1759  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1760 {
1761  return prod(std::begin(c), std::end(c));
1762 }
1763 
1776 template<typename Derived>
1778  const Eigen::MatrixBase<Derived>& A)
1779 {
1780  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1781 
1782  // EXCEPTION CHECKS
1783  // check zero-size
1785  throw Exception("qpp::rho2pure()", Exception::Type::ZERO_SIZE);
1786  // check square matrix
1787  if ( !internal::_check_square_mat(rA))
1788  throw Exception("qpp::rho2pure()", Exception::Type::MATRIX_NOT_SQUARE);
1789  // END EXCEPTION CHECKS
1790 
1791  dyn_col_vect<double> tmp_evals = hevals(rA);
1792  cmat tmp_evects = hevects(rA);
1795  // find the non-zero eigenvector
1796  // there is only one, assuming the state is pure
1797  for ( idx k = 0; k < static_cast<idx>(rA.rows()); ++k )
1798  {
1799  if ( std::abs(tmp_evals(k)) > eps )
1800  {
1801  result = tmp_evects.col(k);
1802  break;
1803  }
1804  }
1805 
1806  return result;
1807 }
1808 
1817 template<typename T>
1818 std::vector<T> complement(std::vector<T> subsys, idx N)
1819 {
1820  // EXCEPTION CHECKS
1821 
1822  if ( N < subsys.size())
1823  throw Exception("qpp::complement()", Exception::Type::OUT_OF_RANGE);
1824  // END EXCEPTION CHECKS
1825 
1826  std::vector<T> all(N);
1827  std::vector<T> subsys_bar(N - subsys.size());
1828 
1829  std::iota(std::begin(all), std::end(all), 0);
1830  std::sort(std::begin(subsys), std::end(subsys));
1831  std::set_difference(std::begin(all), std::end(all),
1832  std::begin(subsys), std::end(subsys),
1833  std::begin(subsys_bar));
1834 
1835  return subsys_bar;
1836 }
1837 
1848 template<typename Derived>
1849 std::vector<double> rho2bloch(
1850  const Eigen::MatrixBase<Derived>& A)
1851 {
1852  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1853 
1854  // EXCEPTION CHECKS
1855 
1856  // check qubit matrix
1858  throw Exception("qpp::rho2bloch()", Exception::Type::NOT_QUBIT_MATRIX);
1859  // END EXCEPTION CHECKS
1860 
1861  std::vector<double> result(3);
1862  cmat X(2, 2), Y(2, 2), Z(2, 2);
1863  X << 0, 1, 1, 0;
1864  Y << 0, -1_i, 1_i, 0;
1865  Z << 1, 0, 0, -1;
1866  result[0] = std::real(trace(rA * X));
1867  result[1] = std::real(trace(rA * Y));
1868  result[2] = std::real(trace(rA * Z));
1869 
1870  return result;
1871 }
1872 
1881 inline cmat bloch2rho(const std::vector<double>& r)
1882 {
1883  // EXCEPTION CHECKS
1884 
1885  // check 3-dimensional vector
1886  if ( r.size() != 3 )
1887  throw Exception("qpp::bloch2rho", "r is not a 3-dimensional vector!");
1888  // END EXCEPTION CHECKS
1889 
1890  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1891  X << 0, 1, 1, 0;
1892  Y << 0, -1_i, 1_i, 0;
1893  Z << 1, 0, 0, -1;
1894  Id2 << 1, 0, 0, 1;
1895 
1896  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1897 }
1898 
1899 } /* namespace qpp */
1900 
1901 #endif /* FUNCTIONS_H_ */
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:105
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:74
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)( const typename Derived::Scalar &))
Functor.
Definition: functions.h:861
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:257
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:67
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:1015
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:1849
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:1219
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:84
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:391
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:1777
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1107
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1818
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:831
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:1881
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:1410
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:1439
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:900
bool _check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:213
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:61
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1557
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:981
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:1262
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:112
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:88
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:300
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:1643
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1176
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:1468
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:48
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:125
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:209
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1140
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &Vs)
Gram-Schmidt orthogonalization.
Definition: functions.h:1294
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:363