发个Eigen的测试兼与Fortran的对比

adda
adda 2016-08-16 字数 2624

做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往C++转,底层数组使用Eigen。

由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,用到的数组大小/成员函数等都是静态的,多态也静态实现。求解变量5个的话,大量用到了Array3d和Array5d这类小的数组。

这么做刚开始还可以,但程序大了以后开发效率太低。

一是对编程的人要求高,组内其他人改起来很慢;

二是各种自动化分析工具基本不可用,重构起来太麻烦。

此外,编译时间太长了,随便修改一行,编译的时间都够泡杯茶喝了。

为了提高开发效率,就有了传统的虚函数和动态数组的C++代码(代码3),大量用到了ArrayXd和各种虚函数。

以代码举例的话,为实现一个x->y的函数(x,y为矢量),三种代码的写法分别是:

代码1:

subroutine f(x, y)

real, intent(in)  :: x(:)

real, intent(out) :: y(:)

代码2:

void f(const Ref<const ArrayNd>& x, Ref<ArrayNd> y);    // N是编译期常量

代码3:

ArrayXd f(const Ref<const ArrayXd>& x);

测试的算例是一个简单的翼型绕流,网格及控制参数完全一致。使用显式时间推进,不涉及线性方程组求解。时间推进100步,耗时:

代码1. Fortran,动态数组/虚函数   IFC 12.0,/O3    4.078s

代码2. C++,静态数组/静态函数     VC++2015,/Ox    2.733s

代码3. C++,动态数组/虚函数       VC++2015,/Ox    6.598s

Fortran代码比C++的模板代码慢大约50%,C++的非模板代码又比Fortran代码慢了65%,比模板代码慢了1.4倍。

测试之前没想到代码3比代码2慢一倍多,用诊断工具测了一下,下面是最耗时的几个模块:

代码2            100.00%       2.733s

RiemannSolver     20.37%       0.557s

限制器函数1       12.27%       0.335s

ucrtbase.dll       9.38%       0.256s

面通量计算         8.79%       0.240s

梯度计算           7.32%       0.200s

代码3            100.00%       6.598s

ucrtbase.dll      40.97%       2.703s

RiemannSolver      8.01%       0.528s

面通量计算         6.44%       0.425s

限制器函数2        4.26%       0.281s

梯度计算           3.50%       0.231s

代码3中仅ucrtbase.dll就耗时2.7s,相当于代码2的全部用时了。ucrtbase.dll是C运行时库,堆内存的分配函数都在里面,猜测是ArrayXd这类对象频繁的分配/释放堆内存拖慢了速度。

后续的改进应当是优化ArrayXd相关的代码,考虑自定义或使用第三方的内存池,但Eigen不支持自定义分配器比较麻烦,可能需要通过map类进行包装。最终希望代码3的效率能够达到Fotran代码1的水平,比纯模板的代码2慢50%以内。尽管付出50%的性能代价,但可以换取更高的开发维护效率,并且采用隐式时间推进后性能差距会进一步缩小,还剩下的那点效率差别就用并行解决了。

NumComp 数值计算
63 个回复
hyperLee
老李 2016-08-16

1 不用虚函数

2 关键矩阵和数组用restrict自己轮,性能至少比你现在全用eigen提升一倍。

【 在 adda (adda) 的大作中提到: 】

: 做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往C++转,底层数组使用Eigen。

: 由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,用到的数组大小/成员函数等都是静态的,多态也静态实现。求解变量5个的话,大量用到了Array3d和Array5d这类小的数组。

: 这么做刚开始还可以,但程序大了以后开发效率太低。

: ...................

hyperLee
老李 2016-08-16

还有你这种比较太粗放,没什么意义。

对了,换Intel c++,惊喜会更多。

c代码超过for太容易了

【 在 adda (adda) 的大作中提到: 】

: 做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往C++转,底层数组使用Eigen。

: 由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,用到的数组大小/成员函数等都是静态的,多态也静态实现。求解变量5个的话,大量用到了Array3d和Array5d这类小的数组。

: 这么做刚开始还可以,但程序大了以后开发效率太低。

: ...................

nobeing
nobeing 2016-08-16

学好模板元编程吧,对数值计算帮助很大。

【 在 adda (adda) 的大作中提到: 】

: 做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往

: 由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,�

: 这么做刚开始还可以,但程序大了以后开发效率太低。

: ...................

adda
adda 2016-08-16

又测了一下,用Intel C++ 16.0,版本2、3都会有6%的提升

【 在 hyperLee 的大作中提到: 】

: 还有你这种比较太粗放,没什么意义。

: 对了,换Intel c++,惊喜会更多。

: c代码超过for太容易了

: ...................

hyperLee
老李 2016-08-16

怎么搞得?你实现差异太大了吧

【 在 adda (adda) 的大作中提到: 】

: 又测了一下,用Intel C++ 16.0,版本2用时基本不变,版本3会有6%的提升

hyperLee
老李 2016-08-16

还有你动态分配太多了,cfd固定网格的话,时间步开始前一次分配够,哪有那么多动态分配的?

allocator基本用着,如果临时需要缓冲区直接在栈上分配,或者全局开个临时变量/对象就行了。

【 在 adda (adda) 的大作中提到: 】

: 做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往C++转,底层数组使用Eigen。

: 由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,用到的数组大小/成员函数等都是静态的,多态也静态实现。求解变量5个的话,大量用到了Array3d和Array5d这类小的数组。

