En Python, vous voudrez peut-être effectuer un traitement compliqué en tournant une instruction for sur un tableau multidimensionnel (les données sont souvent représentées par un tableau multidimensionnel dans le champ atmosphère et océan), mais dans ce cas, le temps d'exécution est souvent long Vous avez peut-être des ennuis. Il existe différentes méthodes d'accélération pour résoudre ce problème, mais une solution consiste à écrire un sous-programme en Fortran, qui est souvent utilisé dans le calcul fluide, et à l'appeler depuis Python. Il existe deux manières principales d'appeler une sous-routine Fortran depuis Python: en utilisant numpy.ctypeslib et en utilisant f2py. L'ordre d'accès du tableau multidimensionnel est différent entre les deux, mais je n'ai pas trouvé la page résumée en japonais, j'ai donc écrit un article.
L'environnement d'exécution est le suivant.
C'est une histoire basique, donc si vous la connaissez, vous pouvez sauter cette section.
Les tableaux multidimensionnels sont gérés dans le langage de programmation, mais ils sont disposés côte à côte dans la direction unidimensionnelle en mémoire. À ce moment-là, la façon dont il est stocké dépend de la langue. Le code ci-dessous écrit la valeur 0 dans les éléments du tableau dans l'ordre dans lequel ils ont été placés en mémoire.
loop.py
import numpy as np
nk = 2
nj = 3
ni = 4
array = np.empty((nk, nj, ni), dtype=np.float32, order="C")
for k in range(nk):
for j in range(nj):
for i in range(ni):
array[k, j, i] = 0
loop.f90
integer(4), parameter :: nk = 2, nj = 3, ni = 4
real(4), dimension(ni, nj, nk) :: array
do k = 1, nk
do j = 1, nj
do i = 1, ni
array(i, j, k) = 0
end do
end do
end do
La différence est qu'en Python (pour être exact, dans le cas de order =" C " dans le tableau numpy), l'ordre est [k, j, i], et en Fortran, l'ordre est (i, j, k). (Parfois appelé priorité de ligne et priorité de colonne, respectivement). Le fait d'y accéder dans l'ordre dans lequel ils sont stockés en mémoire ou d'y accéder de manière discrète affecte le temps d'exécution, en particulier lorsque la taille du tableau est grande.
À titre d'exemple, écrivons un code qui passe un tableau bidimensionnel de numpy à Fortran et en renvoie un double (la même chose est vraie pour trois dimensions ou plus). Les détails du code ne sont pas expliqués ici, donc si vous voulez connaître, par exemple, [Linkage with Fortran, C language](https://kyotogeopython.zawawahoge.com/html/%E5%BF%9C%E7%94%A8 % E7% B7% A8 / Fortran,% 20C% E8% A8% 80% E8% AA% 9E% 20% E3% 81% A8% E3% 81% AE% E9% 80% A3% E6% 90% BA. Voir d'autres articles tels que html).
doubling1.py
from ctypes import POINTER, c_int32, c_void_p
import numpy as np
libDoubling = np.ctypeslib.load_library("doubling1.so", ".")
libDoubling.doubling.argtypes = [
np.ctypeslib.ndpointer(dtype=np.float32),
np.ctypeslib.ndpointer(dtype=np.float32),
POINTER(c_int32),
POINTER(c_int32),
]
libDoubling.doubling.restype = c_void_p
nj = 2
ni = 3
vin = np.empty((nj, ni), dtype=np.float32, order="C")
vin[0, 0] = 1
vin[0, 1] = 2
vin[0, 2] = 3
vin[1, 0] = 4
vin[1, 1] = 5
vin[1, 2] = 6
vout = np.empty((nj, ni), dtype=np.float32, order="C")
_nj = POINTER((c_int32))(c_int32(nj))
_ni = POINTER((c_int32))(c_int32(ni))
libDoubling.doubling(vin, vout, _ni, _nj)
print(vin)
# [[1. 2. 3.]
# [4. 5. 6.]]
print(vout)
# [[ 2. 4. 6.]
# [ 8. 10. 12.]]
doubling1.f90
! gfortran -Wall -cpp -fPIC -shared doubling1.f90 -o doubling1.so
subroutine doubling(vin, vout, ni, nj) bind(C)
use, intrinsic :: ISO_C_BINDING
implicit none
integer(c_int32_t), intent(in) :: ni, nj
real(c_float), dimension(ni, nj), intent(in) :: vin
real(c_float), dimension(ni, nj), intent(out) :: vout
integer(c_int32_t) :: i, j
do j = 1, nj
do i = 1, ni
vout(i, j) = vin(i, j) * 2
end do
end do
end subroutine doubling
Un pointeur est passé à l'argument du sous-programme Fortran (bibliothèque partagée). En Python et Fortran, les index sont l'opposé de [j, i] et (i, j), mais cela devrait être simple si vous considérez l'ordre dans lequel les données sont stockées en mémoire. ..
Ensuite, utilisez f2py et écrivez le même code que précédemment. Je n'entrerai pas non plus dans les détails à ce sujet, mais je pense que vous pouvez vous référer à des articles tels que Gestion des tableaux numpy avec f2py.
doubling2.f90
! f2py -c --fcompiler=gfortran -m doubling2 doubling2.f90
subroutine doubling2(vin, vout, ni, nj)
implicit none
integer(4), intent(in) :: ni, nj
real(4), dimension(nj, ni), intent(in) :: vin
real(4), dimension(nj, ni), intent(out) :: vout
integer(4) :: i, j
do j = 1, nj
do i = 1, ni
vout(j, i) = vin(j, i) * 2
end do
end do
end subroutine doubling2
doubling2.py
import numpy as np
import doubling2
nj = 2
ni = 3
vin = np.empty((nj, ni), dtype=np.float32, order="C")
vin[0, 0] = 1
vin[0, 1] = 2
vin[0, 2] = 3
vin[1, 0] = 4
vin[1, 1] = 5
vin[1, 2] = 6
vout = np.empty((nj, ni), dtype=np.float32, order="C")
vout = doubling2.doubling2(vin, ni, nj)
print(vin)
# [[1. 2. 3.]
# [4. 5. 6.]]
print(vout)
# [[ 2. 4. 6.]
# [ 8. 10. 12.]]
Lorsque vous utilisez f2py, vous pouvez simplement utiliser «[j, i]» en Python comme «(j, i)» en Fortran. Il présente l'avantage que le code a la même apparence, mais l'instruction do de doubling2.f90 présente l'inconvénient que l'ordre d'accès du tableau ne correspond pas à l'ordre de stockage des données en mémoire.
Comment puis-je résoudre ça? En d'autres termes, je souhaite accéder au code côté Fortran dans le même ordre que doubling1.f90.
doubling3.f90
! f2py -c --fcompiler=gfortran -m doubling3 doubling3.f90
subroutine doubling3(vin, vout, ni, nj)
implicit none
integer(4), intent(in) :: ni, nj
real(4), dimension(ni, nj), intent(in) :: vin
real(4), dimension(ni, nj), intent(out) :: vout
integer(4) :: i, j
do j = 1, nj
do i = 1, ni
vout(i, j) = vin(i, j) * 2
end do
end do
end subroutine doubling3
Vous pourriez avoir l'idée que vous pouvez inverser l'ordre de stockage en mémoire en définissant l'argument order du tableau numpy sur" F " au lieu de la valeur par défaut"C". Cependant, cela ne résout pas ce problème. Autrement dit, même si vous changez order =" C " enorder = "F"dans doubling2.py, le code Fortran correspondant est toujours en train de doubler2.f90, et l'utilisation de doubling3.f90 entraînera une erreur. est.
La solution consiste à transposer le tableau en utilisant l'attribut T, comme décrit dans f2py - empêcher la réorganisation du tableau --Stack Overflow. C'est.
doubling3.py
import numpy as np
import doubling3
nj = 2
ni = 3
vin = np.empty((nj, ni), dtype=np.float32, order="C")
vin[0, 0] = 1
vin[0, 1] = 2
vin[0, 2] = 3
vin[1, 0] = 4
vin[1, 1] = 5
vin[1, 2] = 6
_vout = np.empty((nj, ni), dtype=np.float32, order="C")
_vout = doubling3.doubling3(vin.T, ni, nj)
vout = _vout.T
print(vin)
# [[1. 2. 3.]
# [4. 5. 6.]]
print(vout)
# [[ 2. 4. 6.]
# [ 8. 10. 12.]]
Par exemple, si vous souhaitez lire GPV en utilisant la bibliothèque netCDF4, le traiter efficacement avec un sous-programme Fortran, puis le dessiner avec matplotlib + cartopy, cette méthode est recommandée. Il est difficile de déplacer pour éviter «Osetsusuke», mais même s'il est inclus, la quantité de description est plus petite que lors de l'utilisation de numpy.ctypeslib, et f2py semble être plus facile à utiliser.
Recommended Posts