Quantum++  v1.0-rc2
A modern 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 namespace qpp
31 {
32 // Eigen function wrappers
40 template<typename Derived>
42  const Eigen::MatrixBase<Derived>& A)
43 {
44  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
45 
46  // check zero-size
48  throw exception::ZeroSize("qpp::transpose()");
49 
50  return rA.transpose();
51 }
52 
60 template<typename Derived>
62  const Eigen::MatrixBase<Derived>& A)
63 {
64  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
65 
66  // check zero-size
68  throw exception::ZeroSize("qpp::conjugate()");
69 
70  return rA.conjugate();
71 }
72 
80 template<typename Derived>
81 dyn_mat<typename Derived::Scalar> adjoint(const Eigen::MatrixBase<Derived>& A)
82 {
83  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
84 
85  // EXCEPTION CHECKS
86 
87  // check zero-size
89  throw exception::ZeroSize("qpp::adjoint()");
90  // END EXCEPTION CHECKS
91 
92  return rA.adjoint();
93 }
94 
102 template<typename Derived>
103 dyn_mat<typename Derived::Scalar> inverse(const Eigen::MatrixBase<Derived>& A)
104 {
105  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
106 
107  // EXCEPTION CHECKS
108 
109  // check zero-size
111  throw exception::ZeroSize("qpp::inverse()");
112  // END EXCEPTION CHECKS
113 
114  return rA.inverse();
115 }
116 
123 template<typename Derived>
124 typename Derived::Scalar trace(const Eigen::MatrixBase<Derived>& A)
125 {
126  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
127 
128  // EXCEPTION CHECKS
129 
130  // check zero-size
132  throw exception::ZeroSize("qpp::trace()");
133  // END EXCEPTION CHECKS
134 
135  return rA.trace();
136 }
137 
145 template<typename Derived>
146 typename Derived::Scalar det(const Eigen::MatrixBase<Derived>& A)
147 {
148  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
149 
150  // EXCEPTION CHECKS
151 
152  // check zero-size
154  throw exception::ZeroSize("qpp::det()");
155  // END EXCEPTION CHECKS
156 
157  return rA.determinant();
158 }
159 
169 template<typename Derived>
170 typename Derived::Scalar logdet(const Eigen::MatrixBase<Derived>& A)
171 {
172  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
173 
174  // EXCEPTION CHECKS
175 
176  // check zero-size
178  throw exception::ZeroSize("qpp::logdet()");
179 
180  // check square matrix
182  throw exception::MatrixNotSquare("qpp::logdet()");
183  // END EXCEPTION CHECKS
184 
185  Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
187  lu.matrixLU().template triangularView<Eigen::Upper>();
188  typename Derived::Scalar result = std::log(U(0, 0));
189 
190  for (idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
191  result += std::log(U(i, i));
192 
193  return result;
194 
195 }
196 
204 template<typename Derived>
205 typename Derived::Scalar sum(const Eigen::MatrixBase<Derived>& A)
206 {
207  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
208 
209  // EXCEPTION CHECKS
210 
211  // check zero-size
213  throw exception::ZeroSize("qpp::sum()");
214  // END EXCEPTION CHECKS
215 
216  return rA.sum();
217 }
218 
226 template<typename Derived>
227 typename Derived::Scalar prod(const Eigen::MatrixBase<Derived>& A)
228 {
229  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
230 
231  // EXCEPTION CHECKS
232 
233  // check zero-size
235  throw exception::ZeroSize("qpp::prod()");
236  // END EXCEPTION CHECKS
237 
238  return rA.prod();
239 }
240 
247 template<typename Derived>
248 double norm(const Eigen::MatrixBase<Derived>& A)
249 {
250  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
251 
252  // EXCEPTION CHECKS
253 
254  // check zero-size
256  throw exception::ZeroSize("qpp::norm()");
257  // END EXCEPTION CHECKS
258 
259  // convert matrix to complex then return its norm
260  return (rA.template cast<cplx>()).norm();
261 }
262 
271 template<typename Derived>
272 std::pair<dyn_col_vect < cplx>, cmat>
273 
274 eig(const Eigen::MatrixBase<Derived>& A)
275 {
276  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
277 
278  // EXCEPTION CHECKS
279 
280  // check zero-size
282  throw exception::ZeroSize("qpp::eig()");
283 
284  // check square matrix
286  throw exception::MatrixNotSquare("qpp::eig()");
287  // END EXCEPTION CHECKS
288 
289  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
290 
291  return std::make_pair(es.eigenvalues(), es.eigenvectors());
292 }
293 
301 template<typename Derived>
302 dyn_col_vect <cplx> evals(const Eigen::MatrixBase<Derived>& A)
303 {
304  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
305 
306  // EXCEPTION CHECKS
307 
308  // check zero-size
310  throw exception::ZeroSize("qpp::evals()");
311 
312  // check square matrix
314  throw exception::MatrixNotSquare("qpp::evals()");
315  // END EXCEPTION CHECKS
316 
317  return eig(rA).first;
318 }
319 
327 template<typename Derived>
328 cmat evects(const Eigen::MatrixBase<Derived>& A)
329 {
330  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
331 
332  // EXCEPTION CHECKS
333 
334  // check zero-size
336  throw exception::ZeroSize("qpp::evects()");
337 
338  // check square matrix
340  throw exception::MatrixNotSquare("qpp::evects()");
341  // END EXCEPTION CHECKS
342 
343  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
344 
345  return eig(rA).second;
346 }
347 
356 template<typename Derived>
357 std::pair<dyn_col_vect < double>, cmat>
358 
359 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::ZeroSize("qpp::heig()");
368 
369  // check square matrix
371  throw exception::MatrixNotSquare("qpp::heig()");
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::ZeroSize("qpp::hevals()");
396 
397  // check square matrix
399  throw exception::MatrixNotSquare("qpp::hevals()");
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::ZeroSize("qpp::hevects()");
422 
423  // check square matrix
425  throw exception::MatrixNotSquare("qpp::hevects()");
426  // END EXCEPTION CHECKS
427 
428  return heig(rA).second;
429 }
430 
440 template<typename Derived>
441 std::tuple<cmat, dyn_col_vect < double>, cmat>
442 
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::ZeroSize("qpp::svd()");
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::ZeroSize("qpp::svals()");
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::ZeroSize("qpp::svdU()");
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::ZeroSize("qpp::svdV()");
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::ZeroSize("qpp::funm()");
555 
556  // check square matrix
558  throw exception::MatrixNotSquare("qpp::funm()");
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::ZeroSize("qpp::sqrtm()");
588 
589  // check square matrix
591  throw exception::MatrixNotSquare("qpp::sqrtm()");
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::ZeroSize("qpp::absm()");
613 
614  // check square matrix
616  throw exception::MatrixNotSquare("qpp::absm()");
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::ZeroSize("qpp::expm()");
638 
639  // check square matrix
641  throw exception::MatrixNotSquare("qpp::expm()");
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::ZeroSize("qpp::logm()");
663 
664  // check square matrix
666  throw exception::MatrixNotSquare("qpp::logm()");
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::ZeroSize("qpp::sinm()");
688 
689  // check square matrix
691  throw exception::MatrixNotSquare("qpp::sinm()");
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::ZeroSize("qpp::cosm()");
713 
714  // check square matrix
716  throw exception::MatrixNotSquare("qpp::cosm()");
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::ZeroSize("qpp::spectralpowm()");
743 
744  // check square matrix
746  throw exception::MatrixNotSquare("qpp::spectralpowm()");
747  // END EXCEPTION CHECKS
748 
749  // Define A^0 = Id, for z IDENTICALLY zero
750  if (real(z) == 0 && imag(z) == 0)
751  return cmat::Identity(rA.rows(), rA.rows());
752 
753  Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
754  cmat evects = es.eigenvectors();
755  cmat evals = es.eigenvalues();
756  for (idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
757  evals(i) = std::pow(evals(i), z);
758 
759  cmat evalsdiag = evals.asDiagonal();
760 
761  return evects * evalsdiag * evects.inverse();
762 }
763 
776 template<typename Derived>
777 dyn_mat<typename Derived::Scalar> powm(const Eigen::MatrixBase<Derived>& A,
778  idx n)
779 {
780  // EXCEPTION CHECKS
781 
782  // check zero-size
784  throw exception::ZeroSize("qpp::powm()");
785 
786  // check square matrix
788  throw exception::MatrixNotSquare("qpp::powm()");
789  // END EXCEPTION CHECKS
790 
791  // if n = 1, return the matrix unchanged
792  if (n == 1)
793  return A;
794 
797  A.rows(), A.rows());
798 
799  // if n = 0, return the identity (as just prepared in result)
800  if (n == 0)
801  return result;
802 
803  dyn_mat<typename Derived::Scalar> cA = A.derived(); // copy
804 
805  // fast matrix power
806  for (; n > 0; n /= 2)
807  {
808  if (n % 2)
809  result = (result * cA).eval();
810  cA = (cA * cA).eval();
811  }
812 
813  return result;
814 }
815 
824 template<typename Derived>
825 double schatten(const Eigen::MatrixBase<Derived>& A, double p)
826 {
827  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
828 
829  // EXCEPTION CHECKS
830 
831  // check zero-size
833  throw exception::ZeroSize("qpp::schatten()");
834  if (p < 1)
835  throw exception::OutOfRange("qpp::schatten()");
836  // END EXCEPTION CHECKS
837 
838  if (p == infty) // infinity norm (largest singular value)
839  return svals(rA)(0);
840 
841  const dyn_col_vect<double> sv = svals(rA);
842  double result = 0;
843  for (idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
844  result += std::pow(sv[i], p);
845 
846  return std::pow(result, 1. / p);
847 }
848 
849 // other functions
850 
859 template<typename OutputScalar, typename Derived>
860 dyn_mat <OutputScalar> cwise(const Eigen::MatrixBase<Derived>& A,
861  OutputScalar (* f)(
862  const typename Derived::Scalar&))
863 {
864  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
865 
866  // EXCEPTION CHECKS
867 
868  // check zero-size
870  throw exception::ZeroSize("qpp::cwise()");
871  // END EXCEPTION CHECKS
872 
873  dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
874 
875 #ifdef WITH_OPENMP_
876 #pragma omp parallel for collapse(2)
877 #endif // WITH_OPENMP_
878  // column major order for speed
879  for (idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
880  for (idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
881  result(i, j) = (*f)(rA(i, j));
882 
883  return result;
884 }
885 
886 // Kronecker product of multiple matrices, preserve return type
887 // variadic template
898 template<typename T>
900 {
901  return head;
902 }
903 
914 template<typename T, typename ... Args>
915 dyn_mat<typename T::Scalar> kron(const T& head, const Args& ... tail)
916 {
917  return internal::kron2(head, kron(tail...));
918 }
919 
929 template<typename Derived>
930 dyn_mat<typename Derived::Scalar> kron(const std::vector<Derived>& As)
931 {
932  // EXCEPTION CHECKS
933 
934  if (As.size() == 0)
935  throw exception::ZeroSize("qpp::kron()");
936 
937  for (auto&& it : As)
939  throw exception::ZeroSize("qpp::kron()");
940  // END EXCEPTION CHECKS
941 
942  dyn_mat<typename Derived::Scalar> result = As[0].derived();
943  for (idx i = 1; i < As.size(); ++i)
944  {
945  result = kron(result, As[i]);
946  }
947 
948  return result;
949 }
950 
951 // Kronecker product of a list of matrices, preserve return type
952 // deduce the template parameters from initializer_list
963 template<typename Derived>
965  const std::initializer_list<Derived>& As)
966 {
967  return kron(std::vector<Derived>(As));
968 }
969 
979 template<typename Derived>
980 dyn_mat<typename Derived::Scalar> kronpow(const Eigen::MatrixBase<Derived>& A,
981  idx n)
982 {
983  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
984 
985  // EXCEPTION CHECKS
986 
987  // check zero-size
989  throw exception::ZeroSize("qpp::kronpow()");
990 
991  // check out of range
992  if (n == 0)
993  throw exception::OutOfRange("qpp::kronpow()");
994  // END EXCEPTION CHECKS
995 
996  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
997 
998  return kron(As);
999 }
1000 
1001 // Direct sum of multiple matrices, preserve return type
1002 // variadic template
1013 template<typename T>
1015 {
1016  return head;
1017 }
1018 
1029 template<typename T, typename ... Args>
1030 dyn_mat<typename T::Scalar> dirsum(const T& head, const Args& ... tail)
1031 {
1032  return internal::dirsum2(head, dirsum(tail...));
1033 }
1034 
1044 template<typename Derived>
1045 dyn_mat<typename Derived::Scalar> dirsum(const std::vector<Derived>& As)
1046 {
1047  // EXCEPTION CHECKS
1048 
1049  if (As.size() == 0)
1050  throw exception::ZeroSize("qpp::dirsum()");
1051 
1052  for (auto&& it : As)
1054  throw exception::ZeroSize("qpp::dirsum()");
1055  // END EXCEPTION CHECKS
1056 
1057  idx total_rows = 0, total_cols = 0;
1058  for (idx i = 0; i < As.size(); ++i)
1059  {
1060  total_rows += static_cast<idx>(As[i].rows());
1061  total_cols += static_cast<idx>(As[i].cols());
1062  }
1064  dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1065 
1066  idx cur_row = 0, cur_col = 0;
1067  for (idx i = 0; i < As.size(); ++i)
1068  {
1069  result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1070  cur_row += static_cast<idx>(As[i].rows());
1071  cur_col += static_cast<idx>(As[i].cols());
1072  }
1073 
1074  return result;
1075 }
1076 
1077 // Direct sum of a list of matrices, preserve return type
1078 // deduce the template parameters from initializer_list
1089 template<typename Derived>
1091  const std::initializer_list<Derived>& As)
1092 {
1093  return dirsum(std::vector<Derived>(As));
1094 }
1095 
1105 template<typename Derived>
1107  const Eigen::MatrixBase<Derived>& A, idx n)
1108 {
1109  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1110 
1111  // EXCEPTION CHECKS
1112 
1113  // check zero-size
1115  throw exception::ZeroSize("qpp::dirsumpow()");
1116 
1117  // check out of range
1118  if (n == 0)
1119  throw exception::OutOfRange("qpp::dirsumpow()");
1120  // END EXCEPTION CHECKS
1121 
1122  std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1123 
1124  return dirsum(As);
1125 }
1126 
1138 template<typename Derived>
1139 dyn_mat<typename Derived::Scalar> reshape(const Eigen::MatrixBase<Derived>& A,
1140  idx rows, idx cols)
1141 {
1142  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1143 
1144  idx Arows = static_cast<idx>(rA.rows());
1145  idx Acols = static_cast<idx>(rA.cols());
1146 
1147  // EXCEPTION CHECKS
1148 
1149  // check zero-size
1151  throw exception::ZeroSize("qpp::reshape()");
1152 
1153  if (Arows * Acols != rows * cols)
1154  throw exception::DimsMismatchMatrix("qpp::reshape()");
1155  // END EXCEPTION CHECKS
1156 
1157  return Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1158  const_cast<typename Derived::Scalar*>(rA.data()), rows, cols);
1159 }
1160 
1173 template<typename Derived1, typename Derived2>
1174 dyn_mat<typename Derived1::Scalar> comm(const Eigen::MatrixBase<Derived1>& A,
1175  const Eigen::MatrixBase<Derived2>& B)
1176 {
1177  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1178  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1179 
1180  // EXCEPTION CHECKS
1181 
1182  // check types
1183  if (!std::is_same<typename Derived1::Scalar,
1184  typename Derived2::Scalar>::value)
1185  throw exception::TypeMismatch("qpp::comm()");
1186 
1187  // check zero-size
1190  throw exception::ZeroSize("qpp::comm()");
1191 
1192  // check square matrices
1194  throw exception::MatrixNotSquare("qpp::comm()");
1195 
1196  // check equal dimensions
1197  if (rA.rows() != rB.rows())
1198  throw exception::DimsNotEqual("qpp::comm()");
1199  // END EXCEPTION CHECKS
1200 
1201  return rA * rB - rB * rA;
1202 }
1203 
1216 template<typename Derived1, typename Derived2>
1218  const Eigen::MatrixBase<Derived1>& A,
1219  const Eigen::MatrixBase<Derived2>& B)
1220 {
1221  const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1222  const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1223 
1224  // EXCEPTION CHECKS
1225 
1226  // check types
1227  if (!std::is_same<typename Derived1::Scalar,
1228  typename Derived2::Scalar>::value)
1229  throw exception::TypeMismatch("qpp::anticomm()");
1230 
1231  // check zero-size
1234  throw exception::ZeroSize("qpp::anticomm()");
1235 
1236  // check square matrices
1238  throw exception::MatrixNotSquare("qpp::anticomm()");
1239 
1240  // check equal dimensions
1241  if (rA.rows() != rB.rows())
1242  throw exception::DimsNotEqual("qpp::anticomm()");
1243  // END EXCEPTION CHECKS
1244 
1245  return rA * rB + rB * rA;
1246 }
1247 
1258 template<typename Derived>
1259 dyn_mat<typename Derived::Scalar> prj(const Eigen::MatrixBase<Derived>& A)
1260 {
1261  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1262 
1263  // EXCEPTION CHECKS
1264 
1265  // check zero-size
1267  throw exception::ZeroSize("qpp::prj()");
1268 
1269  // check column vector
1270  if (!internal::check_cvector(rA))
1271  throw exception::MatrixNotCvector("qpp::prj()");
1272  // END EXCEPTION CHECKS
1273 
1274  double normA = norm(rA);
1275  if (normA > eps)
1276  return rA * adjoint(rA) / (normA * normA);
1277  else
1279  ::Zero(rA.rows(), rA.rows());
1280 
1281 }
1282 
1290 template<typename Derived>
1291 dyn_mat<typename Derived::Scalar> grams(const std::vector<Derived>& As)
1292 {
1293  // EXCEPTION CHECKS
1294 
1295  // check empty list
1297  throw exception::ZeroSize("qpp::grams()");
1298 
1299  for (auto&& it : As)
1301  throw exception::ZeroSize("qpp::grams()");
1302 
1303  // check that As[0] is a column vector
1304  if (!internal::check_cvector(As[0]))
1305  throw exception::MatrixNotCvector("qpp::grams()");
1306 
1307  // now check that all the rest match the size of the first vector
1308  for (auto&& it : As)
1309  if (it.rows() != As[0].rows() || it.cols() != 1)
1310  throw exception::DimsNotEqual("qpp::grams()");
1311  // END EXCEPTION CHECKS
1312 
1315  As[0].rows());
1316 
1318  dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1319 
1320  std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1321  // find the first non-zero vector in the list
1322  idx pos = 0;
1323  for (pos = 0; pos < As.size(); ++pos)
1324  {
1325  if (norm(As[pos]) > eps) // add it as the first element
1326  {
1327  outvecs.push_back(As[pos]);
1328  break;
1329  }
1330  }
1331 
1332  // start the process
1333  for (idx i = pos + 1; i < As.size(); ++i)
1334  {
1335  cut -= prj(outvecs[i - 1 - pos]);
1336  vi = cut * As[i];
1337  outvecs.push_back(vi);
1338  }
1339 
1340  dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1341 
1342  idx cnt = 0;
1343  for (auto&& it : outvecs)
1344  {
1345  double normA = norm(it);
1346  if (normA > eps) // we add only the non-zero vectors
1347  {
1348  result.col(cnt) = it / normA;
1349  cnt++;
1350  }
1351  }
1352 
1353  return result.block(0, 0, As[0].rows(), cnt);
1354 }
1355 
1356 // deduce the template parameters from initializer_list
1364 template<typename Derived>
1366  const std::initializer_list<Derived>& As)
1367 {
1368  return grams(std::vector<Derived>(As));
1369 }
1370 
1378 template<typename Derived>
1379 dyn_mat<typename Derived::Scalar> grams(const Eigen::MatrixBase<Derived>& A)
1380 {
1381  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1382 
1383  // EXCEPTION CHECKS
1384 
1386  throw exception::ZeroSize("qpp::grams()");
1387  // END EXCEPTION CHECKS
1388 
1389  std::vector<dyn_mat<typename Derived::Scalar>> input;
1390 
1391  for (idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1392  input.push_back(rA.col(i));
1393 
1394  return grams<dyn_mat<typename Derived::Scalar>>
1395  (input);
1396 }
1397 
1408 inline std::vector<idx> n2multiidx(idx n, const std::vector<idx>& dims)
1409 {
1410  // EXCEPTION CHECKS
1411 
1412  if (!internal::check_dims(dims))
1413  throw exception::DimsInvalid("qpp::n2multiidx()");
1414 
1415  if (n >= std::accumulate(std::begin(dims), std::end(dims),
1416  static_cast<idx>(1), std::multiplies<idx>()))
1417  throw exception::OutOfRange("qpp::n2multiidx()");
1418  // END EXCEPTION CHECKS
1419 
1420  // double the size for matrices reshaped as vectors
1421  idx result[2 * maxn];
1422  internal::n2multiidx(n, dims.size(), dims.data(), result);
1423 
1424  return std::vector<idx>(result, result + dims.size());
1425 }
1426 
1437 inline idx multiidx2n(const std::vector<idx>& midx,
1438  const std::vector<idx>& dims)
1439 {
1440  // EXCEPTION CHECKS
1441 
1442  if (!internal::check_dims(dims))
1443  throw exception::DimsInvalid("qpp::multiidx2n()");
1444 
1445  for (idx i = 0; i < dims.size(); ++i)
1446  if (midx[i] >= dims[i])
1447  throw exception::OutOfRange("qpp::multiidx2n()");
1448  // END EXCEPTION CHECKS
1449 
1450  return internal::multiidx2n(midx.data(), dims.size(), dims.data());
1451 }
1452 
1465 inline ket mket(const std::vector<idx>& mask,
1466  const std::vector<idx>& dims)
1467 {
1468  idx N = mask.size();
1469 
1470  idx D = std::accumulate(std::begin(dims), std::end(dims),
1471  static_cast<idx>(1), std::multiplies<idx>());
1472 
1473  // EXCEPTION CHECKS
1474 
1475  // check zero size
1476  if (N == 0)
1477  throw exception::ZeroSize("qpp::mket()");
1478  // check valid dims
1479  if (!internal::check_dims(dims))
1480  throw exception::DimsInvalid("qpp::mket()");
1481  // check mask and dims have the same size
1482  if (mask.size() != dims.size())
1483  throw exception::SubsysMismatchDims("qpp::mket()");
1484  // check mask is a valid vector
1485  for (idx i = 0; i < N; ++i)
1486  if (mask[i] >= dims[i])
1487  throw exception::SubsysMismatchDims("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::ZeroSize("qpp::mket()");
1519 
1520  // check valid dims
1521  if (d == 0)
1522  throw exception::DimsInvalid("qpp::mket()");
1523 
1524  // check mask is a valid vector
1525  for (idx i = 0; i < N; ++i)
1526  if (mask[i] >= d)
1527  throw exception::SubsysMismatchDims("qpp::mket()");
1528  // END EXCEPTION CHECKS
1529 
1530  ket result = ket::Zero(D);
1531  std::vector<idx> dims(N, d);
1532  idx pos = multiidx2n(mask, dims);
1533  result(pos) = 1;
1534 
1535  return result;
1536 }
1537 
1552 inline cmat mprj(const std::vector<idx>& mask,
1553  const std::vector<idx>& dims)
1554 {
1555  idx N = mask.size();
1556 
1557  idx D = std::accumulate(std::begin(dims), std::end(dims),
1558  static_cast<idx>(1), std::multiplies<idx>());
1559 
1560  // EXCEPTION CHECKS
1561 
1562  // check zero size
1563  if (N == 0)
1564  throw exception::ZeroSize("qpp::mprj()");
1565  // check valid dims
1566  if (!internal::check_dims(dims))
1567  throw exception::DimsInvalid("qpp::mprj()");
1568  // check mask and dims have the same size
1569  if (mask.size() != dims.size())
1570  throw exception::SubsysMismatchDims("qpp::mprj()");
1571  // check mask is a valid vector
1572  for (idx i = 0; i < N; ++i)
1573  if (mask[i] >= dims[i])
1574  throw exception::SubsysMismatchDims("qpp::mprj()");
1575  // END EXCEPTION CHECKS
1576 
1577  cmat result = cmat::Zero(D, D);
1578  idx pos = multiidx2n(mask, dims);
1579  result(pos, pos) = 1;
1580 
1581  return result;
1582 }
1583 
1598 inline cmat mprj(const std::vector<idx>& mask, idx d = 2)
1599 {
1600  idx N = mask.size();
1601  idx D = static_cast<idx>(std::llround(std::pow(d, N)));
1602 
1603  // EXCEPTION CHECKS
1604 
1605  // check zero size
1606  if (N == 0)
1607  throw exception::ZeroSize("qpp::mprj()");
1608 
1609  // check valid dims
1610  if (d == 0)
1611  throw exception::DimsInvalid("qpp::mprj()");
1612 
1613  // check mask is a valid vector
1614  for (idx i = 0; i < N; ++i)
1615  if (mask[i] >= d)
1616  throw exception::SubsysMismatchDims("qpp::mprj()");
1617  // END EXCEPTION CHECKS
1618 
1619  cmat result = cmat::Zero(D, D);
1620  std::vector<idx> dims(N, d);
1621  idx pos = multiidx2n(mask, dims);
1622  result(pos, pos) = 1;
1623 
1624  return result;
1625 }
1626 
1635 template<typename InputIterator>
1636 std::vector<double> abssq(InputIterator first, InputIterator last)
1637 {
1638  std::vector<double> weights(std::distance(first, last));
1639  std::transform(first, last, std::begin(weights),
1640  [](cplx z) -> double
1641  {
1642  return std::norm(z);
1643  });
1644 
1645  return weights;
1646 }
1647 
1654 template<typename Container>
1655 std::vector<double> abssq(const Container& c,
1656  typename std::enable_if<
1658  >::type* = nullptr)
1659 // we need the std::enable_if to SFINAE out Eigen expressions
1660 // that will otherwise match, instead of matching
1661 // the overload below:
1662 
1663 // template<typename Derived>
1664 // abssq(const Eigen::MatrixBase<Derived>& A)
1665 {
1666  return abssq(std::begin(c), std::end(c));
1667 }
1668 
1675 template<typename Derived>
1676 std::vector<double> abssq(const Eigen::MatrixBase<Derived>& A)
1677 {
1678  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1679 
1680  // EXCEPTION CHECKS
1681 
1682  // check zero-size
1684  throw exception::ZeroSize("qpp::abssq()");
1685  // END EXCEPTION CHECKS
1686 
1687  return abssq(rA.data(), rA.data() + rA.size());
1688 }
1689 
1698 template<typename InputIterator>
1699 typename std::iterator_traits<InputIterator>::value_type
1700 sum(InputIterator first, InputIterator last)
1701 {
1702  using value_type =
1703  typename std::iterator_traits<InputIterator>::value_type;
1704 
1705  return std::accumulate(first, last, static_cast<value_type>(0));
1706 }
1707 
1715 template<typename Container>
1716 typename Container::value_type
1717 sum(const Container& c,
1718  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1719 {
1720  return sum(std::begin(c), std::end(c));
1721 }
1722 
1731 template<typename InputIterator>
1732 typename std::iterator_traits<InputIterator>::value_type
1733 prod(InputIterator first, InputIterator last)
1734 {
1735  using value_type =
1736  typename std::iterator_traits<InputIterator>::value_type;
1737 
1738  return std::accumulate(first, last, static_cast<value_type>(1),
1739  std::multiplies<value_type>());
1740 }
1741 
1749 template<typename Container>
1750 typename Container::value_type
1751 prod(const Container& c,
1752  typename std::enable_if<is_iterable<Container>::value>::type* = nullptr)
1753 {
1754  return prod(std::begin(c), std::end(c));
1755 }
1756 
1769 template<typename Derived>
1771  const Eigen::MatrixBase<Derived>& A)
1772 {
1773  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1774 
1775  // EXCEPTION CHECKS
1776  // check zero-size
1778  throw exception::ZeroSize("qpp::rho2pure()");
1779  // check square matrix
1780  if (!internal::check_square_mat(rA))
1781  throw exception::MatrixNotSquare("qpp::rho2pure()");
1782  // END EXCEPTION CHECKS
1783 
1784  dyn_col_vect<double> tmp_evals = hevals(rA);
1785  cmat tmp_evects = hevects(rA);
1788  // find the non-zero eigenvector
1789  // there is only one, assuming the state is pure
1790  for (idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1791  {
1792  if (std::abs(tmp_evals(k)) > eps)
1793  {
1794  result = tmp_evects.col(k);
1795  break;
1796  }
1797  }
1798 
1799  return result;
1800 }
1801 
1810 template<typename T>
1811 std::vector<T> complement(std::vector<T> subsys, idx N)
1812 {
1813  // EXCEPTION CHECKS
1814 
1815  if (N < subsys.size())
1816  throw exception::OutOfRange("qpp::complement()");
1817  // END EXCEPTION CHECKS
1818 
1819  std::vector<T> all(N);
1820  std::vector<T> subsys_bar(N - subsys.size());
1821 
1822  std::iota(std::begin(all), std::end(all), 0);
1823  std::sort(std::begin(subsys), std::end(subsys));
1824  std::set_difference(std::begin(all), std::end(all),
1825  std::begin(subsys), std::end(subsys),
1826  std::begin(subsys_bar));
1827 
1828  return subsys_bar;
1829 }
1830 
1841 template<typename Derived>
1842 std::vector<double> rho2bloch(const Eigen::MatrixBase<Derived>& A)
1843 {
1844  const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1845 
1846  // EXCEPTION CHECKS
1847 
1848  // check qubit matrix
1850  throw exception::NotQubitMatrix("qpp::rho2bloch()");
1851  // END EXCEPTION CHECKS
1852 
1853  std::vector<double> result(3);
1854  cmat X(2, 2), Y(2, 2), Z(2, 2);
1855  X << 0, 1, 1, 0;
1856  Y << 0, -1_i, 1_i, 0;
1857  Z << 1, 0, 0, -1;
1858  result[0] = std::real(trace(rA * X));
1859  result[1] = std::real(trace(rA * Y));
1860  result[2] = std::real(trace(rA * Z));
1861 
1862  return result;
1863 }
1864 
1873 inline cmat bloch2rho(const std::vector<double>& r)
1874 {
1875  // EXCEPTION CHECKS
1876 
1877  // check 3-dimensional vector
1878  if (r.size() != 3)
1879  throw exception::CustomException("qpp::bloch2rho",
1880  "r is not a 3-dimensional vector!");
1881  // END EXCEPTION CHECKS
1882 
1883  cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1884  X << 0, 1, 1, 0;
1885  Y << 0, -1_i, 1_i, 0;
1886  Z << 1, 0, 0, -1;
1887  Id2 << 1, 0, 0, 1;
1888 
1889  return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
1890 }
1891 
1892 } /* namespace qpp */
1893 
1894 #endif /* FUNCTIONS_H_ */
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:129
Dimensions not equal exception.
Definition: exception.h:304
Dimension(s) mismatch matrix size exception.
Definition: exception.h:322
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:170
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:71
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:293
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:122
Custom exception.
Definition: exception.h:636
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:227
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:64
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:1014
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:77
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:70
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:40
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:1842
Subsystems mismatch dimensions exception.
Definition: exception.h:395
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:50
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1217
Matrix is not a column vector exception.
Definition: exception.h:168
Quantum++ main namespace.
Definition: codes.h:30
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:248
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
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:860
Invalid dimension(s) exception.
Definition: exception.h:287
Checks whether T is compatible with an STL-like iterable container.
Definition: traits.h:63
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:1770
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:1106
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1811
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:328
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:825
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:1873
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:274
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1408
dyn_mat< typename Derived::Scalar > adjoint(const Eigen::MatrixBase< Derived > &A)
Adjoint.
Definition: functions.h:81
idx multiidx2n(const std::vector< idx > &midx, const std::vector< idx > &dims)
Multi-index to non-negative integer index.
Definition: functions.h:1437
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:470
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:443
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:899
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:45
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1552
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:146
dyn_mat< typename Derived::Scalar > transpose(const Eigen::MatrixBase< Derived > &A)
Transpose.
Definition: functions.h:41
dyn_mat< typename Derived::Scalar > kronpow(const Eigen::MatrixBase< Derived > &A, idx n)
Kronecker power.
Definition: functions.h:980
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &A)
Projector.
Definition: functions.h:1259
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:777
Type mismatch exception.
Definition: exception.h:585
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:302
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:61
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:85
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:101
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:124
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:519
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:336
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:704
Parameter out of range exception.
Definition: exception.h:567
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:89
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:1636
Matrix is not 2 x 2 exception.
Definition: exception.h:447
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &As)
Gram-Schmidt orthogonalization.
Definition: functions.h:1291
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1174
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:60
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1465
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:103
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:654
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:142
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:35
Matrix is not square exception.
Definition: exception.h:151
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:205
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1139
Object has zero size exception.
Definition: exception.h:134
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:359
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:249