| Message ID | orzhkiixqj.fsf@lxoliva.fsfla.org |
|---|---|
| State | New |
| Headers | show |
| Series | skip Cholesky decomposition in is>>n_mv_dist | expand |
On Fri, Aug 9, 2019 at 9:50 AM Alexandre Oliva <oliva@adacore.com> wrote: > 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. > > […] > > Tested on x86_64-linux-gnu. Ok to install? > Yes. Thanks.
On 09/08/19 10:20 +0200, Ulrich Drepper wrote: >On Fri, Aug 9, 2019 at 9:50 AM Alexandre Oliva <oliva@adacore.com> wrote: > >> 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. >> >> […] >> > > >> Tested on x86_64-linux-gnu. Ok to install? >> > >Yes. Thanks. If the operator>> is a friend it can just write straight to the array members of the param_type object: diff --git a/libstdc++-v3/include/ext/random.tcc b/libstdc++-v3/include/ext/random.tcc index 31dc33a2555..77abdd9a1de 100644 --- a/libstdc++-v3/include/ext/random.tcc +++ b/libstdc++-v3/include/ext/random.tcc @@ -700,18 +700,15 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION const typename __ios_base::fmtflags __flags = __is.flags(); __is.flags(__ios_base::dec | __ios_base::skipws); - std::array<_RealType, _Dimen> __mean; - for (auto& __it : __mean) + typename normal_mv_distribution<_Dimen, _RealType>::param_type __param; + for (auto& __it : __param._M_mean) __is >> __it; - std::array<_RealType, _Dimen * (_Dimen + 1) / 2> __varcov; - for (auto& __it : __varcov) + for (auto& __it : __param._M_t) __is >> __it; __is >> __x._M_nd; - __x.param(typename normal_mv_distribution<_Dimen, _RealType>:: - param_type(__mean.begin(), __mean.end(), - __varcov.begin(), __varcov.end())); + __x.param(__param); __is.flags(__flags); return __is; The default constructor for param_type() will pointlessly fill the arrays that are about to be overwritten though, so maybe this isn't an improvement.
diff --git a/libstdc++-v3/include/ext/random b/libstdc++-v3/include/ext/random index 41a2962c8f6e5..d5574e02ba02c 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<size_t _Dimen1, typename _RealType1, + typename _CharT, typename _Traits> + 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; }; diff --git a/libstdc++-v3/include/ext/random.tcc b/libstdc++-v3/include/ext/random.tcc index 31dc33a2555ed..a8a49a3a9fa6a 100644 --- a/libstdc++-v3/include/ext/random.tcc +++ b/libstdc++-v3/include/ext/random.tcc @@ -581,7 +581,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION __sum = *__varcovbegin++ - __sum; if (__builtin_expect(__sum <= _RealType(0), 0)) std::__throw_runtime_error(__N("normal_mv_distribution::" - "param_type::_M_init_full")); + "param_type::_M_init_lower")); *__w++ = std::sqrt(__sum); } } @@ -709,9 +709,11 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION __is >> __x._M_nd; + // The param_type temporary is built with a private constructor, + // to skip the Cholesky decomposition that would be performed + // otherwise. __x.param(typename normal_mv_distribution<_Dimen, _RealType>:: - param_type(__mean.begin(), __mean.end(), - __varcov.begin(), __varcov.end())); + param_type(__mean, __varcov)); __is.flags(__flags); return __is; diff --git a/libstdc++-v3/testsuite/ext/random/beta_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/beta_distribution/operators/serialize.cc index b05417156d191..a4925fc1c41be 100644 --- a/libstdc++-v3/testsuite/ext/random/beta_distribution/operators/serialize.cc +++ b/libstdc++-v3/testsuite/ext/random/beta_distribution/operators/serialize.cc @@ -23,6 +23,7 @@ #include <ext/random> #include <sstream> +#include <testsuite_hooks.h> void test01() @@ -35,6 +36,7 @@ test01() str << u; str >> v; + VERIFY( u == v ); } int main() diff --git a/libstdc++-v3/testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc index 9c2cc46ac1ce0..e9077b2c58d65 100644 --- a/libstdc++-v3/testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc +++ b/libstdc++-v3/testsuite/ext/random/hypergeometric_distribution/operators/serialize.cc @@ -38,6 +38,7 @@ test01() str << u; str >> v; + VERIFY( u == v ); } int diff --git a/libstdc++-v3/testsuite/ext/random/normal_mv_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/normal_mv_distribution/operators/serialize.cc index 8d83f9e6966d2..f5fbc42a686f0 100644 --- a/libstdc++-v3/testsuite/ext/random/normal_mv_distribution/operators/serialize.cc +++ b/libstdc++-v3/testsuite/ext/random/normal_mv_distribution/operators/serialize.cc @@ -23,6 +23,7 @@ #include <ext/random> #include <sstream> +#include <testsuite_hooks.h> void test01() @@ -35,6 +36,7 @@ test01() str << u; str >> v; + VERIFY( u == v ); } int main() diff --git a/libstdc++-v3/testsuite/ext/random/triangular_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/triangular_distribution/operators/serialize.cc index cf17fea8b03ff..75e16cf0437cb 100644 --- a/libstdc++-v3/testsuite/ext/random/triangular_distribution/operators/serialize.cc +++ b/libstdc++-v3/testsuite/ext/random/triangular_distribution/operators/serialize.cc @@ -23,6 +23,7 @@ #include <ext/random> #include <sstream> +#include <testsuite_hooks.h> void test01() @@ -35,6 +36,7 @@ test01() str << u; str >> v; + VERIFY( u == v ); } int main() diff --git a/libstdc++-v3/testsuite/ext/random/von_mises_distribution/operators/serialize.cc b/libstdc++-v3/testsuite/ext/random/von_mises_distribution/operators/serialize.cc index f3d7912e314ba..b32a31dee6421 100644 --- a/libstdc++-v3/testsuite/ext/random/von_mises_distribution/operators/serialize.cc +++ b/libstdc++-v3/testsuite/ext/random/von_mises_distribution/operators/serialize.cc @@ -23,6 +23,7 @@ #include <ext/random> #include <sstream> +#include <testsuite_hooks.h> void test01() @@ -35,6 +36,7 @@ test01() str << u; str >> v; + VERIFY( u == v ); } int main()