40 template<
typename Derived>
42 const Eigen::MatrixBase<Derived>& A)
50 return rA.transpose();
60 template<
typename Derived>
62 const Eigen::MatrixBase<Derived>& A)
70 return rA.conjugate();
80 template<
typename Derived>
102 template<
typename Derived>
123 template<
typename Derived>
124 typename Derived::Scalar
trace(
const Eigen::MatrixBase<Derived>& A)
145 template<
typename Derived>
146 typename Derived::Scalar
det(
const Eigen::MatrixBase<Derived>& A)
157 return rA.determinant();
169 template<
typename Derived>
170 typename Derived::Scalar
logdet(
const Eigen::MatrixBase<Derived>& A)
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));
190 for (
idx i = 1; i < static_cast<idx>(rA.rows()); ++i)
191 result += std::log(U(i, i));
204 template<
typename Derived>
205 typename Derived::Scalar
sum(
const Eigen::MatrixBase<Derived>& A)
226 template<
typename Derived>
227 typename Derived::Scalar
prod(
const Eigen::MatrixBase<Derived>& A)
247 template<
typename Derived>
248 double norm(
const Eigen::MatrixBase<Derived>& A)
260 return (rA.template cast<cplx>()).norm();
271 template<
typename Derived>
272 std::pair<dyn_col_vect < cplx>,
cmat>
274 eig(
const Eigen::MatrixBase<Derived>& A)
289 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
291 return std::make_pair(es.eigenvalues(), es.eigenvectors());
301 template<
typename Derived>
317 return eig(rA).first;
327 template<
typename Derived>
343 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
345 return eig(rA).second;
356 template<
typename Derived>
357 std::pair<dyn_col_vect < double>,
cmat>
359 heig(
const Eigen::MatrixBase<Derived>& A)
374 Eigen::SelfAdjointEigenSolver<cmat> es(rA.template cast<cplx>());
376 return std::make_pair(es.eigenvalues(), es.eigenvectors());
386 template<
typename Derived>
402 return heig(rA).first;
412 template<
typename Derived>
428 return heig(rA).second;
440 template<
typename Derived>
441 std::tuple<cmat, dyn_col_vect < double>,
cmat>
443 svd(
const Eigen::MatrixBase<Derived>& A)
454 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
456 Eigen::DecompositionOptions::ComputeFullU |
457 Eigen::DecompositionOptions::ComputeFullV);
459 return std::make_tuple(sv.matrixU(), sv.singularValues(), sv.matrixV());
469 template<
typename Derived>
481 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>> sv(rA);
483 return sv.singularValues();
493 template<
typename Derived>
505 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
506 sv(rA, Eigen::DecompositionOptions::ComputeFullU);
518 template<
typename Derived>
530 Eigen::JacobiSVD<dyn_mat<typename Derived::Scalar>>
531 sv(rA, Eigen::DecompositionOptions::ComputeFullV);
545 template<
typename Derived>
561 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
564 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
567 cmat evalsdiag = evals.asDiagonal();
569 return evects * evalsdiag * evects.inverse();
578 template<
typename Derived>
594 return funm(rA, &std::sqrt);
603 template<
typename Derived>
628 template<
typename Derived>
644 return funm(rA, &std::exp);
653 template<
typename Derived>
669 return funm(rA, &std::log);
678 template<
typename Derived>
694 return funm(rA, &std::sin);
703 template<
typename Derived>
719 return funm(rA, &std::cos);
733 template<
typename Derived>
750 if (real(z) == 0 && imag(z) == 0)
751 return cmat::Identity(rA.rows(), rA.rows());
753 Eigen::ComplexEigenSolver<cmat> es(rA.template cast<cplx>());
756 for (
idx i = 0; i < static_cast<idx>(evals.rows()); ++i)
759 cmat evalsdiag = evals.asDiagonal();
761 return evects * evalsdiag * evects.inverse();
776 template<
typename Derived>
806 for (; n > 0; n /= 2)
809 result = (result * cA).eval();
810 cA = (cA * cA).eval();
824 template<
typename Derived>
825 double schatten(
const Eigen::MatrixBase<Derived>& A,
double p)
843 for (
idx i = 0; i < static_cast<idx>(sv.rows()); ++i)
844 result += std::pow(sv[i], p);
846 return std::pow(result, 1. / p);
859 template<
typename OutputScalar,
typename Derived>
862 const typename Derived::Scalar&))
876 #pragma omp parallel for collapse(2) 877 #endif // WITH_OPENMP_ 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));
914 template<
typename T,
typename ... Args>
929 template<
typename Derived>
943 for (
idx i = 1; i < As.size(); ++i)
945 result =
kron(result, As[i]);
963 template<
typename Derived>
965 const std::initializer_list<Derived>& As)
967 return kron(std::vector<Derived>(As));
979 template<
typename Derived>
996 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1013 template<
typename T>
1029 template<
typename T,
typename ... Args>
1044 template<
typename Derived>
1052 for (
auto&& it : As)
1057 idx total_rows = 0, total_cols = 0;
1058 for (
idx i = 0; i < As.size(); ++i)
1060 total_rows +=
static_cast<idx>(As[i].rows());
1061 total_cols +=
static_cast<idx>(As[i].cols());
1066 idx cur_row = 0, cur_col = 0;
1067 for (
idx i = 0; i < As.size(); ++i)
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());
1089 template<
typename Derived>
1091 const std::initializer_list<Derived>& As)
1093 return dirsum(std::vector<Derived>(As));
1105 template<
typename Derived>
1107 const Eigen::MatrixBase<Derived>& A,
idx n)
1122 std::vector<dyn_mat<typename Derived::Scalar>> As(n, rA);
1138 template<
typename Derived>
1144 idx Arows =
static_cast<idx>(rA.rows());
1145 idx Acols =
static_cast<idx>(rA.cols());
1153 if (Arows * Acols != rows * cols)
1157 return Eigen::Map<dyn_mat<typename Derived::Scalar >>(
1158 const_cast<typename Derived::Scalar*
>(rA.data()), rows, cols);
1173 template<
typename Derived1,
typename Derived2>
1175 const Eigen::MatrixBase<Derived2>& B)
1183 if (!std::is_same<
typename Derived1::Scalar,
1184 typename Derived2::Scalar>::value)
1197 if (rA.rows() != rB.rows())
1201 return rA * rB - rB * rA;
1216 template<
typename Derived1,
typename Derived2>
1218 const Eigen::MatrixBase<Derived1>& A,
1219 const Eigen::MatrixBase<Derived2>& B)
1227 if (!std::is_same<
typename Derived1::Scalar,
1228 typename Derived2::Scalar>::value)
1241 if (rA.rows() != rB.rows())
1245 return rA * rB + rB * rA;
1258 template<
typename Derived>
1274 double normA =
norm(rA);
1276 return rA *
adjoint(rA) / (normA * normA);
1279 ::Zero(rA.rows(), rA.rows());
1290 template<
typename Derived>
1299 for (
auto&& it : As)
1308 for (
auto&& it : As)
1309 if (it.rows() != As[0].rows() || it.cols() != 1)
1320 std::vector<dyn_mat<typename Derived::Scalar>> outvecs;
1323 for (pos = 0; pos < As.size(); ++pos)
1327 outvecs.push_back(As[pos]);
1333 for (
idx i = pos + 1; i < As.size(); ++i)
1335 cut -=
prj(outvecs[i - 1 - pos]);
1337 outvecs.push_back(vi);
1343 for (
auto&& it : outvecs)
1345 double normA =
norm(it);
1348 result.col(cnt) = it / normA;
1353 return result.block(0, 0, As[0].rows(), cnt);
1364 template<
typename Derived>
1366 const std::initializer_list<Derived>& As)
1368 return grams(std::vector<Derived>(As));
1378 template<
typename Derived>
1389 std::vector<dyn_mat<typename Derived::Scalar>> input;
1391 for (
idx i = 0; i < static_cast<idx>(rA.cols()); ++i)
1392 input.push_back(rA.col(i));
1394 return grams<dyn_mat<typename Derived::Scalar>>
1415 if (n >= std::accumulate(std::begin(dims), std::end(dims),
1416 static_cast<idx>(1), std::multiplies<idx>()))
1424 return std::vector<idx>(result, result + dims.size());
1438 const std::vector<idx>& dims)
1445 for (
idx i = 0; i < dims.size(); ++i)
1446 if (midx[i] >= dims[i])
1466 const std::vector<idx>& dims)
1468 idx N = mask.size();
1470 idx D = std::accumulate(std::begin(dims), std::end(dims),
1471 static_cast<idx>(1), std::multiplies<idx>());
1482 if (mask.size() != dims.size())
1485 for (
idx i = 0; i < N; ++i)
1486 if (mask[i] >= dims[i])
1490 ket result = ket::Zero(D);
1511 idx N = mask.size();
1512 idx D =
static_cast<idx>(std::llround(std::pow(d, N)));
1525 for (
idx i = 0; i < N; ++i)
1530 ket result = ket::Zero(D);
1531 std::vector<idx> dims(N, d);
1553 const std::vector<idx>& dims)
1555 idx N = mask.size();
1557 idx D = std::accumulate(std::begin(dims), std::end(dims),
1558 static_cast<idx>(1), std::multiplies<idx>());
1569 if (mask.size() != dims.size())
1572 for (
idx i = 0; i < N; ++i)
1573 if (mask[i] >= dims[i])
1577 cmat result = cmat::Zero(D, D);
1579 result(pos, pos) = 1;
1600 idx N = mask.size();
1601 idx D =
static_cast<idx>(std::llround(std::pow(d, N)));
1614 for (
idx i = 0; i < N; ++i)
1619 cmat result = cmat::Zero(D, D);
1620 std::vector<idx> dims(N, d);
1622 result(pos, pos) = 1;
1635 template<
typename InputIterator>
1636 std::vector<double>
abssq(InputIterator first, InputIterator last)
1638 std::vector<double> weights(std::distance(first, last));
1639 std::transform(first, last, std::begin(weights),
1640 [](
cplx z) ->
double 1654 template<
typename Container>
1655 std::vector<double>
abssq(
const Container& c,
1656 typename std::enable_if<
1666 return abssq(std::begin(c), std::end(c));
1675 template<
typename Derived>
1676 std::vector<double>
abssq(
const Eigen::MatrixBase<Derived>& A)
1687 return abssq(rA.data(), rA.data() + rA.size());
1698 template<
typename InputIterator>
1699 typename std::iterator_traits<InputIterator>::value_type
1700 sum(InputIterator first, InputIterator last)
1703 typename std::iterator_traits<InputIterator>::value_type;
1705 return std::accumulate(first, last, static_cast<value_type>(0));
1715 template<
typename Container>
1716 typename Container::value_type
1720 return sum(std::begin(c), std::end(c));
1731 template<
typename InputIterator>
1732 typename std::iterator_traits<InputIterator>::value_type
1733 prod(InputIterator first, InputIterator last)
1736 typename std::iterator_traits<InputIterator>::value_type;
1738 return std::accumulate(first, last, static_cast<value_type>(1),
1739 std::multiplies<value_type>());
1749 template<
typename Container>
1750 typename Container::value_type
1754 return prod(std::begin(c), std::end(c));
1769 template<
typename Derived>
1771 const Eigen::MatrixBase<Derived>& A)
1790 for (
idx k = 0; k < static_cast<idx>(rA.rows()); ++k)
1792 if (std::abs(tmp_evals(k)) >
eps)
1794 result = tmp_evects.col(k);
1810 template<
typename T>
1815 if (N < subsys.size())
1819 std::vector<T> all(N);
1820 std::vector<T> subsys_bar(N - subsys.size());
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));
1841 template<
typename Derived>
1842 std::vector<double>
rho2bloch(
const Eigen::MatrixBase<Derived>& A)
1853 std::vector<double> result(3);
1854 cmat X(2, 2), Y(2, 2), Z(2, 2);
1856 Y << 0, -1_i, 1_i, 0;
1858 result[0] = std::real(
trace(rA * X));
1859 result[1] = std::real(
trace(rA * Y));
1860 result[2] = std::real(
trace(rA * Z));
1880 "r is not a 3-dimensional vector!");
1883 cmat X(2, 2), Y(2, 2), Z(2, 2), Id2(2, 2);
1885 Y << 0, -1_i, 1_i, 0;
1889 return (Id2 + r[0] * X + r[1] * Y + r[2] * Z) / 2.;
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