爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 50404|回复: 25

[混合编程] python调用fortran

[复制链接]

新浪微博达人勋

发表于 2016-12-13 10:49:59 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 王磊 于 2016-12-13 10:49 编辑

python调用fortran主要目的是节省时间同时利用一下已有的fortran程序,本人对fortran熟悉一些,对python算不上入门,仅仅是对其感兴趣。看到对f2py功能的介绍,便想尝试其对fortran的调用,故写此贴,便于以后查阅,也便于和需要的朋友分享交流。不得不说如果涉及到单个点的循环操作,最好不要用python的for循环,而是考虑用其内部函数,这会节省很多时间。尽管如此,还是赞叹python的字典功能的灵活性,对错误等的处理功能。


python调用python要点总结
(1)f2py将fortran子程序编译链接为动(静态)态库。
(2)python调用fortran需要写好接口,关键点包括:

    (a)数据类型,numpy默认的数据类型为float64也就是双精度浮点数,所以fortran里的浮点数也必须是双精度的,需要声明为real*8,默认的real是单精度的,不会报错但是会造成内存溢出!
    (b)加intent(out)之后会把python对应函数中的相应变量隐藏,也就是不需要输入了!例如fortran中的函数为 subroutingfoo(a,n)
    给n维数组a加了intent(out)之后,表示它只是输出,不接受输入数据。于是在python中调用时,只需要指定维数n就可以了,即:foo(n) 注意:接口中的虚参,不能存在intent(inout)属性,因为python不认可,如果一个虚参具备intent(inout)属性,则会报如下错:“ValueError:failed to initialize intent(inout) array -- input not fortran contiguous --input 'l' not compatible to 'd'”

(c)一般在使用fortran的时候需要注意二维数组的行列关系的问题,即是说fortran里面的数组是先填充一列、再填充一行,而python数组遵循C语法规范、先填充一行,在天聪一列,不过貌似利用f2py进行传递时,已经进行过智能转换,所以不需要管这个,按照相同的方式使用就行

(d)需要注意的是,在进行调用的时候,intent(in)的变量(输入变量)里,跟维数有关的变量必须放置在最后。比如我的函数命名为:f(a,b,c,d,n)这里面a、b、c、d都是普通的变量或数组之类的,但是n必须是一个跟维数有关的变量,用来在fortran程序里定义数组时使用。上面这条里所说的跟维数有关的变量必须放在最后,适用于n也是一个用来传递的变量的情况下。即是说你在使用时希望将维数信息用变量n来进行传递时,那么就得把它放在最后。再换句话说,n就是intent(in)的情况。 如果数组维长放在了最前面,可能出现的一个错误是:
“(shape(u,0)==m)failed for 1st keyword m: vector2:m=0”。

遵循上述几点,大致可以保证程序无误。
参考链接如下:
http://blog.sciencenet.cn/blog-350278-628083.html
http://www.jb51.net/article/82121.htm


下面是自己编写的一个测试程序,仅供参考。意思到位了,但还不够简洁。

(1)fortran程序如下,foo.f90文件内容
subroutine hello(a,b,c,d,e,f)
implicit none
integer,intent(in)::  a,b,c,d,e,f
write(*,*) a,b,c,d,e,f
end subroutine hello

subroutine vector(u,v,vec,m,n)
implicit none
integer,intent(in) ::  m, n
real(8),intent(in)    ::  u(m,n),v(m,n)
real(8),intent(out)    ::  vec(m,n)    !输出参数,python中不需要出现
integer i

vec=sqrt(u*u+v*v)

do i=1,3
    write(*,*)  "i=",i
    write(*,*)  u(i,:5)
    write(*,*)  v(i,:5)
end do
return
end subroutine vector

(2)python程序如下,python.py文件内容。
import numpy as N
def good_magnitudes(u, v,minmag=0.1):

    shape_input = N.shape(u)
    output = N.zeros(shape_input,dtype=u.dtype.char)
    for i in xrange(shape_input[0]):
        for j in xrange(shape_input[1]):
            mag = ((u[i,j]**2) + (v[i,j]**2))**0.5     #千万不要写for循环,慢的舍生忘死了
            if mag >minmag:
                output[i,j] = mag
            else:
                output[i,j] = minmag
    return output


def faster_good_magnitudes(u, v, minmag=0.1):
    mag = ((u**2) + (v**2))**0.5            #快很多哦
    output = N.where(mag >minmag, mag, minmag)
    return output



