y**b
发帖数: 10166 | 1 用openmp并行化一个模拟程序，发现：
计算只要运行足够步长，openmp每次运行给出的结果都不一样，
而串行运算的结果始终一样。这正常吗？
直觉上openmp由于每次运行都采用不同的计算顺序(各线程的
先后顺序是随机的)，从而可能改变误差的积累方式，一般怎么
处理这类问题？谢谢。 |
l********a
发帖数: 1154 | |
z*****a
发帖数: 3809 | 3 This is normal if the size of the differences are about the size of expected
floating point round-off errors. This is because in floating point
computations, (A+B)+C and A+(B+C) can produce different results.
One way to prevent such problems is to have all worker threads store their
partial results in separate locations and then combine the partial results
in a deterministic order after all the worker threads have finished.
【在 y**b 的大作中提到】
: 用openmp并行化一个模拟程序，发现：
: 计算只要运行足够步长，openmp每次运行给出的结果都不一样，
: 而串行运算的结果始终一样。这正常吗？
: 直觉上openmp由于每次运行都采用不同的计算顺序(各线程的
: 先后顺序是随机的)，从而可能改变误差的积累方式，一般怎么
: 处理这类问题？谢谢。
|
y**b
发帖数: 10166 | 4 多谢两位回复！
我这个问题本来就不存在一个deterministic order, 所以串行的结果一样并不能说明
串行的结果就比并行的结果更正确。还好并行的误差并不离谱，而且输出的结果明显
看出是舍入误差的累积，一位一位地放大。看intel的报告很多大型模型如天气模型会
因为舍入误差累积而达到量级差别，吓了一跳！
最近狠狠研究了几天并行和浮点consistency的问题，发现这个问题还不是一般的普遍，
各位做数值尤其是高精度高频率领域比如碰撞和爆炸可能会有些经验，欢迎指点。
我自己的思路是：
1. 先从gcc跳到icc试一下，毕竟icc提供80bits register extended precision和很多
浮点调控选项如-fp-model, 我试验了一下，发现icc不仅比gcc慢，也不能解决这个
consistency问题。(也说明gcc性能还是很牛的)。
2. GCC quad precision math lib, 这玩意比较新，_float128和sinq()等能提供大约
33位十进制进度，而double和long double只能提供15位精度。我的程序非常依赖这些
三角函数，据考证三角函数在角度极小的时候如果精度不够问题非常大、结果非常古怪，
希望用quadmath能有所改进。现在正在想办法wrap这些quadmath函数，使能选择
或者，否则对源码的改动太大、太不灵活。
这里要多说一下long double这个数据类型，用了128bit进行存储，结果跟double
一样只能提供15位十进制精度，这不是害人吗(当然，数据range能大一些)。测试了
一下，long double消耗的时间几乎恰恰就是double的两倍！
3. 其他
High Precision Arithmetic Library
gmp
mpc
mpfr
mpfr c++
这些库调用起来非常麻烦，暂时不考虑。
expected
【在 z*****a 的大作中提到】
: This is normal if the size of the differences are about the size of expected
: floating point round-off errors. This is because in floating point
: computations, (A+B)+C and A+(B+C) can produce different results.
: One way to prevent such problems is to have all worker threads store their
: partial results in separate locations and then combine the partial results
: in a deterministic order after all the worker threads have finished.
|
m*********t
发帖数: 527 | 5 很多非线性系统本来就是 chaotic 的，因为数值舍入产生的那些积累误差我觉得不代
结果是不对的吧。而且用来做数值的 model 本身也就是各种理想化, 即使那个模型能
够有一个解析解来给出预测也不一定和实际能对得上。。。 |
y**b
发帖数: 10166 | 6 更新一下，用了GCC quad-precision math libraray，初步结果显示openmp每次运行的
结果完全一致(在原来double输出精度的意义上），而原来double或long double在相同
运算下有明显偏差。
没有白折腾。感叹64位计算还没普及，128位计算已经颇有需求了，很多高精度库恐怕就
是例证。遗憾的是quadmath库目前很慢，我的计算显示大约慢30倍，够慢。
遍，
【在 y**b 的大作中提到】
: 多谢两位回复！
: 我这个问题本来就不存在一个deterministic order, 所以串行的结果一样并不能说明
: 串行的结果就比并行的结果更正确。还好并行的误差并不离谱，而且输出的结果明显
: 看出是舍入误差的累积，一位一位地放大。看intel的报告很多大型模型如天气模型会
: 因为舍入误差累积而达到量级差别，吓了一跳！
: 最近狠狠研究了几天并行和浮点consistency的问题，发现这个问题还不是一般的普遍，
: 各位做数值尤其是高精度高频率领域比如碰撞和爆炸可能会有些经验，欢迎指点。
: 我自己的思路是：
: 1. 先从gcc跳到icc试一下，毕竟icc提供80bits register extended precision和很多
: 浮点调控选项如-fp-model, 我试验了一下，发现icc不仅比gcc慢，也不能解决这个
|
y***d
发帖数: 2330 | 7 有两个方面建议考虑一下，
1，计算本身是不是稳定的？比如多体问题，舍入误差很快会导致不同的轨迹
2，有没有大量相减得小量的情况，或者极端的函数参数？比如有没有计算 sin(x)/x，
或者 cos(x)-1？
遍，
【在 y**b 的大作中提到】
: 多谢两位回复！
: 我这个问题本来就不存在一个deterministic order, 所以串行的结果一样并不能说明
: 串行的结果就比并行的结果更正确。还好并行的误差并不离谱，而且输出的结果明显
: 看出是舍入误差的累积，一位一位地放大。看intel的报告很多大型模型如天气模型会
: 因为舍入误差累积而达到量级差别，吓了一跳！
: 最近狠狠研究了几天并行和浮点consistency的问题，发现这个问题还不是一般的普遍，
: 各位做数值尤其是高精度高频率领域比如碰撞和爆炸可能会有些经验，欢迎指点。
: 我自己的思路是：
: 1. 先从gcc跳到icc试一下，毕竟icc提供80bits register extended precision和很多
: 浮点调控选项如-fp-model, 我试验了一下，发现icc不仅比gcc慢，也不能解决这个
|
m*********t
发帖数: 527 | 8 你需要问自己的是，这个精度到底多少 make sense。比如，你用来模拟的输入参数测
量精度是多少位有效数字？你的 model 本身用了多少 assumption ？ 单纯的追求浮点
数舍入误差我个人觉得没啥太大意义。
怕就
【在 y**b 的大作中提到】
: 更新一下，用了GCC quad-precision math libraray，初步结果显示openmp每次运行的
: 结果完全一致(在原来double输出精度的意义上），而原来double或long double在相同
: 运算下有明显偏差。
: 没有白折腾。感叹64位计算还没普及，128位计算已经颇有需求了，很多高精度库恐怕就
: 是例证。遗憾的是quadmath库目前很慢，我的计算显示大约慢30倍，够慢。
:
: 遍，
|
|
y**b
发帖数: 10166 | 9 嗯，这个比数值误差本身重要多了。我对模拟中的数据精度做了限制，
比如可测精度后面的数位全部舍弃，避免了不必要的舍入误差积累，
计算可重复性好多了。多谢。
【在 m*********t 的大作中提到】
: 你需要问自己的是，这个精度到底多少 make sense。比如，你用来模拟的输入参数测
: 量精度是多少位有效数字？你的 model 本身用了多少 assumption ？ 单纯的追求浮点
: 数舍入误差我个人觉得没啥太大意义。
:
: 怕就
|
y**b
发帖数: 10166 | 10 1. 多谢提醒，计算本身的稳定性是前提条件，不稳定的话进行比较就没有意义。
我的测试里面稳定和不稳定的情况都有(汗！), 这下问题清楚多了。因为非线性
模型的缘故，实际步长比用线性近似估算的步长要小不少才能稳定。除了从后处理
图像和能量角度检测稳定性之外，有没有其他什么办法？
2. sin/asin一类函数用得很频繁，这类函数在角度极小的情况下似乎比较敏感。
其他情况倒没有。
【在 y***d 的大作中提到】
: 有两个方面建议考虑一下，
: 1，计算本身是不是稳定的？比如多体问题，舍入误差很快会导致不同的轨迹
: 2，有没有大量相减得小量的情况，或者极端的函数参数？比如有没有计算 sin(x)/x，
: 或者 cos(x)-1？
:
: 遍，
|
y***d
发帖数: 2330 | 11 不知道你的系统是哪种类型，以及数值解法是什么？数值解不对劲的话，要看看高阶修
正。有的数值解法是不对的。
比如我以前做的牛顿力学的系统，数值解法应该符合时间反演不变的要求；因而直接用
x(t+dt) = x(t) + v(t)*dt, v(t+dt) = v(t) + a(x)*dt 肯定是不对的（不是不精确，
而是不对），应该采用的是 verlet 解法 http://en.wikipedia.org/wiki/Verlet_integration。
你可以用看看你的系统里面的守恒量是不是真的守恒。正确的解法加足够的精度，守恒
量只会有舍入造成的轻微噪音，而不正确的解法或者不够的精度，守恒量的数值会偏向
一方越跑越远。如果守恒量有问题，那你的解法可能有问题。
【在 y**b 的大作中提到】
: 1. 多谢提醒，计算本身的稳定性是前提条件，不稳定的话进行比较就没有意义。
: 我的测试里面稳定和不稳定的情况都有(汗！), 这下问题清楚多了。因为非线性
: 模型的缘故，实际步长比用线性近似估算的步长要小不少才能稳定。除了从后处理
: 图像和能量角度检测稳定性之外，有没有其他什么办法？
: 2. sin/asin一类函数用得很频繁，这类函数在角度极小的情况下似乎比较敏感。
: 其他情况倒没有。
|
y**b
发帖数: 10166 | 12 谢谢。忙的一直没空去研究这个verlet解法 :-) 在中心差分的基础上进一步提高
阶数，思路简洁漂亮。
结构动力学和地震工程在时域上普遍采用的是Newmark积分方法(也有用中心差分法)，
既没有人用Euler方法，也极少有人采用高阶Runge-kutta方法(代价太大，可能
只偶尔用于复杂弹塑性模型自身的积分)，多步法就更罕见了。离散元方面(同分子
动力学有类似之处)，verlet方法应该是比较合适的，目前普遍采用的还是中心差分法。
确，
【在 y***d 的大作中提到】
: 不知道你的系统是哪种类型，以及数值解法是什么？数值解不对劲的话，要看看高阶修
: 正。有的数值解法是不对的。
: 比如我以前做的牛顿力学的系统，数值解法应该符合时间反演不变的要求；因而直接用
: x(t+dt) = x(t) + v(t)*dt, v(t+dt) = v(t) + a(x)*dt 肯定是不对的（不是不精确，
: 而是不对），应该采用的是 verlet 解法 http://en.wikipedia.org/wiki/Verlet_integration。
: 你可以用看看你的系统里面的守恒量是不是真的守恒。正确的解法加足够的精度，守恒
: 量只会有舍入造成的轻微噪音，而不正确的解法或者不够的精度，守恒量的数值会偏向
: 一方越跑越远。如果守恒量有问题，那你的解法可能有问题。
|
r****t
发帖数: 10904 | 13 implicit/backward Euler 也没人用么？中心差分是不是就是 crank-nicolson? 那个
其实也有半步是 forward euler。。 verlet 一般用在哪些方面多？
Newmark，中心差分，verlet，高阶Runge-kutta。。。
法。
【在 y**b 的大作中提到】
: 谢谢。忙的一直没空去研究这个verlet解法 :-) 在中心差分的基础上进一步提高
: 阶数，思路简洁漂亮。
: 结构动力学和地震工程在时域上普遍采用的是Newmark积分方法(也有用中心差分法)，
: 既没有人用Euler方法，也极少有人采用高阶Runge-kutta方法(代价太大，可能
: 只偶尔用于复杂弹塑性模型自身的积分)，多步法就更罕见了。离散元方面(同分子
: 动力学有类似之处)，verlet方法应该是比较合适的，目前普遍采用的还是中心差分法。
:
: 确，
|
