44 template <
typename Derived>
45 dyn_mat<typename Derived::Scalar>
46 transpose(
const Eigen::MatrixBase<Derived>& A) {
47 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
51 throw exception::ZeroSize(
"qpp::transpose()");
53 return rA.transpose();
63 template <
typename Derived>
64 dyn_mat<typename Derived::Scalar>
65 conjugate(
const Eigen::MatrixBase<Derived>& A) {
66 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
70 throw exception::ZeroSize(
"qpp::conjugate()");
72 return rA.conjugate();
82 template <
typename Derived>
83 dyn_mat<typename Derived::Scalar> adjoint(
const Eigen::MatrixBase<Derived>& A) {
84 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
90 throw exception::ZeroSize(
"qpp::adjoint()");
103 template <
typename Derived>
104 dyn_mat<typename Derived::Scalar> inverse(
const Eigen::MatrixBase<Derived>& A) {
105 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
111 throw exception::ZeroSize(
"qpp::inverse()");
123 template <
typename Derived>
124 typename Derived::Scalar trace(
const Eigen::MatrixBase<Derived>& A) {
125 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
131 throw exception::ZeroSize(
"qpp::trace()");
144 template <
typename Derived>
145 typename Derived::Scalar det(
const Eigen::MatrixBase<Derived>& A) {
146 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
152 throw exception::ZeroSize(
"qpp::det()");
155 return rA.determinant();
167 template <
typename Derived>
168 typename Derived::Scalar logdet(
const Eigen::MatrixBase<Derived>& A) {
169 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
175 throw exception::ZeroSize(
"qpp::logdet()");
179 throw exception::MatrixNotSquare(
"qpp::logdet()");
182 Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
183 dyn_mat<typename Derived::Scalar> U =
184 lu.matrixLU().template triangularView<Eigen::Upper>();
185 typename Derived::Scalar result = std::log(U(0, 0));
187 for (
idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
188 result += std::log(U(i, i));
200 template <
typename Derived>
201 typename Derived::Scalar sum(
const Eigen::MatrixBase<Derived>& A) {
202 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
208 throw exception::ZeroSize(
"qpp::sum()");
221 template <
typename Derived>
222 typename Derived::Scalar prod(
const Eigen::MatrixBase<Derived>& A) {
223 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
229 throw exception::ZeroSize(
"qpp::prod()");
241 template <
typename Derived>
242 double norm(
const Eigen::MatrixBase<Derived>& A) {
243 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
249 throw exception::ZeroSize(
"qpp::norm()");
253 return (rA.template cast<cplx>()).norm();
264 template <
typename Derived>
265 std::pair<dyn_col_vect<cplx>,
cmat>
267 eig(
const Eigen::MatrixBase<Derived>& A) {
268 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
274 throw exception::ZeroSize(
"qpp::eig()");
278 throw exception::MatrixNotSquare(
"qpp::eig()");
281 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
283 return std::make_pair(es.eigenvalues(), es.eigenvectors());
293 template <
typename Derived>
294 dyn_col_vect<cplx> evals(
const Eigen::MatrixBase<Derived>& A) {
295 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
301 throw exception::ZeroSize(
"qpp::evals()");
305 throw exception::MatrixNotSquare(
"qpp::evals()");
308 return eig(rA).first;
318 template <
typename Derived>
319 cmat evects(
const Eigen::MatrixBase<Derived>& A) {
320 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
326 throw exception::ZeroSize(
"qpp::evects()");
330 throw exception::MatrixNotSquare(
"qpp::evects()");
333 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
335 return eig(rA).second;
346 template <
typename Derived>
347 std::pair<dyn_col_vect<double>,
cmat>
349 heig(
const Eigen::MatrixBase<Derived>& A) {
350 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
356 throw exception::ZeroSize(
"qpp::heig()");
360 throw exception::MatrixNotSquare(
"qpp::heig()");
363 Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
365 return std::make_pair(es.eigenvalues(), es.eigenvectors());
375 template <
typename Derived>
376 dyn_col_vect<double> hevals(
const Eigen::MatrixBase<Derived>& A) {
377 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
383 throw exception::ZeroSize(
"qpp::hevals()");
387 throw exception::MatrixNotSquare(
"qpp::hevals()");
390 return heig(rA).first;
400 template <
typename Derived>
401 cmat hevects(
const Eigen::MatrixBase<Derived>& A) {
402 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
408 throw exception::ZeroSize(
"qpp::hevects()");
412 throw exception::MatrixNotSquare(
"qpp::hevects()");
415 return heig(rA).second;
427 template <
typename Derived>
428 std::tuple<cmat, dyn_col_vect<double>,
cmat>
430 svd(
const Eigen::MatrixBase<Derived>& A) {
431 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
437 throw exception::ZeroSize(
"qpp::svd()");
440 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
441 rA, Eigen::DecompositionOptions::ComputeFullU |
442 Eigen::DecompositionOptions::ComputeFullV);
444 return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
454 template <
typename Derived>
455 dyn_col_vect<double> svals(
const Eigen::MatrixBase<Derived>& A) {
456 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
462 throw exception::ZeroSize(
"qpp::svals()");
465 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
467 return sv.singularValues();
477 template <
typename Derived>
478 cmat svdU(
const Eigen::MatrixBase<Derived>& A) {
479 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
485 throw exception::ZeroSize(
"qpp::svdU()");
488 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
489 rA, Eigen::DecompositionOptions::ComputeFullU);
501 template <
typename Derived>
502 cmat svdV(
const Eigen::MatrixBase<Derived>& A) {
503 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
509 throw exception::ZeroSize(
"qpp::svdV()");
512 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(
513 rA, Eigen::DecompositionOptions::ComputeFullV);
527 template <
typename Derived>
528 cmat funm(
const Eigen::MatrixBase<Derived>& A,
cplx (*f)(
const cplx&)) {
529 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
535 throw exception::ZeroSize(
"qpp::funm()");
539 throw exception::MatrixNotSquare(
"qpp::funm()");
542 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
543 cmat evects = es.eigenvectors();
544 cmat evals = es.eigenvalues();
545 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
546 evals(i) = (*f)(evals(i));
548 cmat evalsdiag = evals.asDiagonal();
550 return evects * evalsdiag * evects.inverse();
559 template <
typename Derived>
560 cmat sqrtm(
const Eigen::MatrixBase<Derived>& A) {
561 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
567 throw exception::ZeroSize(
"qpp::sqrtm()");
571 throw exception::MatrixNotSquare(
"qpp::sqrtm()");
574 return funm(rA, &std::sqrt);
583 template <
typename Derived>
584 cmat absm(
const Eigen::MatrixBase<Derived>& A) {
585 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
591 throw exception::ZeroSize(
"qpp::absm()");
595 throw exception::MatrixNotSquare(
"qpp::absm()");
598 return sqrtm(adjoint(rA) * rA);
607 template <
typename Derived>
608 cmat expm(
const Eigen::MatrixBase<Derived>& A) {
609 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
615 throw exception::ZeroSize(
"qpp::expm()");
619 throw exception::MatrixNotSquare(
"qpp::expm()");
622 return funm(rA, &std::exp);
631 template <
typename Derived>
632 cmat logm(
const Eigen::MatrixBase<Derived>& A) {
633 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
639 throw exception::ZeroSize(
"qpp::logm()");
643 throw exception::MatrixNotSquare(
"qpp::logm()");
646 return funm(rA, &std::log);
655 template <
typename Derived>
656 cmat sinm(
const Eigen::MatrixBase<Derived>& A) {
657 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
663 throw exception::ZeroSize(
"qpp::sinm()");
667 throw exception::MatrixNotSquare(
"qpp::sinm()");
670 return funm(rA, &std::sin);
679 template <
typename Derived>
680 cmat cosm(
const Eigen::MatrixBase<Derived>& A) {
681 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
687 throw exception::ZeroSize(
"qpp::cosm()");
691 throw exception::MatrixNotSquare(
"qpp::cosm()");
694 return funm(rA, &std::cos);
708 template <
typename Derived>
709 cmat spectralpowm(
const Eigen::MatrixBase<Derived>& A,
const cplx z) {
710 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
716 throw exception::ZeroSize(
"qpp::spectralpowm()");
720 throw exception::MatrixNotSquare(
"qpp::spectralpowm()");
724 if (real(z) == 0 && imag(z) == 0)
725 return cmat::Identity(rA.rows(), rA.rows());
727 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
728 cmat evects = es.eigenvectors();
729 cmat evals = es.eigenvalues();
730 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
731 evals(i) = std::pow(evals(i), z);
733 cmat evalsdiag = evals.asDiagonal();
735 return evects * evalsdiag * evects.inverse();
750 template <
typename Derived>
751 dyn_mat<typename Derived::Scalar> powm(
const Eigen::MatrixBase<Derived>& A,
757 throw exception::ZeroSize(
"qpp::powm()");
761 throw exception::MatrixNotSquare(
"qpp::powm()");
768 dyn_mat<typename Derived::Scalar> result =
769 dyn_mat<typename Derived::Scalar>::Identity(A.rows(), A.rows());
775 dyn_mat<typename Derived::Scalar> cA = A.derived();
778 for (; n > 0; n /= 2) {
780 result = (result * cA).eval();
781 cA = (cA * cA).eval();
795 template <
typename Derived>
796 double schatten(
const Eigen::MatrixBase<Derived>& A,
double p) {
797 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
803 throw exception::ZeroSize(
"qpp::schatten()");
805 throw exception::OutOfRange(
"qpp::schatten()");
811 const dyn_col_vect<double> sv = svals(rA);
813 for (
idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
814 result += std::pow(sv[i], p);
816 return std::pow(result, 1. / p);
829 template <
typename OutputScalar,
typename Derived>
830 dyn_mat<OutputScalar>
831 cwise(
const Eigen::MatrixBase<Derived>& A,
832 OutputScalar (*f)(
const typename Derived::Scalar&)) {
833 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
839 throw exception::ZeroSize(
"qpp::cwise()");
842 dyn_mat<OutputScalar> result(rA.rows(), rA.cols());
845 #pragma omp parallel for collapse(2) 846 #endif // WITH_OPENMP_ 848 for (
idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
849 for (
idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
850 result(i, j) = (*f)(rA(i, j));
867 template <
typename T>
868 dyn_mat<typename T::Scalar> kron(
const T& head) {
882 template <
typename T,
typename... Args>
883 dyn_mat<typename T::Scalar> kron(
const T& head,
const Args&... tail) {
896 template <
typename Derived>
897 dyn_mat<typename Derived::Scalar> kron(
const std::vector<Derived>& As) {
901 throw exception::ZeroSize(
"qpp::kron()");
905 throw exception::ZeroSize(
"qpp::kron()");
908 dyn_mat<typename Derived::Scalar> result = As[0].derived();
909 for (
idx i = 1; i < As.size(); ++i) {
910 result = kron(result, As[i]);
928 template <
typename Derived>
929 dyn_mat<typename Derived::Scalar>
930 kron(
const std::initializer_list<Derived>& As) {
931 return kron(std::vector<Derived>(As));
943 template <
typename Derived>
944 dyn_mat<typename Derived::Scalar> kronpow(
const Eigen::MatrixBase<Derived>& A,
946 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
952 throw exception::ZeroSize(
"qpp::kronpow()");
956 throw exception::OutOfRange(
"qpp::kronpow()");
959 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
976 template <
typename T>
977 dyn_mat<typename T::Scalar> dirsum(
const T& head) {
991 template <
typename T,
typename... Args>
992 dyn_mat<typename T::Scalar> dirsum(
const T& head,
const Args&... tail) {
1005 template <
typename Derived>
1006 dyn_mat<typename Derived::Scalar> dirsum(
const std::vector<Derived>& As) {
1010 throw exception::ZeroSize(
"qpp::dirsum()");
1012 for (
auto&& it : As)
1014 throw exception::ZeroSize(
"qpp::dirsum()");
1017 idx total_rows = 0, total_cols = 0;
1018 for (
idx i = 0; i < As.size(); ++i) {
1019 total_rows +=
static_cast<idx>(As[i].rows());
1020 total_cols +=
static_cast<idx>(As[i].cols());
1022 dyn_mat<typename Derived::Scalar> result =
1023 dyn_mat<typename Derived::Scalar>::Zero(total_rows, total_cols);
1025 idx cur_row = 0, cur_col = 0;
1026 for (
idx i = 0; i < As.size(); ++i) {
1027 result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
1028 cur_row +=
static_cast<idx>(As[i].rows());
1029 cur_col +=
static_cast<idx>(As[i].cols());
1047 template <
typename Derived>
1048 dyn_mat<typename Derived::Scalar>
1049 dirsum(
const std::initializer_list<Derived>& As) {
1050 return dirsum(std::vector<Derived>(As));
1062 template <
typename Derived>
1063 dyn_mat<typename Derived::Scalar> dirsumpow(
const Eigen::MatrixBase<Derived>& A,
1065 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1071 throw exception::ZeroSize(
"qpp::dirsumpow()");
1075 throw exception::OutOfRange(
"qpp::dirsumpow()");
1078 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1094 template <
typename Derived>
1095 dyn_mat<typename Derived::Scalar> reshape(
const Eigen::MatrixBase<Derived>& A,
1097 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1099 idx Arows =
static_cast<idx>(rA.rows());
1100 idx Acols =
static_cast<idx>(rA.cols());
1106 throw exception::ZeroSize(
"qpp::reshape()");
1108 if (Arows * Acols != rows * cols)
1109 throw exception::DimsMismatchMatrix(
"qpp::reshape()");
1112 return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1113 const_cast<typename Derived::Scalar*
>(rA.data()), rows, cols);
1128 template <
typename Derived1,
typename Derived2>
1129 dyn_mat<typename Derived1::Scalar> comm(
const Eigen::MatrixBase<Derived1>& A,
1130 const Eigen::MatrixBase<Derived2>& B) {
1131 const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1132 const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1137 if (!std::is_same<
typename Derived1::Scalar,
1138 typename Derived2::Scalar>::value)
1139 throw exception::TypeMismatch(
"qpp::comm()");
1143 throw exception::ZeroSize(
"qpp::comm()");
1147 throw exception::MatrixNotSquare(
"qpp::comm()");
1150 if (rA.rows() != rB.rows())
1151 throw exception::DimsNotEqual(
"qpp::comm()");
1154 return rA * rB - rB * rA;
1169 template <
typename Derived1,
typename Derived2>
1170 dyn_mat<typename Derived1::Scalar>
1171 anticomm(
const Eigen::MatrixBase<Derived1>& A,
1172 const Eigen::MatrixBase<Derived2>& B) {
1173 const dyn_mat<typename Derived1::Scalar>& rA = A.derived();
1174 const dyn_mat<typename Derived2::Scalar>& rB = B.derived();
1179 if (!std::is_same<
typename Derived1::Scalar,
1180 typename Derived2::Scalar>::value)
1181 throw exception::TypeMismatch(
"qpp::anticomm()");
1185 throw exception::ZeroSize(
"qpp::anticomm()");
1189 throw exception::MatrixNotSquare(
"qpp::anticomm()");
1192 if (rA.rows() != rB.rows())
1193 throw exception::DimsNotEqual(
"qpp::anticomm()");
1196 return rA * rB + rB * rA;
1209 template <
typename Derived>
1210 dyn_mat<typename Derived::Scalar> prj(
const Eigen::MatrixBase<Derived>& A) {
1211 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1217 throw exception::ZeroSize(
"qpp::prj()");
1221 throw exception::MatrixNotCvector(
"qpp::prj()");
1224 double normA = norm(rA);
1226 return rA * adjoint(rA) / (normA * normA);
1228 return dyn_mat<typename Derived::Scalar>::Zero(rA.rows(), rA.rows());
1238 template <
typename Derived>
1239 dyn_mat<typename Derived::Scalar> grams(
const std::vector<Derived>& As) {
1244 throw exception::ZeroSize(
"qpp::grams()");
1246 for (
auto&& it : As)
1248 throw exception::ZeroSize(
"qpp::grams()");
1252 throw exception::MatrixNotCvector(
"qpp::grams()");
1255 for (
auto&& it : As)
1256 if (it.rows() != As[0].rows() || it.cols() != 1)
1257 throw exception::DimsNotEqual(
"qpp::grams()");
1260 dyn_mat<typename Derived::Scalar> cut =
1261 dyn_mat<typename Derived::Scalar>::Identity(As[0].rows(), As[0].rows());
1263 dyn_mat<typename Derived::Scalar> vi =
1264 dyn_mat<typename Derived::Scalar>::Zero(As[0].rows(), 1);
1266 std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1269 for (pos = 0; pos < As.size(); ++pos) {
1270 if (norm(As[pos]) >
eps)
1272 outvecs.push_back(As[pos]);
1278 for (
idx i = pos + 1; i < As.size(); ++i) {
1279 cut -= prj(outvecs[i - 1 - pos]);
1281 outvecs.push_back(vi);
1284 dyn_mat<typename Derived::Scalar> result(As[0].rows(), outvecs.size());
1287 for (
auto&& it : outvecs) {
1288 double normA = norm(it);
1291 result.col(cnt) = it / normA;
1296 return result.block(0, 0, As[0].rows(), cnt);
1307 template <
typename Derived>
1308 dyn_mat<typename Derived::Scalar>
1309 grams(
const std::initializer_list<Derived>& As) {
1310 return grams(std::vector<Derived>(As));
1320 template <
typename Derived>
1321 dyn_mat<typename Derived::Scalar> grams(
const Eigen::MatrixBase<Derived>& A) {
1322 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1327 throw exception::ZeroSize(
"qpp::grams()");
1330 std::vector<dyn_mat<typename Derived::Scalar>> input;
1332 for (
idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1333 input.push_back(rA.col(i));
1335 return grams<dyn_mat<typename Derived::Scalar>>(input);
1348 inline std::vector<idx>
n2multiidx(
idx n,
const std::vector<idx>& dims) {
1352 throw exception::DimsInvalid(
"qpp::n2multiidx()");
1354 if (n >= std::accumulate(std::begin(dims), std::end(dims),
1355 static_cast<idx>(1), std::multiplies<idx>()))
1356 throw exception::OutOfRange(
"qpp::n2multiidx()");
1363 return std::vector<idx>(result, result + dims.size());
1377 const std::vector<idx>& dims) {
1381 throw exception::DimsInvalid(
"qpp::multiidx2n()");
1383 for (
idx i = 0; i < dims.size(); ++i)
1384 if (midx[i] >= dims[i])
1385 throw exception::OutOfRange(
"qpp::multiidx2n()");
1403 inline ket mket(
const std::vector<idx>& mask,
const std::vector<idx>& dims) {
1404 idx N = mask.size();
1406 idx D = std::accumulate(std::begin(dims), std::end(dims),
1407 static_cast<idx>(1), std::multiplies<idx>());
1413 throw exception::ZeroSize(
"qpp::mket()");
1416 throw exception::DimsInvalid(
"qpp::mket()");
1418 if (mask.size() != dims.size())
1419 throw exception::SubsysMismatchDims(
"qpp::mket()");
1421 for (
idx i = 0; i < N; ++i)
1422 if (mask[i] >= dims[i])
1423 throw exception::SubsysMismatchDims(
"qpp::mket()");
1426 ket result = ket::Zero(D);
1445 inline ket mket(
const std::vector<idx>& mask,
idx d = 2) {
1446 idx N = mask.size();
1447 idx D =
static_cast<idx>(std::llround(std::pow(d, N)));
1453 throw exception::ZeroSize(
"qpp::mket()");
1457 throw exception::DimsInvalid(
"qpp::mket()");
1460 for (
idx i = 0; i < N; ++i)
1462 throw exception::SubsysMismatchDims(
"qpp::mket()");
1465 ket result = ket::Zero(D);
1466 std::vector<idx> dims(N, d);
1483 template <
char... Bits>
1484 ket operator"" _q() {
1485 constexpr
idx n =
sizeof...(Bits);
1486 constexpr
char bits[n + 1] = {Bits...,
'\0'};
1487 qpp::ket q = qpp::ket::Zero(std::pow(2, n));
1493 for (
idx i = 0; i < n; ++i) {
1494 if (bits[i] !=
'0' && bits[i] !=
'1')
1495 throw exception::OutOfRange(R
"xxx("qpp::operator "" _q())xxx"); 1499 pos = std::stoi(bits,
nullptr, 2);
1519 inline cmat mprj(
const std::vector<idx>& mask,
const std::vector<idx>& dims) {
1520 idx N = mask.size();
1522 idx D = std::accumulate(std::begin(dims), std::end(dims),
1523 static_cast<idx>(1), std::multiplies<idx>());
1529 throw exception::ZeroSize(
"qpp::mprj()");
1532 throw exception::DimsInvalid(
"qpp::mprj()");
1534 if (mask.size() != dims.size())
1535 throw exception::SubsysMismatchDims(
"qpp::mprj()");
1537 for (
idx i = 0; i < N; ++i)
1538 if (mask[i] >= dims[i])
1539 throw exception::SubsysMismatchDims(
"qpp::mprj()");
1542 cmat result = cmat::Zero(D, D);
1544 result(pos, pos) = 1;
1563 inline cmat mprj(
const std::vector<idx>& mask,
idx d = 2) {
1564 idx N = mask.size();
1565 idx D =
static_cast<idx>(std::llround(std::pow(d, N)));
1571 throw exception::ZeroSize(
"qpp::mprj()");
1575 throw exception::DimsInvalid(
"qpp::mprj()");
1578 for (
idx i = 0; i < N; ++i)
1580 throw exception::SubsysMismatchDims(
"qpp::mprj()");
1583 cmat result = cmat::Zero(D, D);
1584 std::vector<idx> dims(N, d);
1586 result(pos, pos) = 1;
1599 template <
typename InputIterator>
1600 std::vector<double> abssq(InputIterator first, InputIterator last) {
1601 std::vector<double> weights(std::distance(first, last));
1602 std::transform(first, last, std::begin(weights),
1603 [](
cplx z) ->
double {
return std::norm(z); });
1614 template <
typename Container>
1616 abssq(
const Container& c,
1617 typename std::enable_if<is_iterable<Container>::value>::type* =
nullptr)
1625 return abssq(std::begin(c), std::end(c));
1634 template <
typename Derived>
1635 std::vector<double> abssq(
const Eigen::MatrixBase<Derived>& A) {
1636 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1642 throw exception::ZeroSize(
"qpp::abssq()");
1645 return abssq(rA.data(), rA.data() + rA.size());
1656 template <
typename InputIterator>
1657 typename std::iterator_traits<InputIterator>::value_type
1658 sum(InputIterator first, InputIterator last) {
1659 using value_type =
typename std::iterator_traits<InputIterator>::value_type;
1661 return std::accumulate(first, last, static_cast<value_type>(0));
1671 template <
typename Container>
1672 typename Container::value_type
1673 sum(
const Container& c,
1674 typename std::enable_if<is_iterable<Container>::value>::type* =
nullptr) {
1675 return sum(std::begin(c), std::end(c));
1686 template <
typename InputIterator>
1687 typename std::iterator_traits<InputIterator>::value_type
1688 prod(InputIterator first, InputIterator last) {
1689 using value_type =
typename std::iterator_traits<InputIterator>::value_type;
1691 return std::accumulate(first, last, static_cast<value_type>(1),
1692 std::multiplies<value_type>());
1702 template <
typename Container>
1703 typename Container::value_type
1704 prod(
const Container& c,
1705 typename std::enable_if<is_iterable<Container>::value>::type* =
nullptr) {
1706 return prod(std::begin(c), std::end(c));
1721 template <
typename Derived>
1722 dyn_col_vect<typename Derived::Scalar>
1723 rho2pure(
const Eigen::MatrixBase<Derived>& A) {
1724 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1729 throw exception::ZeroSize(
"qpp::rho2pure()");
1732 throw exception::MatrixNotSquare(
"qpp::rho2pure()");
1735 dyn_col_vect<double> tmp_evals = hevals(rA);
1736 cmat tmp_evects = hevects(rA);
1737 dyn_col_vect<typename Derived::Scalar> result =
1738 dyn_col_vect<typename Derived::Scalar>::Zero(rA.rows());
1741 for (
idx k = 0; k < static_cast<idx>(rA.rows()); ++k) {
1742 if (std::abs(tmp_evals(k)) >
eps) {
1743 result = tmp_evects.col(k);
1759 template <
typename T>
1760 std::vector<T> complement(std::vector<T> subsys,
idx N) {
1763 if (N < subsys.size())
1764 throw exception::OutOfRange(
"qpp::complement()");
1767 std::vector<T> all(N);
1768 std::vector<T> subsys_bar(N - subsys.size());
1770 std::iota(std::begin(all), std::end(all), 0);
1771 std::sort(std::begin(subsys), std::end(subsys));
1772 std::set_difference(std::begin(all), std::end(all), std::begin(subsys),
1773 std::end(subsys), std::begin(subsys_bar));
1788 template <
typename Derived>
1789 std::vector<double> rho2bloch(
const Eigen::MatrixBase<Derived>& A) {
1790 const dyn_mat<typename Derived::Scalar>& rA = A.derived();
1796 throw exception::NotQubitMatrix(
"qpp::rho2bloch()");
1799 std::vector<double> result(3);
1800 cmat X(2, 2), Y(2, 2), Z(2, 2);
1802 Y << 0, -1_i, 1_i, 0;
1804 result[0] = std::real(trace(rA * X));
1805 result[1] = std::real(trace(rA * Y));
1806 result[2] = std::real(trace(rA * Z));
1819 inline cmat bloch2rho(
const std::vector<double>& r) {
1824 throw exception::CustomException(
"qpp::bloch2rho",
1825 "r is not a 3-dimensional vector!");
1828 cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1830 Y << 0, -1_i, 1_i, 0;
1834 return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
bool check_nonzero_size(const T &x) noexcept
Definition: util.h:123
constexpr idx maxn
Maximum number of allowed qubits/qudits (subsystems)
Definition: constants.h:73
dyn_mat< typename Derived1::Scalar > kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:273
bool check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:117
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:66
idx multiidx2n(const idx *const midx, idx numdims, const idx *const dims) noexcept
Definition: util.h:71
void n2multiidx(idx n, idx numdims, const idx *const dims, idx *result) noexcept
Definition: util.h:43
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:54
Quantum++ main namespace.
Definition: codes.h:35
std::complex< double > cplx
Complex number in double precision.
Definition: types.h:49
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:87
bool check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:99
std::size_t idx
Definition: experimental.h:43
dyn_mat< typename Derived1::Scalar > dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:317
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:64
bool check_dims(const std::vector< idx > &dims)
Definition: util.h:134
bool check_qubit_matrix(const Eigen::MatrixBase< Derived > &A) noexcept
Definition: util.h:236