f2py, проблемы прохождение функции Python для Fortran

Question

У меня есть этот простой Fortran код (stack.f90):

 subroutine fortran_sum(f,xs,nf,nxs) 
     integer nf,nxs 
     double precision xs,result 
     dimension xs(nxs),result(nf) 
     external f 
     result = 0.0 
     do I = 1,nxs 
     result = result + f(xs(I)) 
     print *,xs(I),f(xs(I)) 
     enddo 
     return 
     end 

Что я компиляция с помощью:

f2py -c --compiler=mingw32 -m stack2 stack2.f90 

А затем проверить с помощью этого Python скрипта (stack.py):

import numpy as np 
from numpy import cos, sin , exp 
from stack import fortran_sum 
def func(x): 
    return x**2 

if __name__ == '__main__': 
    xs = np.linspace(0.,10.,10) 
    ans = fortran_sum(func,xs,1) 
    print 'Fortran:',ans 
    print 'Python:',func(xs).sum() 

Когда я бегу, используя "python stack.py" это дает:

0.0000000000000000  0.00000000 
    1.1111111111111112    Infinity 
    2.2222222222222223    Infinity 
    3.3333333333333335  9.19089638E-26 
    4.4444444444444446    Infinity 
    5.5555555555555554  0.00000000 
    6.6666666666666670  9.19089638E-26 
    7.7777777777777786  1.60398502E+09 
    8.8888888888888893    Infinity 
    10.000000000000000  0.00000000 
Fortran: None 
Python: 351.851851852 

Мои вопросы:

• почему функция не оценивается правильно?

• как вернуть result в Python?

• можно оценить массив xs сразу в Фортране?

Спасибо!

---

EDIT: С больших советами от @SethMMorton я пришел заканчивал со следующим:

 subroutine fortran_sum(f,xs,nxs,result) 
     implicit none 
     integer :: I 
     integer, intent(in) :: nxs 
     double precision, intent(in) :: xs(nxs) 
     double precision, intent(out) :: result 
     double precision :: f 
     external :: f 
     ! "FIX" will be added here 
     result = 0.0 
     do I = 1,nxs 
     result = result + f(xs(I)) 
     print *,xs(I),f(xs(I)) 
     enddo 
     return 
     end 

Запуска stack.py с этой командой модифицированной: ans = fortran_sum(func,xs); дает:

0.0000000000000000  0.0000000000000000 
    1.1111111111111112  3.8883934247189009E+060 
    2.2222222222222223  3.8883934247189009E+060 
    3.3333333333333335  9.1908962428537221E-026 
    4.4444444444444446  3.8883934247189009E+060 
    5.5555555555555554  5.1935286092977251E-060 
    6.6666666666666670  9.1908962428537221E-026 
    7.7777777777777786  1603984978.1728516 
    8.8888888888888893  3.8883934247189009E+060 
    10.000000000000000  0.0000000000000000 
Fortran: 1.55535736989e+61 
Python: 351.851851852 

Это не так. Это странное поведение не происходит, если я добавляю промежуточную переменную x=x(I) И вызывается функция, используя эту переменную f(x). Смешная вещь заключается в том, что если я позвоню f(x) один раз, то и нужный звонок f(x(I)) также работает. После применения этого «FIX»:

double precision :: x, f_dummy 
x = xs(I) 
f_dummy = f(x) 

затем компиляции и запуска, получается правильный результат:

0.0000000000000000  0.0000000000000000 
    1.1111111111111112  1.2345679
    2.2222222222222223  4.9382716049382722 
    3.3333333333333335  11.111111111111112 
    4.4444444444444446  19.753086419753089 
    5.5555555555555554  30.864197530864196 
    6.6666666666666670  44.444444444444450 
    7.7777777777777786  60.493827160493836 
    8.8888888888888893  79.
    10.000000000000000  100.00000000000000 
Fortran: 351.851851852 
Python: 351.851851852 

Было бы неплохо, если бы кто-то может объяснить, почему?

Answer 1

Почему функция не оценивается правильно?

Вы отправляете результат функции обратного вызова f непосредственно на метод print. Поскольку Fortran не имеет представления о том, какой тип данных f будет возвращен (поскольку он не объявлен и не является функцией Fortran), print не может правильно интерпретировать данные для правильной печати.если присвоить результат переменной, а затем распечатать переменную вы должны увидеть ожидаемый результат:

double precision x 
.... 
x = f(xs(1)) 
print *, x 

Как вернуть result в Python?

Вам необходимо сообщить f2py о намерениях переменных. На самом деле это можно сделать со стандартным синтаксисом Fortran 90, который я настоятельно рекомендую делать, поскольку он делает ваш код более четким (я знаю, что вы используете Fortran 90, потому что вы можете печатать весь массив, не зацикливая на нем).

subroutine stack2(f,xs,result,nf,nxs) 
    implicit none 
    integer,   intent(in) :: nf, nxs 
    double precision, intent(in) :: xs(nxs) 
    double precision, intent(out) :: result(nf) 
    external f 
    ... 
end subroutine stack2 

Теперь ясно, какие переменные входят в рутину и которые выходят. Обратите внимание, что result должен быть аргументом для подпрограммы, чтобы ее можно было отключить.

f2py теперь поймет, что result должен быть возвращен с функции stack2.

Можно ли оценить массив xs сразу в Fortran

Я не уверен, что вы имеете в виду именно, но я предполагаю, что вы хотите знать, если вы можете сделать это на весь массив в точке а времени, а не в цикле do.

Возможно

Поскольку xs массив, вы должны быть в состоянии просто передать весь массив f и должны работать на весь массив. Это может не сработать, потому что Fortran требует, чтобы функция была ELEMENTAL, и это может быть не в пределах f2py. Если f2py не выполняет функцию ELEMENTAL, вы должны сделать это в цикле.

Вопрос вы можете обратиться

Линия result = result + f(xs(I)) является назначение такое же значение для каждого элемента result. Неясно, действительно ли это то, что вы хотите, но это может быть проблемой, если это не так.

код на основе резюме

Ниже приведен пример версии кода, используя эти новые предложения.

subroutine stack2(f,xs,result,nf,nxs) 
    implicit none 
    integer,   intent(in) :: nf, nxs 
    double precision, intent(in) :: xs(nxs) 
    double precision, intent(out) :: result(nf) 
    double precision :: x 
    integer   :: I 
    external f 
    result = 0.0d0 ! Make this a double constant since result is double 
    do I = 1,nxs 
     x = f(xs(I)) 
     result = result + x 
     print *, x 
    end do 
    print *, xs 
    return ! Not needed in Fortran 90 
end subroutine stack2 

Обратите внимание, что f2py будет правильно выяснить длину вашего входного массива, поэтому при вызове его нужно только определить длину result:

import numpy as np 
from stack2 import stack2 
def func(x): 
    return x**2 

if __name__ == '__main__': 
    xs = np.linspace(0.,10.,10) 
    ans = stack2(func,xs,5) # This was changed 
    print 'ans:',ans