n=100000000
u=N.reshape(N.arange(n),(1000,n//1000))
v=N.reshape(N.arange(n),(1000,n//1000))
vec=N.reshape(N.arange(n),(1000,n//1000))

print u[0][:10]
print v[0][:10]

print u[1][:10]
print v[1][:10]
print u[2][:10]
print v[2][:10]




import foo                       #导入fortran子程序
import time

foo.hello(1,2,3,4,5,6)        #参数的个数不受限制


begin_time0 = time.time()
foo.vector(u,v,1000,100000)     #千万注意,维数长度的量,一定要放在最后面,当然你用变量存储维长也是可以的,也建议用变量存储
print time.time() -begin_time0

begin_time1 = time.time()  
vect2=faster_good_magnitudes(u, v, minmag=0.1)
print time.time() -begin_time1



begin_time2 = time.time()
vect1=good_magnitudes(u, v, minmag=0.1)
print time.time() -begin_time2

(3)操作步骤:


    f2py -m foo -c foo.f90
python  python.py


(4)结果如下:


[0 1 2 3 4 5 6 7 8 9]                                    对应python中:   print u[0][:10]
[0 1 2 3 4 5 6 7 8 9]                                    对应python中:   print v[0][:10]
[100000 100001 100002 100003 100004 100005 100006 100007 100008 100009]     对应python中:  print u[1][:10]
[100000 100001 100002 100003 100004 100005 100006 100007 100008 100009]     对应python中:  print v[1][:10]
[200000 200001 200002 200003 200004 200005 200006 200007 200008 200009]     对应python中:  print u[2][:10]

[200000 200001 200002 200003 200004 200005 200006 200007 200008 200009]     对应python中:  print v[2][:10]
           1           2           3           4           5           6

下面输出对应fortran中的:
do i=1,3
    write(*,*)  "i=",i
    write(*,*)  u(i,:5)
    write(*,*)  v(i,:5)
end do

i=           1
   0.0000000000000000        1.0000000000000000        2.0000000000000000        3.0000000000000000        4.0000000000000000     
   0.0000000000000000        1.0000000000000000        2.0000000000000000        3.0000000000000000        4.0000000000000000     
i=           2
   100000.00000000000        100001.00000000000        100002.00000000000        100003.00000000000        100004.00000000000     
   100000.00000000000        100001.00000000000        100002.00000000000        100003.00000000000        100004.00000000000     
i=           3
   200000.00000000000        200001.00000000000        200002.00000000000        200003.00000000000        200004.00000000000     
   200000.00000000000        200001.00000000000        200002.00000000000        200003.00000000000        200004.00000000000     

输出时间:
混合fortran片段:4.3228559494s

直接python:10.6145708561s
单点for循环: 490s+
由此输出可以发现,使用相同的行列规则既可以正确的使用数据。同时,数据量小的时候,直接python程序即可(使用一些技巧提升速度),不需要混合调用,除非非得使用单次循环。

评分

参与人数 1金钱 +15 贡献 +5 收起 理由
mofangbao + 15 + 5

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-5-6 23:32:50 | 显示全部楼层
有个问题想请教一下!我们大气海洋中的计算模式f90文件常常都是嵌套的,也就是一个module里面要use 另一个module或者subroutine,这种情况下,是不是必须解除这种依赖性?也就是说,用python调用fortran计算,由python出图,这种想法只能适用于单独的一个subroutine计算,该subroutine不能牵连到其他f90文件中的subroutine。有没有解决的方法呢?最近在弄这个,这一点有点困惑,我想的解决方法是把所有关联的f90文件都整合在一个f90文件中,让其变成一个单独的subroutine。但是鉴于有一些计算模式f90文件相当多,适用范围太窄,有没有更好的解决方法呢?
密码修改失败请联系微信:mofangbao
回复 支持 2 反对 0

使用道具 举报

新浪微博达人勋

发表于 2017-2-23 16:54:18 | 显示全部楼层
还没仔细看,我现在也在想是用python调用Fortran呢?还是自己写python的代码?收藏了!谢谢!
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2016-12-13 15:02:29 | 显示全部楼层
细节总结的很赞!不知这位王磊是不是本科的同学    总之受益 收藏
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-12-13 16:55:10 | 显示全部楼层
小伙子,干得不错!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-19 11:45:49 | 显示全部楼层
学习了,谢谢分享!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-19 12:26:47 | 显示全部楼层
我明明就定义了 inout..
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-9 19:40:11 | 显示全部楼层
收藏了,感谢分享
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-10 10:55:00 | 显示全部楼层
有些细节可能和编译器及系统平台相关。
关于
  1. 跟维数有关的变量必须放置在最后
复制代码
这条,在windows上编译时,位置是无关的。但在linux上运行时,维数就必须放到相关变量的前面,否则就会出错。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-18 11:26:52 | 显示全部楼层
谢谢!我自己还要手动实践一下!再细细琢磨下!
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表