@@ -32,9 +32,8 @@ namespace libmpdataxx
3232 static_assert (n_dims > 0 , " n_dims <= 0" );
3333 static_assert (n_tlev > 0 , " n_tlev <= 0" );
3434
35- // TODO: T_sumtype (perhaps worh using double even if summing floats?)
3635 std::unique_ptr<blitz::Array<real_t , 1 >> xtmtmp;
37- std::unique_ptr<blitz::Array<real_t , 1 >> sumtmp;
36+ std::unique_ptr<blitz::Array<double , 1 >> sumtmp;
3837
3938 protected:
4039
@@ -93,12 +92,12 @@ namespace libmpdataxx
9392 throw std::runtime_error (" number of subdomains greater than number of gridpoints" );
9493
9594 if (n_dims != 1 )
96- sumtmp.reset (new blitz::Array<real_t , 1 >(grid_size[0 ]));
95+ sumtmp.reset (new blitz::Array<double , 1 >(grid_size[0 ]));
9796 xtmtmp.reset (new blitz::Array<real_t , 1 >(size));
9897 }
9998
10099 // / @brief concurrency-aware summation of array elements
101- real_t sum (const arr_t &arr, const idx_t <n_dims> &ijk, const bool sum_khn)
100+ double sum (const arr_t &arr, const idx_t <n_dims> &ijk, const bool sum_khn)
102101 {
103102 // doing a two-step sum to reduce numerical error
104103 // and make parallel results reproducible
@@ -114,7 +113,7 @@ namespace libmpdataxx
114113 (*sumtmp)(c) = blitz::sum (arr (slice_idx));
115114 }
116115 barrier ();
117- real_t result;
116+ double result;
118117 if (sum_khn)
119118 result = blitz::kahan_sum (*sumtmp);
120119 else
@@ -124,7 +123,7 @@ namespace libmpdataxx
124123 }
125124
126125 // / @brief concurrency-aware summation of a (element-wise) product of two arrays
127- real_t sum (const arr_t &arr1, const arr_t &arr2, const idx_t <n_dims> &ijk, const bool sum_khn)
126+ double sum (const arr_t &arr1, const arr_t &arr2, const idx_t <n_dims> &ijk, const bool sum_khn)
128127 {
129128 // doing a two-step sum to reduce numerical error
130129 // and make parallel results reproducible
@@ -140,7 +139,7 @@ namespace libmpdataxx
140139 (*sumtmp)(c) = blitz::sum (arr1 (slice_idx) * arr2 (slice_idx));
141140 }
142141 barrier ();
143- real_t result;
142+ double result;
144143 if (sum_khn)
145144 result = blitz::kahan_sum (*sumtmp);
146145 else
0 commit comments