Почему мои программы f2py медленнее, чем программы python

Question

Недавно я написал трудоемкую программу с python и решил переписать наиболее трудоемкую часть с fortran.

Однако код fortran, завернутый в f2py, медленнее, чем код python. Может ли кто-нибудь сказать мне, как найти, что здесь происходит?

Для справки, вот функция Python:

def iterative_method(alpha0, beta0, epsilon0, epsilons0, omega, smearing=0.01, precision=0.01, max_step=20, flag=0): 
    # alpha0, beta0, epsilon0, epsilons0 are numpy arrays 
    m, n = np.shape(epsilon0) 
    Omega = np.eye(m, dtype=np.complex) * (omega + smearing * 1j) 
    green = LA.inv(Omega - epsilon0) # LA is numpy.linalg 
    alpha = np.dot(alpha0, np.dot(green, alpha0)) 
    beta = np.dot(beta0, np.dot(green, beta0)) 
    epsilon = epsilon0 + np.dot(alpha0, np.dot(green, beta0)) + np.dot(beta0, np.dot(green, alpha0)) 
    epsilons = epsilons0 + np.dot(alpha0, np.dot(green, beta0)) 

    while np.max(np.abs(alpha0)) > precision and np.max(np.abs(beta0)) > precision and flag < max_step: 
     flag += 1 
     return iterative_method(alpha, beta, epsilon, epsilons, omega, smearing, precision, min_step, max_step, flag) 
return epsilon, epsilons, flag 

Соответствующий Fortran код

SUBROUTINE iterate(eout, esout, alpha, beta, e, es, omega, smearing, prec, max_step, rank) 
    INTEGER, PARAMETER :: dp = kind(1.0d0) 
    REAL(kind=dp) :: omega, smearing, prec 
    INTEGER :: max_step, step, rank, cnt 
    COMPLEX(kind=dp) :: alpha(rank,rank), beta(rank,rank), omega_mat(rank, rank),& 
    green(rank, rank), e(rank,rank), es(rank,rank) 
    COMPLEX(kind=dp), INTENT(out) :: eout(rank, rank), esout(rank, rank) 
    step = 0 
    omega_mat = 0 
    DO cnt=1, rank 
     omega_mat(cnt, cnt) = 1.0_dp 
    ENDDO 
    omega_mat = omega_mat * (omega + (0.0_dp, 1.0_dp) * smearing) 
    DO WHILE (maxval(abs(alpha)) .gt. prec .or. maxval(abs(beta)) .gt. prec .and. step .lt. max_step) 
     green = zInverse(rank, omega_mat - e) ! zInverse is calling lapack to compute inverse of the matrix 
     e = e + matmul(alpha, matmul(green, beta)) + matmul(beta, matmul(green, alpha)) 
     es = es + matmul(alpha, matmul(green, beta)) 
     alpha = matmul(alpha, matmul(green, alpha)) 
     beta = matmul(beta, matmul(green, beta)) 
     step = step + 1 
    ENDDO 
    print *, step 
    eout = e 
    esout = es 
END SUBROUTINE iterate 

В тесте, Python код, используемый около 5 секунд, пока Fortran код, используемый около 7 секунд, что вряд ли приемлемо. Кроме того, я вряд ли вижу накладные расходы в коде fortran. Виноват ли обертка?

Редактировать: Я не использовал BlAS для matmul. После использования BLAS, fortran и python performace составляют около 5 секунд.

Answer 1

Прежде всего, сделайте this на код python, чтобы вы точно знали, как он проводит свое время. Затем вы можете сделать аналогичную вещь в коде Fortran с помощью отладчика, если хотите.

Я подозреваю, что в основном все времени переходит в матричные операции, поэтому любая разница в скорости происходит из-за математической библиотеки, а не от языка, который ее вызывает. This post передает некоторые из моих опытов. Часто подпрограммы, выполняющие такие вещи, как преобразование матриц, обратное или преобразование Холецкого, предназначены для работы на больших матрицах, но не на небольших.

Например, подпрограмма LAPACK-матричного умножения DGEMM имеет два символьных аргумента: TRANSA и TRANSB, которые могут быть в верхнем или нижнем регистре, указав, переносится ли каждая входная матрица. Чтобы проверить значение этих аргументов, он вызывает функцию LSAME. Я обнаружил, что если я трачу большую часть своего времени на умножение малых матриц, например, 4x4, программа фактически тратит почти все время на вызов LSAME и очень мало времени на самом деле умножает матрицы. Вы можете видеть, как легко это исправить.