: 这么做刚开始还可以,但程序大了以后开发效率太低。

: ...................

adda
adda 2016-08-16

搞错了,都有6%的提升

【 在 hyperLee 的大作中提到: 】

: 怎么搞得?你实现差异太大了吧

adda
adda 2016-08-16

现在只是实现了原型,后续要支持化学反应流动的计算,组分数/反应系数都从文件中读入,存储组分密度的数组只能是动态的

主要的数组肯定是时间步开始前分配好的了,但计算过程中的临时变量还是需要动态分配,比如根据格心值计算单元边界值,边界值用完就丢掉了

【 在 hyperLee 的大作中提到: 】

: 还有你动态分配太多了,cfd固定网格的话,时间步开始前一次分配够,哪有那么多动态分配的?

: allocator基本用着,如果临时需要缓冲区直接在栈上分配,或者全局开个临时变量/对象就行了。

: ...................

hyperLee
老李 2016-08-16

这也太少了吧。

【 在 adda (adda) 的大作中提到: 】

: 搞错了,都有6%的提升

hyperLee
老李 2016-08-16

这种临时变量在栈上分配就行了,或者用匿名空间存个全局变量也行,都比你局部动态分配删除好。

【 在 adda (adda) 的大作中提到: 】

: 现在只是实现了原型,后续要支持化学反应流动的计算,组分数/反应系数都从文件中读入,存储组分密度的数组只能是动态的

: 主要的数组肯定是时间步开始前分配好的了,但计算过程中的临时变量还是需要动态分配,比如根据格心值计算单元边界值,边界值用完就丢掉了

adda
adda 2016-08-16

vc和icc都开了最高优化,真就这么多提高了

相对其他编译器,icc的优势还真是没ifc大,或者说,非intel的c++编译器优化做得比fortran编译器好?

【 在 hyperLee 的大作中提到: 】

: 这也太少了吧。

adda
adda 2016-08-16

匿名空间存个全局变量这个好,还真是没想到这种做法

【 在 hyperLee 的大作中提到: 】

: 这种临时变量在栈上分配就行了,或者用匿名空间存个全局变量也行,都比你局部动态分配删除好。

laser2000
掌上智能版都是大厦挂 2016-08-16

好高端

我只敢用C,开辟类似 (*x)[5] 这种东西,再用宏实现ijk和一维id之间的换算

pasuka
青芝坞 2016-08-16

编译慢的话,拆分成几个模块,上动态或者静态链接库呗

ivf不顺手,还有gfortran、pgi等等

甚至新版MATLAB的Coder功能也可以试试,直接m脚本翻译成人类可以阅读的c

【 在 adda (adda) 的大作中提到: 】

: 做CFD的FVM/DG方法,原有一套十年历史的Fortran代码(代码1),尽管该代码模块化程度还可以,也引入了不少Fortran95/2003的技术,但程序的灵活性/重构的便利性/第三方库的支持上还是与C++差的太远,而且最新Fortran标准的编译器支持不好,用IFC中经常踩坑,因此最近在往

: 由于程序可以计算任意多个组分,对FVM,每个单元求解变量的个数5+;对DG,求解变量个数更多;单元测试中又用到了波动方程,每个单元的求解变量个数为1。为了高效的处理这些情况,首先实现了一套基于模板的代码(代码2),求解变量数、DG的阶数等作为模板的非类型参数,�

: 这么做刚开始还可以,但程序大了以后开发效率太低。

: ...................

pasuka
青芝坞 2016-08-16

大神,用宏实现会不会踩雷呢?

【 在 laser2000 (掌上智能版都是大厦挂) 的大作中提到: 】

: 好高端

: 我只敢用C,开辟类似 (*x)[5] 这种东西,再用宏实现ijk和一维id之间的换算

hyperLee
老李 2016-08-16

你这种需求,用模板实现一下,包裹个动态数组,一内联跟你的宏没有任何区别。

而且代码看着还清爽,没有雷。

你这种宏中少不了一堆括号。搞得编译器尴尬无比。

【 在 laser2000 (掌上智能版都是大厦挂) 的大作中提到: 】

: 好高端

: 我只敢用C,开辟类似 (*x)[5] 这种东西,再用宏实现ijk和一维id之间的换算

MTB
山地车 2016-08-16

Fortran 动态数组比固定数组计算慢多少?

【 在 hyperLee () 的大作中提到: 】

: 你这种需求,用模板实现一下,包裹个动态数组,一内联跟你的宏没有任何区别。

: 而且代码看着还清爽,没有雷。

: 你这种宏中少不了一堆括号。搞得编译器尴尬无比。

MTB
山地车 2016-08-16

试了一下,fortran动态数组比  静态数组慢了好几个量级,吓一跳……

【 在 hyperLee () 的大作中提到: 】

: 你这种需求,用模板实现一下,包裹个动态数组,一内联跟你的宏没有任何区别。

: 而且代码看着还清爽,没有雷。

: 你这种宏中少不了一堆括号。搞得编译器尴尬无比。

chineserice
遵白后援会de书记 2016-08-16

大神,请问学这些fortran高深知识

那本书比较合适。。。

【 在 hyperLee 的大作中提到: 】

: 这种临时变量在栈上分配就行了,或者用匿名空间存个全局变量也行,都比你局部动态分配删除好。