43 template<
typename Derived>
45 const Eigen::MatrixBase<Derived>& A)
53 return rA.transpose();
63 template<
typename Derived>
65 const Eigen::MatrixBase<Derived>& A)
73 return rA.conjugate();
83 template<
typename Derived>
102 template<
typename Derived>
120 template<
typename Derived>
121 typename Derived::Scalar
trace(
const Eigen::MatrixBase<Derived>& A)
139 template<
typename Derived>
140 typename Derived::Scalar
det(
const Eigen::MatrixBase<Derived>& A)
148 return rA.determinant();
160 template<
typename Derived>
161 typename Derived::Scalar
logdet(
const Eigen::MatrixBase<Derived>& A)
174 Eigen::PartialPivLU<dyn_mat<typename Derived::Scalar>> lu(rA);
176 lu.matrixLU().template triangularView<Eigen::Upper>();
177 typename Derived::Scalar result = std::log(U(0, 0));
179 for (
idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
180 result += std::log(U(i, i));
193 template<
typename Derived>
194 typename Derived::Scalar
sum(
const Eigen::MatrixBase<Derived>& A)
212 template<
typename Derived>
213 typename Derived::Scalar
prod(
const Eigen::MatrixBase<Derived>& A)
230 template<
typename Derived>
231 double norm(
const Eigen::MatrixBase<Derived>& A)
240 return (rA.template cast<cplx>()).norm();
251 template<
typename Derived>
252 std::pair<dyn_col_vect<cplx>,
cmat>
eig(
const Eigen::MatrixBase<Derived>& A)
264 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
266 return std::make_pair(es.eigenvalues(), es.eigenvectors());
276 template<
typename Derived>
289 return eig(rA).first;
299 template<
typename Derived>
312 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
314 return eig(rA).second;
325 template<
typename Derived>
327 const Eigen::MatrixBase<Derived>& A)
339 Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
341 return std::make_pair(es.eigenvalues(), es.eigenvectors());
351 template<
typename Derived>
364 return heig(rA).first;
374 template<
typename Derived>
388 return heig(rA).second;
400 template<
typename Derived>
401 std::tuple<cmat, dyn_col_vect<double>,
cmat>
402 svd(
const Eigen::MatrixBase<Derived>& A)
410 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
412 Eigen::DecompositionOptions::ComputeFullU |
413 Eigen::DecompositionOptions::ComputeFullV);
415 return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
425 template<
typename Derived>
434 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
436 return sv.singularValues();
446 template<
typename Derived>
455 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
456 sv(rA, Eigen::DecompositionOptions::ComputeFullU);
468 template<
typename Derived>
477 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
478 sv(rA, Eigen::DecompositionOptions::ComputeFullV);
492 template<
typename Derived>
505 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
508 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
511 cmat evalsdiag = evals.asDiagonal();
513 return evects * evalsdiag * evects.inverse();
522 template<
typename Derived>
535 return funm(rA, &std::sqrt);
544 template<
typename Derived>
566 template<
typename Derived>
579 return funm(rA, &std::exp);
588 template<
typename Derived>
601 return funm(rA, &std::log);
610 template<
typename Derived>
623 return funm(rA, &std::sin);
632 template<
typename Derived>
645 return funm(rA, &std::cos);
659 template<
typename Derived>
674 if (real(z) == 0 && imag(z) == 0)
675 return cmat::Identity(rA.rows(), rA.rows());
677 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
680 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
683 cmat evalsdiag = evals.asDiagonal();
685 return evects * evalsdiag * evects.inverse();
700 template<
typename Derived>
716 cA.rows(), cA.rows());
719 for (; n > 0; n /= 2)
722 result = (result * cA).eval();
723 cA = (cA * cA).eval();
737 template<
typename Derived>
738 double schatten(
const Eigen::MatrixBase<Derived>& A,
double p)
764 template<
typename OutputScalar,
typename Derived>
766 OutputScalar (* f)(
const typename Derived::Scalar&))
776 #pragma omp parallel for collapse(2)
778 for (
idx j = 0; j < static_cast<idx>(rA.cols()); ++j)
779 for (
idx i = 0; i < static_cast<idx>(rA.rows()); ++i)
780 result(i, j) = (*f)(rA(i, j));
813 template<
typename T,
typename ... Args>
828 template<
typename Derived>
839 for (
idx i = 1; i < As.size(); ++i)
841 result =
kron(result, As[i]);
859 template<
typename Derived>
861 const std::initializer_list<Derived>& As)
863 return kron(std::vector<Derived>(As));
875 template<
typename Derived>
889 std::vector <dyn_mat<typename Derived::Scalar>> As(n, rA);
922 template<
typename T,
typename ... Args>
937 template<
typename Derived>
947 idx total_rows = 0, total_cols = 0;
948 for (
idx i = 0; i < As.size(); ++i)
950 total_rows +=
static_cast<idx>(As[i].rows());
951 total_cols +=
static_cast<idx>(As[i].cols());
956 idx cur_row = 0, cur_col = 0;
957 for (
idx i = 0; i < As.size(); ++i)
959 result.block(cur_row, cur_col, As[i].rows(), As[i].cols()) = As[i];
960 cur_row +=
static_cast<idx>(As[i].rows());
961 cur_col +=
static_cast<idx>(As[i].cols());
979 template<
typename Derived>
981 const std::initializer_list<Derived>& As)
983 return dirsum(std::vector<Derived>(As));
995 template<
typename Derived>
997 const Eigen::MatrixBase<Derived>& A,
idx n)
1009 std::vector <dyn_mat<typename Derived::Scalar>> As(n, rA);
1025 template<
typename Derived>
1031 idx Arows =
static_cast<idx>(rA.rows());
1032 idx Acols =
static_cast<idx>(rA.cols());
1038 if (Arows * Acols != rows * cols)
1042 return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1043 const_cast<typename Derived::Scalar*
>(rA.data()), rows, cols);
1058 template<
typename Derived1,
typename Derived2>
1060 const Eigen::MatrixBase<Derived2>& B)
1068 if (!std::is_same<
typename Derived1::Scalar,
1069 typename Derived2::Scalar>::value)
1082 if (rA.rows() != rB.rows())
1085 return rA * rB - rB * rA;
1100 template<
typename Derived1,
typename Derived2>
1102 const Eigen::MatrixBase<Derived1>& A,
1103 const Eigen::MatrixBase<Derived2>& B)
1111 if (!std::is_same<
typename Derived1::Scalar,
1112 typename Derived2::Scalar>::value)
1126 if (rA.rows() != rB.rows())
1129 return rA * rB + rB * rA;
1142 template<
typename Derived>
1155 double normV =
norm(rV);
1157 return rV *
adjoint(rV) / (normV * normV);
1160 ::Zero(rV.rows(), rV.rows());
1171 template<
typename Derived>
1178 for (
auto&& it : Vs)
1187 for (
auto&& it : Vs)
1188 if (it.rows() != Vs[0].rows() || it.cols() != 1)
1198 std::vector <dyn_mat<typename Derived::Scalar>> outvecs;
1201 for (pos = 0; pos < Vs.size(); ++pos)
1205 outvecs.push_back(Vs[pos]);
1211 for (
idx i = pos + 1; i < Vs.size(); ++i)
1213 cut -=
prj(outvecs[i - 1 - pos]);
1215 outvecs.push_back(vi);
1221 for (
auto&& it : outvecs)
1223 double normV =
norm(it);
1226 result.col(cnt) = it / normV;
1231 return result.block(0, 0, Vs[0].rows(), cnt);
1242 template<
typename Derived>
1244 const std::initializer_list<Derived>& Vs)
1246 return grams(std::vector<Derived>(Vs));
1256 template<
typename Derived>
1264 std::vector <dyn_mat<typename Derived::Scalar>> input;
1266 for (
idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1267 input.push_back(rA.col(i));
1269 return grams<dyn_mat<typename Derived::Scalar>>(input);
1287 if (n >= std::accumulate(std::begin(dims), std::end(dims),
1288 static_cast<idx>(1), std::multiplies<idx>()))
1295 return std::vector<idx>(result, result + dims.size());
1309 const std::vector <idx>& dims)
1314 for (
idx i = 0; i < dims.size(); ++i)
1315 if (midx[i] >= dims[i])
1335 const std::vector <idx>& dims)
1337 idx n = mask.size();
1339 idx D = std::accumulate(std::begin(dims), std::end(dims),
1340 static_cast<idx>(1), std::multiplies<idx>());
1349 if (mask.size() != dims.size())
1352 for (
idx i = 0; i < n; ++i)
1353 if (mask[i] >= dims[i])
1357 ket result = ket::Zero(D);
1378 idx n = mask.size();
1379 idx D =
static_cast<idx>(std::llround(std::pow(d, n)));
1388 for (
idx i = 0; i < n; ++i)
1393 ket result = ket::Zero(D);
1394 std::vector <idx> dims(n, d);
1416 const std::vector <idx>& dims)
1418 idx n = mask.size();
1420 idx D = std::accumulate(std::begin(dims), std::end(dims),
1421 static_cast<idx>(1), std::multiplies<idx>());
1430 if (mask.size() != dims.size())
1433 for (
idx i = 0; i < n; ++i)
1434 if (mask[i] >= dims[i])
1438 cmat result = cmat::Zero(D, D);
1440 result(pos, pos) = 1;
1461 idx n = mask.size();
1462 idx D =
static_cast<idx>(std::llround(std::pow(d, n)));
1471 for (
idx i = 0; i < n; ++i)
1476 cmat result = cmat::Zero(D, D);
1477 std::vector <idx> dims(n, d);
1479 result(pos, pos) = 1;
1491 template<
typename InputIterator>
1492 std::vector <double>
abssq(InputIterator first, InputIterator last)
1494 std::vector <double> weights(std::distance(first, last));
1495 std::transform(first, last, std::begin(weights),
1496 [](
cplx z) ->
double
1510 template<
typename Derived>
1511 std::vector <double>
abssq(
const Eigen::MatrixBase<Derived>& V)
1523 return abssq(rV.data(), rV.data() + rV.rows());
1534 template<
typename InputIterator>
1535 typename std::iterator_traits<InputIterator>::value_type
1536 sum(InputIterator first, InputIterator last)
1539 typename std::iterator_traits<InputIterator>::value_type;
1541 return std::accumulate(first, last, static_cast<value_type>(0));
1551 template<
typename Container>
1552 typename Container::value_type
1555 using value_type =
typename Container::value_type;
1557 return std::accumulate(std::begin(c), std::end(c),
1558 static_cast<value_type>(0));
1569 template<
typename InputIterator>
1570 typename std::iterator_traits<InputIterator>::value_type
1571 prod(InputIterator first, InputIterator last)
1574 typename std::iterator_traits<InputIterator>::value_type;
1576 return std::accumulate(first, last, static_cast<value_type>(1),
1577 std::multiplies<value_type>());
1588 template<
typename Container>
1589 typename Container::value_type
1592 using value_type =
typename Container::value_type;
1594 return std::accumulate(std::begin(c), std::end(c),
1595 static_cast<value_type>(1),
1596 std::multiplies<value_type>());
1611 template<
typename Derived>
1613 const Eigen::MatrixBase<Derived>& A)
1632 for (
idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1634 if (std::abs(tmp_evals(k)) >
eps)
1636 result = tmp_evects.col(k);
1652 template<
typename T>
1655 if (N < subsys.size())
1658 std::vector <T> all(N);
1659 std::vector <T> subsys_bar(N - subsys.size());
1661 std::iota(std::begin(all), std::end(all), 0);
1662 std::sort(std::begin(subsys), std::end(subsys));
1663 std::set_difference(std::begin(all), std::end(all),
1664 std::begin(subsys), std::end(subsys),
1665 std::begin(subsys_bar));
1680 template<
typename Derived>
1682 const Eigen::MatrixBase<Derived>& A)
1690 std::vector <double> result(3);
1691 cmat X(2, 2), Y(2, 2), Z(2, 2);
1693 Y << 0, -1_i, 1_i, 0;
1695 result[0] = std::real(
trace(rA * X));
1696 result[1] = std::real(
trace(rA * Y));
1697 result[2] = std::real(
trace(rA * Z));
1714 throw Exception(
"qpp::bloch2rho",
"r is not a 3-dimensional vector!");
1716 cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1718 Y << 0, -1_i, 1_i, 0;
1722 return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
bool _check_cvector(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:104
std::pair< dyn_col_vect< double >, cmat > heig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition of Hermitian expression.
Definition: functions.h:326
Derived::Scalar logdet(const Eigen::MatrixBase< Derived > &A)
Logarithm of the determinant.
Definition: functions.h:161
constexpr idx maxn
Maximum number of allowed qu(d)its (subsystems)
Definition: constants.h:82
cmat absm(const Eigen::MatrixBase< Derived > &A)
Matrix absolut value.
Definition: functions.h:545
dyn_mat< typename Derived1::Scalar > _kron2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:265
cmat expm(const Eigen::MatrixBase< Derived > &A)
Matrix exponential.
Definition: functions.h:567
Derived::Scalar prod(const Eigen::MatrixBase< Derived > &A)
Element-wise product of A.
Definition: functions.h:213
constexpr double eps
Used to decide whether a number or expression in double precision is zero or not. ...
Definition: constants.h:75
dyn_mat< typename T::Scalar > dirsum(const T &head)
Direct sum.
Definition: functions.h:907
Eigen::Matrix< Scalar, Eigen::Dynamic, Eigen::Dynamic > dyn_mat
Dynamic Eigen matrix over the field specified by Scalar.
Definition: types.h:84
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:1681
Eigen::VectorXcd ket
Complex (double precision) dynamic Eigen column vector.
Definition: types.h:57
dyn_mat< typename Derived1::Scalar > anticomm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Anti-commutator.
Definition: functions.h:1101
Quantum++ main namespace.
Definition: codes.h:30
std::pair< dyn_col_vect< cplx >, cmat > eig(const Eigen::MatrixBase< Derived > &A)
Full eigen decomposition.
Definition: functions.h:252
double norm(const Eigen::MatrixBase< Derived > &A)
Frobenius norm.
Definition: functions.h:231
cmat spectralpowm(const Eigen::MatrixBase< Derived > &A, const cplx z)
Matrix power.
Definition: functions.h:660
bool _check_square_mat(const Eigen::MatrixBase< Derived > &A)
Definition: util.h:77
dyn_col_vect< double > hevals(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvalues.
Definition: functions.h:352
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:765
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:1612
dyn_mat< typename Derived::Scalar > dirsumpow(const Eigen::MatrixBase< Derived > &A, idx n)
Direct sum power.
Definition: functions.h:996
std::vector< T > complement(std::vector< T > subsys, idx N)
Constructs the complement of a subsystem vector.
Definition: functions.h:1653
cmat evects(const Eigen::MatrixBase< Derived > &A)
Eigenvectors.
Definition: functions.h:300
cmat funm(const Eigen::MatrixBase< Derived > &A, cplx(*f)(const cplx &))
Functional calculus f(A)
Definition: functions.h:493
double schatten(const Eigen::MatrixBase< Derived > &A, double p)
Schatten matrix norm.
Definition: functions.h:738
cmat svdU(const Eigen::MatrixBase< Derived > &A)
Left singular vectors.
Definition: functions.h:447
cmat bloch2rho(const std::vector< double > &r)
Computes the density matrix corresponding to the 3-dimensional real Bloch vector r.
Definition: functions.h:1710
std::vector< idx > n2multiidx(idx n, const std::vector< idx > &dims)
Non-negative integer index to multi-index.
Definition: functions.h:1282
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:1308
dyn_col_vect< double > svals(const Eigen::MatrixBase< Derived > &A)
Singular values.
Definition: functions.h:426
dyn_mat< typename T::Scalar > kron(const T &head)
Kronecker product.
Definition: functions.h:798
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:52
idx _multiidx2n(const idx *midx, idx numdims, const idx *dims) noexcept
Definition: util.h:54
cmat mprj(const std::vector< idx > &mask, const std::vector< idx > &dims)
Projector onto multi-partite qudit ket.
Definition: functions.h:1415
Generates custom exceptions, used when validating function parameters.
Definition: exception.h:39
Derived::Scalar det(const Eigen::MatrixBase< Derived > &A)
Determinant.
Definition: functions.h:140
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:876
dyn_mat< typename Derived::Scalar > powm(const Eigen::MatrixBase< Derived > &A, idx n)
Matrix power.
Definition: functions.h:701
dyn_mat< typename Derived::Scalar > prj(const Eigen::MatrixBase< Derived > &V)
Projector.
Definition: functions.h:1143
dyn_col_vect< cplx > evals(const Eigen::MatrixBase< Derived > &A)
Eigenvalues.
Definition: functions.h:277
cmat hevects(const Eigen::MatrixBase< Derived > &A)
Hermitian eigenvectors.
Definition: functions.h:375
dyn_mat< typename Derived::Scalar > conjugate(const Eigen::MatrixBase< Derived > &A)
Complex conjugate.
Definition: functions.h:64
bool _check_nonzero_size(const T &x) noexcept
Definition: util.h:113
constexpr double infty
Used to denote infinity in double precision.
Definition: constants.h:96
dyn_mat< typename Derived1::Scalar > _dirsum2(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Definition: util.h:305
Derived::Scalar trace(const Eigen::MatrixBase< Derived > &A)
Trace.
Definition: functions.h:121
cmat svdV(const Eigen::MatrixBase< Derived > &A)
Right singular vectors.
Definition: functions.h:469
cmat cosm(const Eigen::MatrixBase< Derived > &A)
Matrix cos.
Definition: functions.h:633
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 > dyn_col_vect
Dynamic Eigen column vector over the field specified by Scalar.
Definition: types.h:96
cmat sinm(const Eigen::MatrixBase< Derived > &A)
Matrix sin.
Definition: functions.h:611
std::vector< double > abssq(InputIterator first, InputIterator last)
Computes the absolute values squared of a range of complex numbers.
Definition: functions.h:1492
dyn_mat< typename Derived1::Scalar > comm(const Eigen::MatrixBase< Derived1 > &A, const Eigen::MatrixBase< Derived2 > &B)
Commutator.
Definition: functions.h:1059
Eigen::MatrixXcd cmat
Complex (double precision) dynamic Eigen matrix.
Definition: types.h:67
ket mket(const std::vector< idx > &mask, const std::vector< idx > &dims)
Multi-partite qudit ket.
Definition: functions.h:1334
void _n2multiidx(idx n, idx numdims, const idx *dims, idx *result) noexcept
Definition: util.h:41
dyn_mat< typename Derived::Scalar > inverse(const Eigen::MatrixBase< Derived > &A)
Inverse.
Definition: functions.h:103
cmat logm(const Eigen::MatrixBase< Derived > &A)
Matrix logarithm.
Definition: functions.h:589
cmat sqrtm(const Eigen::MatrixBase< Derived > &A)
Matrix square root.
Definition: functions.h:523
std::size_t idx
Non-negative integer index.
Definition: types.h:36
bool _check_dims(const std::vector< idx > &dims)
Definition: util.h:119
Derived::Scalar sum(const Eigen::MatrixBase< Derived > &A)
Element-wise sum of A.
Definition: functions.h:194
dyn_mat< typename Derived::Scalar > reshape(const Eigen::MatrixBase< Derived > &A, idx rows, idx cols)
Reshape.
Definition: functions.h:1026
std::tuple< cmat, dyn_col_vect< double >, cmat > svd(const Eigen::MatrixBase< Derived > &A)
Full singular value decomposition.
Definition: functions.h:402
dyn_mat< typename Derived::Scalar > grams(const std::vector< Derived > &Vs)
Gram-Schmidt orthogonalization.
Definition: functions.h:1172