Convertire codice Fortran in Python – 2

Continuo da qui.

Questione preliminare: perché non usare le liste?
Risposta: le liste sono per elementi eterogenei e sono lente e avide di memoria. Se i dati sono omogenei (p.es. tutti int o tutti float o casi simili) ci sono gli array mono- e pluri-dimensionali (matrici).

Python ha un modulo per la gestione degli arrays, ma c’è di meglio: NumPy.

Ogni trattazione di questo argomento inizia con un avvertimento: attenzione alla gestione degli indici, e continua con l’esempio Fortran “per colonna” invece del più comune “per riga”. OK, da verificare. C’è poi un altro aspetto cui prestare attenzione, anzi due. In Fortran il primo indice è 1; in Pthon 0; lo sanno tutti, è banale ma è facile fare confusione, esempio:

$ py3
>>> import numpy as np
>>> ar = np.arange(10)
>>> ar
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> ar.shape
(10,)
>>> ar.size
10
>>> ar[0]
0
>>> ar[5]
5
>>> ar[ar.size]
Traceback (most recent call last):
  File "", line 1, in 
IndexError: index 10 is out of bounds for axis 0 with size 10
>>> ar[ar.size - 1]
9

OK? Quando si ha una matrice attenzione ala gestione degli indici, dicono tutti. Da verificare in Fortran (ebbene sì! ma è semplice e c’è una piacevole sorpresa).

Gestione indici per per riga, file row_col.f

      program row_col

      integer row, col
      real m(1000,1000)

      do row = 1, 1000
          do col = 1, 1000
              m(row, col) = col + 10 * row
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      do row = 1, 1000
          do col = 1, 1000
              m(row, col) = m(row, col) / 2
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      end

Compilo ed eseguo:

$ gfortran -o row-col row-col.f
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,016s
user	0m0,012s
sys	0m0,004s
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,018s
user	0m0,014s
sys	0m0,004s
$ time ./row-col
   11.0000000       220.000000       10001.0000
   5.50000000       110.000000       5000.50000

real	0m0,019s
user	0m0,012s
sys	0m0,008s

Uh! per un milione di dati float, real*4 in Fortran 77, il tempo di esecuzione è troppo piccolo per essere significativo. E questa è la versione lenta, ecco quella corretta, col-row.f

      program col_row

      integer row, col
      real m(1000,1000)

      do row = 1, 1000
          do col = 1, 1000
              m(col, row) = col + 10 * row
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      do row = 1, 1000
          do col = 1, 1000
              m(col, row) = m(col, row) / 2
          end do
      end do

      print*, m(1, 1), m(2, 200), m(1000, 1)

      end

Compilo ed eseguo:

$ gfortran -o col-row col-row.f
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,018s
user	0m0,018s
sys	0m0,000s
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,017s
user	0m0,005s
sys	0m0,013s
$ time ./col-row
   11.0000000       2002.00000       1010.00000
   5.50000000       1001.00000       505.000000

real	0m0,017s
user	0m0,013s
sys	0m0,004s

sì, migliore! di quanto? pochissimo!

Non siamo più ai tempi del PDP-11, nemmeno a quelli del VAX suo successore e la RAM oggi è chip e abbondante. Ecco perché il Web è pieno di foto di gatti, mica come ao tempi di Lenna!

OK, questo per il Fortran, entri Python, con Numpy.

Prima versione col-row.py

#!/usr/bin/python3

import numpy as np

M = np.arange(1000*1000, dtype=float).reshape(1000,1000)
for row in range(1000):
    for col in range(1000):
        M[col, row] = col + 1 + 10 * (row + 1)

print(M[0, 0], M[1, 199], M[999, 0])

for row in range(1000):
    for col in range(1000):
        M[col, row] = M[col, row] / 2

print(M[0, 0], M[1, 199], M[999, 0])

Eseguo:

$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,811s
user	0m0,798s
sys	0m0,012s
$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,807s
user	0m0,781s
sys	0m0,020s
$ time py3 col-row.py
11.0 2002.0 1010.0
5.5 1001.0 505.0

real	0m0,798s
user	0m0,776s
sys	0m0,016s

Seconda versione row-col.py

#!/usr/bin/python3

import numpy as np

M = np.arange(1000*1000, dtype=float).reshape(1000,1000)
for row in range(1000):
    for col in range(1000):
        M[row, col] = col + 1 + 10 * (row + 1)

print(M[0, 0], M[1, 199], M[999, 0])

for row in range(1000):
    for col in range(1000):
        M[row, col] = M[row, col] / 2

print(M[0, 0], M[1, 199], M[999, 0])
$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,771s
user	0m0,755s
sys	0m0,016s
$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,823s
user	0m0,800s
sys	0m0,024s
$ time py3 row-col.py
11.0 220.0 10001.0
5.5 110.0 5000.5

real	0m0,794s
user	0m0,774s
sys	0m0,020s

sì, Python è lento ma ampiamente accettabile.

OK, ho tutto quello che serve per la conversione; con Python sarà semplice aggiornare e aggiungere nuove funzionalità, i grafici, ma non solo. Per esempio capita (ahemm, sempre) di dover fare la matrice trasposta, quella con righe e colonne scambiate. A me che sono vecchio vedere come fa NumPy vengono le vertigini:

$ py3
>>> import numpy as np
>>> M = np.arange(12).reshape(3,4)
>>> M
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> M.T
array([[ 0,  4,  8],
       [ 1,  5,  9],
       [ 2,  6, 10],
       [ 3,  7, 11]])

Il manuale di NumPy è qui.

OK? vero che Python è (quasi) l’unico linguaggio in città. O almeno quello più sexy 🤩

Posta un commento o usa questo indirizzo per il trackback.

Rispondi

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo di WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione /  Modifica )

Google photo

Stai commentando usando il tuo account Google. Chiudi sessione /  Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione /  Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione /  Modifica )

Connessione a %s...

Questo sito utilizza Akismet per ridurre lo spam. Scopri come vengono elaborati i dati derivati dai commenti.

%d blogger hanno fatto clic su Mi Piace per questo: