returning values from FORTRAN with f2py

[basegmez]

Chris Fonnesbeck wrote in message ...
>I am wondering how to persuade some fortran code that I have compiled
>with f2py to return its calculated values to me in python.  I have
>some code that does gaussian elimination (mostly from Numerical
>Reciopes), and is supposed to return an inverted matrix so that I can
>solve a system of linear equations.  However, when run in python, the
>function runs without error but returns no values:
>
>>>> a = [[2,1,3],[1,-1,1],[3,6,-1]]
>>>> b = [7,2,10]
>>>> import gauss
>>>> gauss.gauss(a,3,b,1)
>>>>
>
>Is there something that I have to change in my fortran code to make
>this work?
>
>Many thanks,
>cjf

I just used Numeric to do something similar.  I used Newton-Raphson to solve
a set of non-linear equations.

from __future__ import nested_scopes
from Numeric import array
from Numeric import matrixmultiply as mm
from Matrix import inverse
g = 386.08

''' Given Ix, Iy, Iz, and weight, calculates X, Y and Z dims of a
rectangular box
    uses Newton-Raphson to solve three simultaneous nonlinear equations'''
def dim_box(w, Ix, Iy, Iz, tol = .0001):
    m = w/g
    mc = m/12
    def f(x, y, z):   # Function
        return array([(y**2 + z**2)*mc - Ix, (x**2 + z**2)*mc - Iy, (x**2 +
y**2)*mc - Iz])
    def J(x, y, z):    # Jacobian
        return array([[0, 2*y*mc, 2*z*mc], [2*x*mc, 0, 2*z*mc], [2*x*mc,
2*y*mc, 0]])
    x, y, z = 1.8, 4, 5
    for i in range(200):
        res = mm(inverse(J(x, y, z)), -f(x, y, z)) + array([x, y, z])
        x, y, z = res[0], res[1], res[2]
        check = f(x, y, z)
        if abs(check[0]**2 + check[1]**2 + check[2]**2) <= tol:
            return x, y, z
    print 'No solution!', '\n', f(x, y, z), '\n', x, y, z

# Example:
# >>> dim_box(1000, 1600, 700, 1800)
# (45.6599607331 79.0853715939 34.0329252342)


Alternatively, you can use the solve_linear_equations from
LinearAlgebraModule.

from Numeric import array
from LinearAlgebra import solve_linear_equations as solve

#Example:
# a + b = 5
# 2a - b = 7
#should give
# [4., 1.]

LHS = array([[1, 1], [2, -1]])
RHS = array([5, 7])
print solve(LHS, RHS)


Fahri Basegmez