In Python möchten Sie möglicherweise eine komplizierte Verarbeitung durchführen, indem Sie eine for-Anweisung für ein mehrdimensionales Array aktivieren (Daten werden häufig durch ein mehrdimensionales Array in der Atmosphäre und im Ozeanfeld dargestellt). In diesem Fall ist die Ausführungszeit jedoch häufig lang Sie könnten in Schwierigkeiten sein. Es gibt verschiedene Beschleunigungsmethoden, um dies zu lösen. Eine Lösung besteht jedoch darin, eine Unterroutine in Fortran zu schreiben, die häufig in der Flüssigkeitsberechnung verwendet wird, und sie von Python aus aufzurufen. Es gibt zwei Möglichkeiten, eine Fortran-Subroutine aus Python aufzurufen: Verwenden von numpy.ctypeslib und Verwenden von f2py. Die Zugriffsreihenfolge des mehrdimensionalen Arrays unterscheidet sich zwischen den beiden, aber ich konnte die auf Japanisch zusammengefasste Seite nicht finden, also habe ich einen Artikel geschrieben.
Die Ausführungsumgebung ist wie folgt.
Es ist eine grundlegende Geschichte. Wenn Sie sie kennen, können Sie diesen Abschnitt überspringen.
Obwohl die Programmiersprache mehrdimensionale Arrays verarbeitet, sind sie im Speicher nebeneinander in eindimensionaler Richtung angeordnet. Zu diesem Zeitpunkt hängt es von der Sprache ab, wie es gespeichert wird. Im folgenden Code wird der Wert 0 in der Reihenfolge in die Elemente des Arrays geschrieben, in der sie gespeichert werden.
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
Der Unterschied besteht darin, dass in Python (um genau zu sein, im Fall von "order =" C "im numpy-Array) die Reihenfolge" [k, j, i] "und in Fortran" (i, j, k) "lautet. (Wird manchmal als Zeilenpriorität bzw. Spaltenpriorität bezeichnet.) Ob auf sie in der Reihenfolge zugegriffen wird, in der sie im Speicher gespeichert sind, oder auf diskrete Weise, wirkt sich auf die Ausführungszeit aus, insbesondere wenn die Größe des Arrays groß ist.
Schreiben wir als Beispiel einen Code, der ein zweidimensionales Array von numpy an Fortran übergibt und ein Doppel davon zurückgibt (dasselbe gilt für drei oder mehr Dimensionen). Die Details des Codes werden hier nicht erläutert. Wenn Sie also beispielsweise Verknüpfung mit Fortran, C-Sprache kennen möchten % E7% B7% A8 / Fortran,% 20C% E8% A8% 80% E8% AA% 9E% 20% E3% 81% A8% E3% 81% AE% E9% 80% A3% E6% 90% BA. Siehe andere Artikel wie 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
Ein Zeiger wird an das Argument der Fortran-Subroutine (Shared Library) übergeben. In Python und Fortran sind die Indizes das Gegenteil von "[j, i]" und "(i, j)", aber es sollte einfach sein, wenn Sie die Reihenfolge berücksichtigen, in der die Daten im Speicher gespeichert werden. ..
Verwenden Sie als Nächstes f2py und schreiben Sie denselben Code wie zuvor. Ich werde auch nicht auf Details eingehen, aber ich denke, Sie können auf Artikel wie Umgang mit Numpy-Arrays mit f2py verweisen.
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.]]
Wenn Sie f2py verwenden, können Sie in Python einfach "[j, i]" als "(j, i)" in Fortran verwenden. Es hat den Vorteil, dass der Code gleich aussieht, aber die do-Anweisung in doubling2.f90 hat den Nachteil, dass die Zugriffsreihenfolge des Arrays nicht der Datenspeicherreihenfolge im Speicher entspricht.
Wie kann ich das lösen? Mit anderen Worten, ich möchte auf den Code auf der Fortran-Seite in derselben Reihenfolge wie doubling1.f90 zugreifen.
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
Sie könnten denken, dass Sie die Speicherreihenfolge im Speicher umkehren können, indem Sie das Argument "order" des numpy-Arrays auf "F" anstelle des Standardarguments "C" setzen. Dies löst dieses Problem jedoch nicht. Das heißt, selbst wenn Sie order =" C " inorder = "F"in doubling2.py ändern, lautet der entsprechende Fortran-Code immer noch doubling2.f90, und die Verwendung von doubling3.f90 führt zu einem Fehler. ist.
Die Lösung besteht darin, das Array mithilfe des Attributs T zu transponieren, wie in f2py - Array-Neuordnung verhindern - Stapelüberlauf beschrieben. Das ist.
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.]]
Wenn Sie beispielsweise GPV mithilfe der netCDF4-Bibliothek lesen möchten, diese effizient mit einer Fortran-Subroutine verarbeiten und dann mit matplotlib + cartopy zeichnen möchten, wird diese Methode empfohlen. Es ist mühsam, umzuziehen, um "Osetsusuke" zu vermeiden, aber selbst wenn es enthalten ist, ist der Umfang der Beschreibung im Vergleich zur Verwendung von numpy.ctypeslib gering, und ich denke, dass f2py einfacher zu verwenden ist.
Recommended Posts