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>
105 template<
typename Derived>
126 template<
typename Derived>
127 typename Derived::Scalar
trace(
const Eigen::MatrixBase<Derived>& A)
148 template<
typename Derived>
149 typename Derived::Scalar
det(
const Eigen::MatrixBase<Derived>& A)
160 return rA.determinant();
172 template<
typename Derived>
173 typename Derived::Scalar
logdet(
const Eigen::MatrixBase<Derived>& A)
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));
194 for (
idx i = 1; i < static_cast<idx>(rA.rows()); ++i )
195 result += std::log(U(i, i));
208 template<
typename Derived>
209 typename Derived::Scalar
sum(
const Eigen::MatrixBase<Derived>& A)
230 template<
typename Derived>
231 typename Derived::Scalar
prod(
const Eigen::MatrixBase<Derived>& A)
251 template<
typename Derived>
252 double norm(
const Eigen::MatrixBase<Derived>& A)
264 return (rA.template cast<cplx>()).norm();
275 template<
typename Derived>
276 std::pair<dyn_col_vect < cplx>,
cmat>
278 eig(
const Eigen::MatrixBase<Derived>& A)
293 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
295 return std::make_pair(es.eigenvalues(), es.eigenvectors());
305 template<
typename Derived>
321 return eig(rA).first;
331 template<
typename Derived>
347 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
349 return eig(rA).second;
360 template<
typename Derived>
361 std::pair<dyn_col_vect < double>,
cmat>
363 heig(
const Eigen::MatrixBase<Derived>& A)
378 Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
380 return std::make_pair(es.eigenvalues(), es.eigenvectors());
390 template<
typename Derived>
406 return heig(rA).first;
416 template<
typename Derived>
433 return heig(rA).second;
445 template<
typename Derived>
446 std::tuple<cmat, dyn_col_vect < double>,
cmat>
448 svd(
const Eigen::MatrixBase<Derived>& A)
459 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
461 Eigen::DecompositionOptions::ComputeFullU |
462 Eigen::DecompositionOptions::ComputeFullV);
464 return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
474 template<
typename Derived>
486 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
488 return sv.singularValues();
498 template<
typename Derived>
510 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
511 sv(rA, Eigen::DecompositionOptions::ComputeFullU);
523 template<
typename Derived>
535 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
536 sv(rA, Eigen::DecompositionOptions::ComputeFullV);
550 template<
typename Derived>
566 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
569 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i )
572 cmat evalsdiag = evals.asDiagonal();
574 return evects * evalsdiag * evects.inverse();
583 template<
typename Derived>
599 return funm(rA, &std::sqrt);
608 template<
typename Derived>
633 template<
typename Derived>
649 return funm(rA, &std::exp);
658 template<
typename Derived>
674 return funm(rA, &std::log);
683 template<
typename Derived>
699 return funm(rA, &std::sin);
708 template<
typename Derived>
724 return funm(rA, &std::cos);
738 template<
typename Derived>
756 if ( real(z) == 0 && imag(z) == 0 )
757 return cmat::Identity(rA.rows(), rA.rows());
759 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
762 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i )
765 cmat evalsdiag = evals.asDiagonal();
767 return evects * evalsdiag * evects.inverse();
782 template<
typename Derived>
812 for ( ; n > 0; n /= 2 )
815 result = (result * cA).eval();
816 cA = (cA * cA).eval();
830 template<
typename Derived>
831 double schatten(
const Eigen::MatrixBase<Derived>& A,
double p)
860 template<
typename OutputScalar,
typename Derived>
863 const typename Derived::Scalar&))
877 #pragma omp parallel for collapse(2)
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));
915 template<
typename T,
typename ... Args>
930 template<
typename Derived>
935 if ( As.size() == 0 )
938 for (
auto&& it : As )
944 for (
idx i = 1; i < As.size(); ++i )
946 result =
kron(result, As[i]);
964 template<
typename Derived>
966 const std::initializer_list<Derived>& As)
968 return kron(std::vector<Derived>(As));
980 template<
typename Derived>
997 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1014 template<
typename T>
1030 template<
typename T,
typename ... Args>
1045 template<
typename Derived>
1050 if ( As.size() == 0 )
1053 for (
auto&& it : As )
1058 idx total_rows = 0, total_cols = 0;
1059 for (
idx i = 0; i < As.size(); ++i )
1061 total_rows +=
static_cast<idx>(As[i].rows());
1062 total_cols +=
static_cast<idx>(As[i].cols());
1067 idx cur_row = 0, cur_col = 0;
1068 for (
idx i = 0; i < As.size(); ++i )
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());
1090 template<
typename Derived>
1092 const std::initializer_list<Derived>& As)
1094 return dirsum(std::vector<Derived>(As));
1106 template<
typename Derived>
1108 const Eigen::MatrixBase<Derived>& A,
idx n)
1123 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1139 template<
typename Derived>
1145 idx Arows =
static_cast<idx>(rA.rows());
1146 idx Acols =
static_cast<idx>(rA.cols());
1154 if ( Arows * Acols != rows * cols )
1159 return Eigen::Map<dyn_mat<typename Derived::Scalar>>(
1160 const_cast<typename Derived::Scalar*
>(rA.data()), rows, cols);
1175 template<
typename Derived1,
typename Derived2>
1177 const Eigen::MatrixBase<Derived2>& B)
1185 if ( !std::is_same<
typename Derived1::Scalar,
1186 typename Derived2::Scalar>::value )
1199 if ( rA.rows() != rB.rows())
1203 return rA * rB - rB * rA;
1218 template<
typename Derived1,
typename Derived2>
1220 const Eigen::MatrixBase<Derived1>& A,
1221 const Eigen::MatrixBase<Derived2>& B)
1229 if ( !std::is_same<
typename Derived1::Scalar,
1230 typename Derived2::Scalar>::value )
1244 if ( rA.rows() != rB.rows())
1248 return rA * rB + rB * rA;
1261 template<
typename Derived>
1277 double normV =
norm(rV);
1279 return rV *
adjoint(rV) / (normV * normV);
1282 ::Zero(rV.rows(), rV.rows());
1293 template<
typename Derived>
1302 for (
auto&& it : Vs )
1311 for (
auto&& it : Vs )
1312 if ( it.rows() != Vs[0].rows() || it.cols() != 1 )
1323 std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1326 for ( pos = 0; pos < Vs.size(); ++pos )
1330 outvecs.push_back(Vs[pos]);
1336 for (
idx i = pos + 1; i < Vs.size(); ++i )
1338 cut -=
prj(outvecs[i - 1 - pos]);
1340 outvecs.push_back(vi);
1346 for (
auto&& it : outvecs )
1348 double normV =
norm(it);
1351 result.col(cnt) = it / normV;
1356 return result.block(0, 0, Vs[0].rows(), cnt);
1367 template<
typename Derived>
1369 const std::initializer_list<Derived>& Vs)
1371 return grams(std::vector<Derived>(Vs));
1381 template<
typename Derived>
1392 std::vector<dyn_mat<typename Derived::Scalar>> input;
1394 for (
idx i = 0; i < static_cast<idx>(rA.cols()); ++i )
1395 input.push_back(rA.col(i));
1397 return grams<dyn_mat<typename Derived::Scalar>>(input);
1417 if ( n >= std::accumulate(std::begin(dims), std::end(dims),
1418 static_cast<idx>(1), std::multiplies<idx>()))
1426 return std::vector<idx>(result, result + dims.size());
1440 const std::vector<idx>& dims)
1447 for (
idx i = 0; i < dims.size(); ++i )
1448 if ( midx[i] >= dims[i] )
1469 const std::vector<idx>& dims)
1471 idx n = mask.size();
1473 idx D = std::accumulate(std::begin(dims), std::end(dims),
1474 static_cast<idx>(1), std::multiplies<idx>());
1485 if ( mask.size() != dims.size())
1488 for (
idx i = 0; i < n; ++i )
1489 if ( mask[i] >= dims[i] )
1494 ket result = ket::Zero(D);
1515 idx n = mask.size();
1516 idx D =
static_cast<idx>(std::llround(std::pow(d, n)));
1529 for (
idx i = 0; i < n; ++i )
1535 ket result = ket::Zero(D);
1536 std::vector<idx> dims(n, d);
1558 const std::vector<idx>& dims)
1560 idx n = mask.size();
1562 idx D = std::accumulate(std::begin(dims), std::end(dims),
1563 static_cast<idx>(1), std::multiplies<idx>());
1574 if ( mask.size() != dims.size())
1577 for (
idx i = 0; i < n; ++i )
1578 if ( mask[i] >= dims[i] )
1583 cmat result = cmat::Zero(D, D);
1585 result(pos, pos) = 1;
1606 idx n = mask.size();
1607 idx D =
static_cast<idx>(std::llround(std::pow(d, n)));
1620 for (
idx i = 0; i < n; ++i )
1626 cmat result = cmat::Zero(D, D);
1627 std::vector<idx> dims(n, d);
1629 result(pos, pos) = 1;
1642 template<
typename InputIterator>
1643 std::vector<double>
abssq(InputIterator first, InputIterator last)
1645 std::vector<double> weights(std::distance(first, last));
1646 std::transform(first, last, std::begin(weights),
1647 [](
cplx z) ->
double
1661 template<
typename Container>
1662 std::vector<double>
abssq(
const Container& c,
1663 typename std::enable_if<
1673 return abssq(std::begin(c), std::end(c));
1682 template<
typename Derived>
1683 std::vector<double>
abssq(
const Eigen::MatrixBase<Derived>& A)
1694 return abssq(rA.data(), rA.data() + rA.size());
1705 template<
typename InputIterator>
1706 typename std::iterator_traits<InputIterator>::value_type
1707 sum(InputIterator first, InputIterator last)
1710 typename std::iterator_traits<InputIterator>::value_type;
1712 return std::accumulate(first, last, static_cast<value_type>(0));
1722 template<
typename Container>
1723 typename Container::value_type
1727 return sum(std::begin(c), std::end(c));
1738 template<
typename InputIterator>
1739 typename std::iterator_traits<InputIterator>::value_type
1740 prod(InputIterator first, InputIterator last)
1743 typename std::iterator_traits<InputIterator>::value_type;
1745 return std::accumulate(first, last, static_cast<value_type>(1),
1746 std::multiplies<value_type>());
1756 template<
typename Container>
1757 typename Container::value_type
1761 return prod(std::begin(c), std::end(c));
1776 template<
typename Derived>
1778 const Eigen::MatrixBase<Derived>& A)
1797 for (
idx k = 0; k < static_cast<idx>(rA.rows()); ++k )
1799 if ( std::abs(tmp_evals(k)) >
eps )
1801 result = tmp_evects.col(k);
1817 template<
typename T>
1822 if ( N < subsys.size())
1826 std::vector<T> all(N);
1827 std::vector<T> subsys_bar(N - subsys.size());
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));
1848 template<
typename Derived>
1850 const Eigen::MatrixBase<Derived>& A)
1861 std::vector<double> result(3);
1862 cmat X(2, 2), Y(2, 2), Z(2, 2);
1864 Y << 0, -1_i, 1_i, 0;
1866 result[0] = std::real(
trace(rA * X));
1867 result[1] = std::real(
trace(rA * Y));
1868 result[2] = std::real(
trace(rA * Z));
1886 if ( r.size() != 3 )
1887 throw Exception(
"qpp::bloch2rho",
"r is not a 3-dimensional vector!");
1890 cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1892 Y << 0, -1_i, 1_i, 0;
1896 return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
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
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
dyn_mat< OutputScalar > cwise(const Eigen::MatrixBase< Derived > &A, OutputScalar(*f)(const typename Derived::Scalar &))
Functor.
Definition: functions.h:861
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