From c787deb0124b667802d8519bc285894bb6d771d7 Mon Sep 17 00:00:00 2001 From: Alexandre Oliva Date: Fri, 9 Aug 2019 09:20:58 +0000 Subject: skip Cholesky decomposition in is>>n_mv_dist normal_mv_distribution maintains the variance-covariance matrix param in Cholesky-decomposed form. Existing param_type constructors, when taking a full or lower-triangle varcov matrix, perform Cholesky decomposition to convert it to the internal representation. This internal representation is visible both in the varcov() result, and in the streamed-out representation of a normal_mv_distribution object. The problem is that when that representation is streamed back in, the read-back decomposed varcov matrix is used as a lower-triangle non-decomposed varcov matrix, and it undergoes Cholesky decomposition again. So, each cycle of stream-out/stream-in changes the varcov matrix to its "square root", instead of restoring the original params. This patch includes Corentin's changes that introduce verification in testsuite/ext/random/normal_mv_distribution/operators/serialize.cc and other similar tests that the object read back in compares equal to the written-out object: the modified tests pass only if (u == v). This patch also fixes the error exposed by his change, introducing an alternate private constructor for param_type, used only by operator>>. for libstdc++-v3/ChangeLog * include/ext/random (normal_mv_distribution::param_type::param_type): New private ctor taking a decomposed varcov matrix, for use by... (operator>>): ... this, befriended. * include/ext/random.tcc (operator>>): Use it. (normal_mv_distribution::param_type::_M_init_lower): Adjust member function name in exception message. for libstdc++-v3/ChangeLog from Corentin Gay * testsuite/ext/random/beta_distribution/operators/serialize.cc, testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc, testsuite/ext/random/normal_mv_distribution/operators/serialize.cc, testsuite/ext/random/triangular_distribution/operators/serialize.cc, testsuite/ext/random/von_mises_distribution/operators/serialize.cc: Add call to `VERIFY`. From-SVN: r274233 --- libstdc++-v3/include/ext/random | 15 +++++++++++++++ 1 file changed, 15 insertions(+) (limited to 'libstdc++-v3/include/ext/random') diff --git a/libstdc++-v3/include/ext/random b/libstdc++-v3/include/ext/random index 41a2962..d5574e0 100644 --- a/libstdc++-v3/include/ext/random +++ b/libstdc++-v3/include/ext/random @@ -752,6 +752,21 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION _InputIterator2 __varbegin, _InputIterator2 __varend); + // param_type constructors apply Cholesky decomposition to the + // varcov matrix in _M_init_full and _M_init_lower, but the + // varcov matrix output ot a stream is already decomposed, so + // we need means to restore it as-is when reading it back in. + template + friend std::basic_istream<_CharT, _Traits>& + operator>>(std::basic_istream<_CharT, _Traits>& __is, + __gnu_cxx::normal_mv_distribution<_Dimen1, _RealType1>& + __x); + param_type(std::array<_RealType, _Dimen> const &__mean, + std::array<_RealType, _M_t_size> const &__varcov) + : _M_mean (__mean), _M_t (__varcov) + {} + std::array<_RealType, _Dimen> _M_mean; std::array<_RealType, _M_t_size> _M_t; }; -- cgit v1.1