NoviceCoder
NoviceCoder

Reputation: 474

Can I use libraries like numpy and scipy in fortran?

Im trying to speed up my python code by porting a bunch of my nested loops over to fortran and calling them as subroutines.

But alot of my loops call numpy, and special functions from scipy like bessel functions.

Before I try and use fortran I was wondering if it was possible to import scipy and numpy to my fortran subroutine and call the modules for bessel functions?

Else would I have to create the bessel function in fortran in order to use it?

Ideally, I would create some sort of subroutine that would optimize this code below. This is just a snippet of my entire project to give you an idea of what I'm trying to accomplish.

I understand that there are other practices I should implement to improve the speed, but for now I was investigating the benefits of calling fortran subroutines in my main python program.

    for m in range(self.MaxNum_Eigen):
        #looping throught the eigenvalues for the given maximum number of eigenvalues allotted
        bm = self.beta[m]


        #not sure
        #*note: rprime = r. BUT tprime ~= t.

        #K is a list of 31 elements for this particular case

        K = (bm / math.sqrt( (self.H2**2) + (bm**2) ))*(math.sqrt(2) / self.b)*((scipy.special.jv(0, bm * self.r))/ (scipy.special.jv(0, bm * self.b))) # Kernel, K0(bm, r).

        #initial condition
        F = [37] * (self.n1)

        # Integral transform of the initial condition
        #Fbar = (np.trapz(self.r,self.r*K*F))

        '''
            matlab syntax trapz(X,Y), x ethier spacing or vector
            matlab:     trapz(r,r.*K.*F)                trapz(X,Y)
            python:     np.trapz(self.r*K*F, self.r)    trapz(Y,X)

        '''

        #*(np.trapz(self.r,self.r*K*F))
        Fbar = np.ones((self.n1,self.n2))*(np.trapz(self.r*K*F, self.r))

        #steady state condition: integral is in steady state
        SS = np.zeros((sz[0],sz[1]))

        coeff = 5000000*math.exp(-(10**3)) #defining value outside of loop with higher precision

        for i in range(sz[0]):
            for j in range(sz[1]):

                '''
                    matlab reshape(Array, size1, size2) takes multiple arguments the item its resizeing and the new desired shape

                    create self variables and so we are not re-initializing them over and over agaian?

                    using generators? How to use generators

                '''
                s = np.reshape(tau[i,j,:],(1,n3))


                # will be used for rprime and tprime in Ozisik solution.
                [RR,TT] = np.meshgrid(self.r,s)



                '''
                    ##### ERROR DUE TO ROUNDING OF HEAT SOURCE ####

                    error in rounding  5000000*math.exp(-(10**3)) becomes zero

                    #log10(e−10000)=−10000∗(0.4342944819)=−4342.944819

                    #e−1000=10−4342.944819=10−4343100.05518=1.13548386531×10−4343

                '''

                #g = 5000000*math.exp(-(10**3)) #*(RR - self.c*TT)**2) #[W / m^2] heat source.
                g = coeff * (RR - self.c*TT)**2


                K = (bm/math.sqrt(self.H2**2 + bm**2))*(math.sqrt(2)/self.b)*((scipy.special.jv(0,bm*RR))/(scipy.special.jv(0,bm*self.b)))

                #integral transform of heat source
                gbar = np.trapz(RR*K*g, self.r, 2) #trapz(Y,X,dx (spacing) )
                gbar = gbar.transpose()



                #boundary condition. BE SURE TO WRITE IN TERMS OF s!!!
                f2 = self.h2 * 37

                A = (self.alpha/self.k)*gbar + ((self.alpha*self.b)/self.k2)*((bm/math.sqrt(self.H2**2 + bm**2))*(math.sqrt(2)/self.b)*((scipy.special.jv(0,bm*self.b))/(scipy.special.jv(0,bm*self.b))))*f2
                #A is essentially a constant is this correct all the time?
                #What does A represent

                SS[i, j] = np.trapz(np.exp( (-self.alpha*bm**2)*(T[i,j] - s) )*A, s)

        #INSIDE M loop
        K = (bm / math.sqrt((self.H2 ** 2) + (bm ** 2)))*(math.sqrt(2) /self.b)*((scipy.special.jv(0, bm * R))/ (scipy.special.jv(0, bm * self.b)))


        U[:,:, m] = np.exp(-self.alpha * bm ** 2 * T)* K* Fbar + K* SS

        #print(['Eigenvalue ' num2str(m) ', found at time ' num2str(toc) ' seconds'])

Upvotes: 2

Views: 881

Answers (1)

NoviceCoder
NoviceCoder

Reputation: 474

Compilation of answers given in the comments

Answers specific to my code: As vorticity mentioned my code in itself was not using the numpy, and scipy packages to the fullest extent.

In regards to Bessel, function 'royvib' mentions using using .jo from scipy rather than .jv. Calling the special Bessel function jv. is much more computationally expensive, especially since I knew that I would be using a zeroth order bessel function for many of my declarations the minor change from jv -> j0 solved speed up the process.

In addition, I declared variables outside the loop to prevent expensive calls to searching for my appropriate functions. Example below.

Before

for i in range(SomeLength):
    some_var = scipy.special.jv(1,coeff)

After

Bessel = scipy.special.jv
for i in range(SomeLength):
    some_var = Bessel(1,coeff)

Storing the function saved time by not using the dot ('.') the command to look through the libraries every single loop. However keep in mind this does make python less readable, which is the main reason I choose to do this project in python. I do not have an exact amount of time this step cut from my process.

Fortran specific: Since I was able to improve my python code I did not go this route an lack of specifics, but the general answer as stated by 'High Performance Mark' is that yes there are libraries that have been made to handle Bessel functions in Fortran.

If I do port my code over to Fortran or use f2py to mix Fortran and python I will update this answer accordingly.

Upvotes: 1

Related Questions