I have finally managed to wrap some Fortran into Python using numpy's f2py. I am going to need to write a small library of subroutines that I can then use for a larger project, but I am already having trouble with the example subroutine I just compiled.
This is a simple algorithm to solve linear systems of equations in Fortran 90:
SUBROUTINE solve_lin_sys(A, c, x, n)
! =====================================================
! Uses gauss elimination and backwards substitution
! to reduce a linear system and solve it.
! Problem (0-div) can arise if there are 0s on main diag.
! =====================================================
IMPLICIT NONE
INTEGER:: i, j, k
REAL*8::fakt, summ
INTEGER, INTENT(in):: n
REAL*8, INTENT(inout):: A(n,n)
REAL*8, INTENT(inout):: c(n)
REAL*8, INTENT(out):: x(n)
DO i = 1, n-1 ! pick the var to eliminate
DO j = i 1, n ! pick the row where to eliminate
fakt = A(j,i) / A(i,i) ! elimination factor
DO k = 1, n ! eliminate
A(j,k) = A(j,k) - A(i,k)*fakt
END DO
c(j)=c(j)-c(i)*fakt ! iterate on known terms
END DO
END DO
! Actual solving:
x(n) = c(n) / A(n,n) ! last variable being solved
DO i = n-1, 1, -1
summ = 0.d0
DO j = i 1, n
summ = summ A(i,j)*x(j)
END DO
x(i) = (c(i) - summ) / A(i,i)
END DO
END SUBROUTINE solve_lin_sys
I successfully made it into a Python module and I can import it, but I can't figure out what kind of arguments I should pass it into python.
This is what I've tried:
import numpy as np
import fortran_operations as ops
sys = [[3, -0.1, -0.2],
[0.1, 7, -0.3],
[0.3, -0.2, 10]]
known = [7.85, -19.3, 71.4]
A = np.array(sys)
c = np.array(known)
x = ops.solve_lin_sys(A, c, 3)
I tried using both the normal list and the numpy array as an argument for the Fortran function, but in all cases I get the error:
x = ops.solve_lin_sys(A, c, 3) ------ ValueError: failed in converting 1st argument `a' of fortran_operations.solve_lin_sys to C/Fortran array
I thought you were supposed to simply use numpy arrays. Printing solve_lin_sys.doc also states that I should use rank 2 arrays for A and rank 1 for C. Isn't that already the shape of the arrays when I defined them like that?
Sorry if this has a stupid solution, I'm a beginner with all of this.
Thank you in advance
CodePudding user response:
Arrays of dimension 2 and higher must be continuous with 'Fortran' order of elements.
import ops
sys = [[3, -0.1, -0.2],
[0.1, 7, -0.3],
[0.3, -0.2, 10]]
known = [7.85, -19.3, 71.4]
A = np.array(sys, dtype='d', order='F')
c = np.array(known, dtype='d')
x = ops.solve_lin_sys(A, c, 3)