요즘엔 포트란보다 파이썬(눔피)으로 주로 계산하는데 주목적은 메스메티카로 하면 쉽고 금방 하지만 시간낭비하고 싶을때 많이 쓴다.
파이썬의 장점은 포트란보다 짜기 편하다는데 있지만 막상 나처럼 쪼렙들이 프로그램을 짜면 포트란에 비해 엄청 느리다. 그러니까 파이썬, 포트란 둘다 못하는 사람은 그냥 닥치고 포트란을 배우는게 났다...ㅜㅜ 가 아니라 그래도 조금이라도 공부한 의리가 있으니 나중을 위해 기록...
파이썬에선 베셀함수를 쓰려면 그냥 유니버설 평션인 sph_jn (
링크) 을 쓰면 되지만 이걸쓰다보면 어이없는게 함수처럼 쓸수있는것도 아니면서 n 번째 오더까지 함수값과 미분값을 같이 준다. 따라서 쓸데없이 계산시간만 길어지는데 그래서 생각해낸 방법이 f2py 란걸 이용해서 그냥 새로 만들어 쓰는법. f2py 는 포트란 코드를 파이썬 모듈처럼 쓸수 있게 해주는 것이다. (
링크) 나도 아직 제대로 쓸줄은 모르는데 심심할때마다 공부중... 아무튼 어찌어찌 성공해서 프로그램을 돌려보니 시간 차이가 많이 난다.
| 1 | import numpy as np |
| 2 | import time |
| 3 | from SBESS import sbess |
| 4 | from scipy.special import sph_jn |
| 5 | |
| 6 | n = 10 |
| 7 | x = 2 |
| 8 | |
| 9 | print sph_jn(n,x)[0][n] |
| 10 | print np.nan_to_num(sbess(x,n)) |
| 11 | |
| 12 | startTime = time.time() |
| 13 | |
| 14 | for i in range(10000): |
| 15 | sph_jn(n,x)[0][n] |
| 16 | |
| 17 | elapsedTime = time.time() - startTime |
| 18 | |
| 19 | |
| 20 | print ' Calculate Run Time %5.8f' % elapsedTime |
| 21 | |
| 22 | |
| 23 | startTime = time.time() |
| 24 | |
| 25 | for i in range(10000): |
| 26 | sbess(x,n) |
| 27 | |
| 28 | elapsedTime = time.time() - startTime |
| 29 | |
| 30 | |
| 31 | print ' Calculate Run Time %5.8f' % elapsedTime |
| 32 |
 |
스페리컬 베설 제이 n=10을 x=2 일때 값을 불러오는 코드 이걸 10000번 반복한다. sph_jn 은 원래 싸이파이(눔피) 내장된 함수이고 sbess는 f2py로 만든것. nan_to_num 은 소수 16번째 이하는 0을 안주고 NaN 이라고 값을 주길래 써줌. 이런 계산을 만번 반복할 일은 거의 없겠지만 막상 해봤을때 속도 차이가 20배정도 난다. 1초 걸릴 프로그램이 20초 걸린다면 기다리다 못해 늙어 죽을수도 있으므로 나름 성공이라고 할수 있다.
f2py 를 쓰기 위해선 당연히 잘 돌아가는 포트란 코드가 필요하다. 그리고 명령어 몇번 치면 된다. f2py는 우분투에선 그냥 시냅틱 패키지 관리자에서 받을수 있다. 아래는 그냥 내가 쓰고있는 스페리컬 베셀 포트란 코드. 나도 어찌어찌 받은건데 이런거 막 올려도 되는지 잘 모르겠다. (안되면 자삭) sbess.f90 이 들어있는 폴더로 이동한 뒤 터미널에서 아래 첫번쨰 명령어를 실행하면 sbess.pyf 가 생성되고 두번째 명령어를 실행시키면 sbess.so 가 생성된다. 그러면 위 프로그램처럼 그냥 import 로 불러 올수 있다.
! f2py -m SBESS -h SBESS.pyf SBESS.f90
! f2py -c SBESS.pyf SBESS.f
sbess.f90
|
| 1 | |
| 2 | MODULE nrtype |
| 3 | INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9) |
| 4 | INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4) |
| 5 | INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) |
| 6 | INTEGER, PARAMETER :: SP = KIND(1.0) |
| 7 | INTEGER, PARAMETER :: DP = KIND(1.0D0) |
| 8 | INTEGER, PARAMETER :: SPC = KIND((1.0,1.0)) |
| 9 | INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0)) |
| 10 | INTEGER, PARAMETER :: LGT = KIND(.true.) |
| 11 | REAL(SP), PARAMETER :: PI=3.141592653589793238462643383279502884197_sp |
| 12 | REAL(SP), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_sp |
| 13 | REAL(SP), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_sp |
| 14 | REAL(SP), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_sp |
| 15 | REAL(SP), PARAMETER :: EULER=0.5772156649015328606065120900824024310422_sp |
| 16 | REAL(DP), PARAMETER :: PI_D=3.141592653589793238462643383279502884197_dp |
| 17 | REAL(DP), PARAMETER :: PIO2_D=1.57079632679489661923132169163975144209858_dp |
| 18 | REAL(DP), PARAMETER :: TWOPI_D=6.283185307179586476925286766559005768394_dp |
| 19 | TYPE sprs2_sp |
| 20 | INTEGER(I4B) :: n,len |
| 21 | REAL(SP), DIMENSION(:), POINTER :: val |
| 22 | INTEGER(I4B), DIMENSION(:), POINTER :: irow |
| 23 | INTEGER(I4B), DIMENSION(:), POINTER :: jcol |
| 24 | END TYPE sprs2_sp |
| 25 | TYPE sprs2_dp |
| 26 | INTEGER(I4B) :: n,len |
| 27 | REAL(DP), DIMENSION(:), POINTER :: val |
| 28 | INTEGER(I4B), DIMENSION(:), POINTER :: irow |
| 29 | INTEGER(I4B), DIMENSION(:), POINTER :: jcol |
| 30 | END TYPE sprs2_dp |
| 31 | END MODULE nrtype |
| 32 | |
| 33 | !=====================================================================! |
| 34 | FUNCTION Sbess(X,N) |
| 35 | !=====================================================================! |
| 36 | use nrtype |
| 37 | |
| 38 | implicit none |
| 39 | real(dp), intent(in) :: x |
| 40 | integer(i4b), intent(in) :: n |
| 41 | real(dp) :: sbess |
| 42 | |
| 43 | LOGICAL(lgt) Down,Noacc |
| 44 | |
| 45 | real(dp), PARAMETER :: Eps=1.e-9_dp |
| 46 | real(dp), PARAMETER :: P1=0.1_dp, P35=0.35_dp, & |
| 47 | P5=0.5_dp, C48=48._dp |
| 48 | ! |
| 49 | real(dp), DIMENSION(0:999) :: a |
| 50 | ! |
| 51 | integer(i4b) :: i, j, nu, nu1, nu2, np, lam, mu |
| 52 | real(dp) :: fj, fj1, fj2, ap, b, Alog2e, u, alpha, xabs |
| 53 | !----------------------------------------------------------------------- |
| 54 | ! |
| 55 | Noacc=.FALSE. |
| 56 | xabs=ABS(X) |
| 57 | IF(xabsTHEN |
| 58 | IF(N==0) THEN |
| 59 | Fj=1._dp |
| 60 | ELSE |
| 61 | Fj=0._dp |
| 62 | END IF |
| 63 | ELSE |
| 64 | IF(N>250) THEN |
| 65 | print*,"*** ORDER OF BESSELFUNCTION TOO LARGE ***" |
| 66 | STOP |
| 67 | END IF |
| 68 | IF(xabs*xabs>=real(N*(N+1))) THEN |
| 69 | Fj=SIN(xabs)/xabs |
| 70 | IF(N>0) THEN |
| 71 | Fj1=Fj |
| 72 | Fj=Fj1/xabs-COS(xabs)/xabs |
| 73 | IF(N>1) THEN |
| 74 | Fj2=Fj |
| 75 | DO J=2,N |
| 76 | Fj=real(2*J-1)*Fj2/xabs-Fj1 |
| 77 | Fj1=Fj2 |
| 78 | Fj2=Fj |
| 79 | enddo |
| 80 | END IF |
| 81 | END IF |
| 82 | ELSE |
| 83 | Ap=P1 |
| 84 | B=P35 |
| 85 | Alog2e=C48 |
| 86 | U = 2._dp * xabs/real(2*N+1) |
| 87 | Nu1 = N + INT(Alog2e*(Ap+B*U*(2._dp-U*U)/(2._dp*(1._dp-U*U)))) |
| 88 | Np=INT(xabs-P5+SQRT(B*xabs*Alog2e)) |
| 89 | IF(N>Np) THEN |
| 90 | U=2._dp*xabs/real(2*NP+1) |
| 91 | Nu2=Np+INT(Alog2e*(Ap+B*U*(2._dp-U*U)/(2._dp*(1._dp-U*U)))) |
| 92 | ELSE |
| 93 | Nu2=100000 |
| 94 | END IF |
| 95 | Nu=MIN(Nu1,Nu2,997) |
| 96 | ! IF(Nu>=997) THEN |
| 97 | ! WRITE(6,*) '*** AccURACY NOT HIGH ENOUGH ***' |
| 98 | ! Noacc=.TRUE. |
| 99 | ! END IF |
| 100 | |
| 101 | Down=.TRUE. |
| 102 | A(Nu+1)=0._dp |
| 103 | DO I=Nu,1,-1 |
| 104 | IF(Down) THEN |
| 105 | A(I)=xabs/(real(2*I+1)-xabs*A(I+1)) |
| 106 | IF(A(I)>1._dp) THEN |
| 107 | Down=.FALSE. |
| 108 | Lam=I |
| 109 | A(I-1)=1._dp |
| 110 | END IF |
| 111 | ELSE |
| 112 | A(I-1)=real(2*I+1)*A(I)/xabs-A(I+1) |
| 113 | END IF |
| 114 | enddo |
| 115 | IF(Down) THEN |
| 116 | A(0)=SIN(xabs)/xabs |
| 117 | Lam=0 |
| 118 | ELSE |
| 119 | Alpha=1._dp/((A(0)-xabs*A(1))*COS(xabs)+xabs*A(0)*SIN(xabs)) |
| 120 | Mu=MIN(Lam,N)+1 |
| 121 | DO J=0,Mu-1 |
| 122 | A(J)=Alpha*A(J) |
| 123 | enddo |
| 124 | END IF |
| 125 | DO J=Lam+1,N |
| 126 | A(J)=A(J)*A(J-1) |
| 127 | enddo |
| 128 | Fj=A(N) |
| 129 | END IF |
| 130 | END IF |
| 131 | ! |
| 132 | Sbess=Fj |
| 133 | |
| 134 | RETURN |
| 135 | END FUNCTION Sbess |
| 136 | |
| 137 | |
| 138 | |
끝. |
댓글
댓글 쓰기