Category Archives: Fortran

Il linguaggio di programmazione Fortran

L’avvenire del Fortran e i suoi successori

fortran

Forse mi sono lasciato prendere la mano, non so se sia così necessario questo post. Ma siccome l’ho scritto lo pubblico :???:
Tutto è partito da qui: Why Scientists Are Still Using FORTRAN in 2014 che rimanda a Ars technica, qui: Scientific computing’s future: Can any coding language top a 1950s behemoth? dove troviamo

But almost universally, the language in which these simulation codes are written is Fortran, a relic from the 1950s.

e

Wherever you see giant simulations of the type that run for days on the world’s most massive supercomputers, you are likely to see Fortran code.

Poi dice che subito dopo viene inventato il Lisp. Che non ha successo per la sua stranezza:

Weirdness, because the prefix notation used by Lisp made expressions in that language look a lot less like normal mathematics than did math rendered in Fortran. (Fortran stood, after all, for FORmula TRANslator.) A normal chemist or engineer is far more comfortable with y = (a + b)/c than with (setf y ((/ (+ a b) c))).

Ok, l’ho sentita tante volte che mi verrebbe voglia di ignorarla. Invece no! Perché è (quasi) tutto lì. Davvero, ve ne faccio il riassunto.

(setf y (/ (+ a b) c))

Parto dalla parentesi più interna: (+ a b); questa è una lista il cui primo elemento è una funzione. Questa funzione opera sui successivi elementi della lista e li somma. Questa somma (chiamiamola A+B viene restituita: Lisp restituisce sempre l’ultima espressione calcolata.
Quindi la nostra espressione diventa:

(setf y (/ A+B c))

dove (/ A+B c) è una lista in cui il primo elemento è la funzione / che restituisce il rapporto del secondo elemento della lista per tutti i successivi.

Il valore restituito viene assegnato alla variabile y, tramite setf (setf sta per set (quote e sì, si può scrivere anche setq o set ' –non sono esattamente la stessa cosa ma il risultato è lo stesso).

Uh! visto? il Lisp è uniforme. Si comporta sempre allo stesso modo: il primo elemento della lista (chiamato CAR o first) indica come trattare il resto della lista (CDR o rest) e restituisce il risultato.
Per cui il Lisp è fatto da liste dentro liste dentro liste ancora.  Tutto lì (ok, sto semplificando). Non so voi ma a me sembra molto semplice.

Per contro il Fortran riproduce la sintassi matematica che si usa a scuola:

y = (a + b) / c

Sì, lo ammetto, siamo abituati così e per un calcolo semplice come questo meglio il Fortran. Ma per un caso non banale si vede la differenza: con il Lisp ho una concatenazione di liste e con il Fortran una serie di assegnazioni (attenzione il segno = si deve leggere “diventa“, dice il Wirth (che scrive :=)).

11c

Poi una cosa ancora: io sono vecchio e non so nemmeno se si usano ancora ma quando frequentavo il Poli, roba anni ’70, i migliori usavano le calcolatrici HP, con RPNreverse polish notation– quasi come il Lisp (cioè no, esattamente il contrario), nel nostro caso si sarebbe scritto

a
b
+
c
/

Con risparmio di digitalizzazione di tasti –non servono le parentesi–, si consumano meno i diti e –sopratutto– i risultati intermedi sono disponibili nello stack (mi dicono che dovrei usare il termine “pila”).

Ma tornando alla programmazione allora quando è uscito il Fortran era una cosa nuova di pakka, un compilatore per linguaggio facile per gli umani, alternativo alle criptiche istruzioni assembly. E sui calcolatori computer c’era solo il Fortran. Molto diverso da quello di oggi, istruzioni molto elementari, quasi da assembly. E tanti GOTO, spesso mascherati come IF.

Ok, fine dello sfogo, torno al post di Ars Technica secondo il quale ci sono tre candidati a sostituire il Fortran:

Haskell, bello, sintassi da matematici, funzionale (cioè odia le variabili), chissà;
Clojure, Lisp & JVM, sì c’è Java sotto; ho scritto qualcosa in passato;
Julia, nuovissimo, ne so molto poco.

Ma ce ne sono altri, secondo me: il C++, Java, Python e altri ancora. Tutti ampiamente utilizzati.
E i nuovissimi Go e Rust. Di Go sono a conoscenza di un programma funzionante, per il Web, Rust promette bene (ne parla spesso robitex su questo blog).
Inoltre proprio ieri vosto questo post: #golang: The Next Great Teaching Language. Il blog è nuovo ma promette bene. E il post è illuminante.

Invece il Fortran non è amato dai giovani (quelli che conosco), lo trovano vecchio, inusuale. E se possono lo sostituiscono con qualcosa di più pratico, per esempio Matlab (a volte basta la versione free Octave). Ma sto parlando di chi usa il computer come un mezzo complementare. E lo usa per AutoCAD, fogli di calcolo Excel e come macchina da scrivere.

Mi mancano i programmatori di una volta, quando il mondo era nuovo ;-)

Precisione iii

Corto circuito

Questo post è un bellissimo corto circuito: conosco Orlando perché è un amico del GuIT che tra l’altro ho conosciuto l’anno scorso di persona al GuIT meeting di Napoli. E proprio Orlando ha commentato il post di Juhan Precisione su questo stesso blog, di qualche giorno fa su alcuni calcoli.
Questo dimostra che gli appassionati non vedono l’ora di scoprire nuove cose sul proprio argomento preferito.

Chiudo il circuito

Eccomi dunque a chiudere il circuito informatico presentrando un programma in Go che esegue le stesse operazioni eseguite in awk ed in Fortran:

package main

import (
    "fmt"
    "math/big"
)

func main() {
    s, u := big.NewInt(2), big.NewInt(1)
    for i := 1; i <= 10; i++ {
        // a = s - 1
        a := big.NewInt(0).Set(s).Sub(s, u)
        // s*a + 1 = s*(s - 1) + 1
        s.Mul(a, s).Add(s, u)
        fmt.Println(i, s)
    }
}

Tecnicamente, vorrei farvi notare che è possibile in Go concatenare le chiamate dei metodi come per esempio si legge alla riga 12 del codice sorgente, se il metodo stesso ritorna lo stesso oggetto.

Ho utilizzato il pacchetto dei grandi numeri big necessario perché dopo la sesta iterazione i numeri superano il limite massimo del tipo int64. Certo avrei potuto utilizzare il tipo uint64 ma potevo chiudere il circuito soltanto esagerando alla grande portando il ciclo fino a 10 iterazioni… ecco di seguito il risultato:

1 3
2 7
3 43
4 1807
5 3263443
6 10650056950807
7 113423713055421844361000443
8 12864938683278671740537145998360961546653259485195807
9 165506647324519964198468195444439180017513152706377497841851388766535868639572406808911988131737645185443
10 27392450308603031423410234291674686281194364367580914627947367941608692026226993634332118404582438634929548737283992369758487974306317730580753883429460344956410077034761330476016739454649828385541500213920807

A questo punto, ci sta propio bene un
Alla Prossima!
R.

Precisione – 2

precisione
Il post di ieri, Precisione è stato visitatissimo! E Orlando mi chiede una cosa che è difficile far stare in un commento. E allora invece di un commento ecco un intero post.

Orlando pubblica un programma in Fortran questo:

program test
  real :: s=2.0
  do i=1,7
    s=s*(s-1)+1
    write(*,10) s
  end do
10 format(F30.1)
end program test

che produce questo risultato:

                           3.0
                           7.0
                          43.0
                        1807.0
                     3263443.0
              10650057179136.0
 113423716646946804535918592.0

simile ma non uguale a quelli degli script da me pubblicati.
Allora, vediamo. Comincio con modificare lo script di bc, facendogli scrivere i risultai dei singoli passi del ciclo, ecco:

s = 2.0
for (i = 1; i <= 7; i++) {
    s = s * (s - 1) + 1
    print s, "\n"
}
quit

3.0
7.0
43.0
1807.0
3263443.0
10650056950807.0
113423713055421844361000443.0

Dove si vede che le difformità sono negli ultimi due passi, questi:

3.0
7.0
43.0
1807.0
3263443.0
* 10650056950807.0
F 10650057179136.0
         ^ 
* 113423713055421844361000443.0
F 113423716646946804535918592.0
          ^

dove * rappresenta l’output di bc, F quello del Fortran e ^ la prima cifra diversa tra i due. Si vede subito che il numero di cifre per i linguaggi normali (anche il C, dice infatti Orlando) sono 7-8. Più che sufficienti in condizioni normali come concludevo nel post di ieri.
Naturalmente si può fare meglio. In particolare per il programma Fortran ditato sopra basta modificare la seconda riga, così:

real(16) :: s=2.0

e ottengo lo stesso risultato per Fortran e bc (e gli altri programmi, presumibilmente, anzi certamente).

In passato (l’anno scorso) ho raccontato un po’ di cose sul Fortran (per via che si stava facendo manutenzione a cose preistoriche), si trova tutto nell’indice di OK, panico, cercare per Fortran.

In particolare, siccome non ricordo certo cosa ho scritto a giugno del 2012 consiglio una scorsa qui e qui.

Resto disponibile per modifiche e chiarimenti, anzi questo post possiamo considerarlo un prepost ;-)
E, un’ultima cosa (Steve J. finiva sempre così), lo sapete che Orlando ha un sito tutto suo? Questo: Orlando Iovino’s Personal Website. Da finire, certo ;-)

f2py – usare codice Fortran in Python

jakeTutto è partito da qui: Numba vs. Cython: Take 2.

Jake VanderPlas è un tipo tosto, fareste bene a seguirlo. Si presenta così: “I’m an NSF post-doctoral fellow at University of Washington, where I research and teach in Astronomy, Astrophysics, and Machine Learning“.

Jake accenna a f2py. Ehi! potrebbe essere utile a chi ha vecchia roba scritta in Fortran, magari si riesce a metterle un’interfaccia grafica o fare altre mescite di quelle che, fa … beh! Chissà, vediamo.

Intanto salta fuori che f2py è un programma a se stante che compila il codice Fortran (usando il compilatore che hai, dev’essercene uno) creando una libreria (.so), quelle che Python chiama moduli.
Ma è anche un modulo di NumPy, cioè Python. Ed è tutto scritto in Python.
Allora la prima cosa da fare (è da tanto che mi dicevo di farlo) è installare NumPy: f2py è lì. Installare NumPy richiede 952MB di spazio su disco! OK, c’è tempo per un caffè lungo.

Fatto, viene installato anche SciPy, ne ho già parlato in passato. Per vedere comincio a riprodurre l’esempio preso dalla documentazione, questa è la subroutine in Fortan 77:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
C FILE: FIB1.F
      SUBROUTINE FIB(A,N)
C
C     CALCULATE FIRST N FIBONACCI NUMBERS
C
      INTEGER N
      REAL*8 A(N)
      DO I=1,N
         IF (I.EQ.1) THEN
            A(I) = 0.0D0
         ELSEIF (I.EQ.2) THEN
            A(I) = 1.0D0
         ELSE 
            A(I) = A(I-1) + A(I-2)
         ENDIF
      ENDDO
      END
C END FILE FIB1.F

Si crea il modulo Python con il comando:

f2py -c fib1.f -m fib1

e questo è il file di Python:

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy as np
import fib1 #quello in Fortran

a = np.zeros(10, 'd')
fib1.fib(a)
print a

f1

OK, vediamo se si può migliorare: Fibonacci restituisce interi, anzi ai tempi di Leonardo Pisano (Fibonacci per i matematti) i float non erano conosciuti. E meno che mai in doppia precisione (come si dice in Fortran).

Le modifiche sono minime:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
C FILE: FIB2.F
      SUBROUTINE FIB(A,N)
C
C     CALCULATE FIRST N FIBONACCI NUMBERS
C versione modificata, uso interi
      INTEGER N
      INTEGER A(N)
      DO I=1,N
         IF (I.EQ.1) THEN
            A(I) = 0.0D0
         ELSEIF (I.EQ.2) THEN
            A(I) = 1.0D0
         ELSE 
            A(I) = A(I-1) + A(I-2)
         ENDIF
      ENDDO
      END
C END FILE FIB2.F

che compilo con

f2-0

f2py riempie il terminale di scritte, roba da preoccuparsi se uno non sapesse che è normale. Ecco lo script Python

#!/usr/bin/python
# -*- coding: utf-8 -*-

#modificato, uso interi

import numpy as np
import fib2 #quello in Fortran

a = np.zeros(10, 'int')
fib2.fib(a)
print a

f2-1

OK. Vediamo qualche caso tipico, non quelli presenti nei tutorial, ma quelli che so che possono capitare, per esempio un “Hello World!” poco funzionale, chissà se funziona.

1
2
3
4
5
      subroutine hw
c scrive il classico saluto

      print *, "Hello World!"
      end

Ho dovuto usare la sintassi vecchia, quella fixed-form, in attesa di leggere il manuale. Ma un passo per volta.

#!/usr/bin/python
# -*- coding: utf-8 -*-

import hw
print "Adesso chiamo il Fortran"
hw.hw()
print "Fatto, OK?"

hw

OK. La pessima abitudine dei programmatori di una volta era di scrivere su file quando capita (visto che mi sto funzionalizzando?), vediamo:

1
2
3
4
5
6
7
8
      subroutine hw2
c scrive il classico saluto

      print *, "Hello World!"
      open(1, file = "prova-f2py.txt")
      write(1, '("Hello World!")')
      close(1)
      end

Per sicurezza ho controllato che con un main in Fortran funzionasse, non si sa mai… È OK.

#!/usr/bin/python
# -*- coding: utf-8 -*-

import hw2
print "Adesso chiamo il Fortran"
hw2.hw2()
print "Fatto, OK?"

hw2

OK, continua a funzionare. E se uno ha una subroutine che ne chiama un’altra in un altro file cosa succede?

1
2
3
4
      subroutine hw3

      call hw2  
      end

La hw2 è la stessa del caso precedente. Il comando per f2py è:

f2py -c hw3.f hw2.f -m hw3

e poi:

#!/usr/bin/python
# -*- coding: utf-8 -*-

import hw3
print "Adesso chiamo il Fortran"
hw2.hw2()
print "Fatto, OK?"

hw3

Chissà se posso passare una stringa al Fortran?

1
2
3
4
5
6
      integer function fs(st)
      character*(*) st

      print *, st
      fs = len(st)
      end
#!/usr/bin/python
# -*- coding: utf-8 -*-

import fs

msg = "f2py rockz!"
print "Adesso chiamo il Fortran"
res = fs.fs(msg)
print res
print "Fatto, OK?"

fs

OK! Un’ultima cosa, poi mi leggo tutto il manuale ma per intanto il free-form, non so se non si possa semplificare, questo è quello che sono riuscito a ottenere finora:

1
2
3
4
5
6
7
8
9
module fhw
  integer i
contains
  subroutine hw
    ! scrive il classico saluto

    print *, "Hello World!"
  end subroutine hw
end module fhw
#!/usr/bin/python
# -*- coding: utf-8 -*-

import hw
print "Adesso chiamo il Fortran"
hw.hw()
print "Fatto, OK?"

hw

No! è solo questione di mettere l’estensione .f90 e funziona tutto:

1
2
3
4
5
subroutine hw
! scrive il classico saluto

print *, "Hello World!"
end subroutine hw
#!/usr/bin/python
# -*- coding: utf-8 -*-

import hw
print "Adesso chiamo il Fortran"
hw.hw()
print "Fatto, OK?"

L’output è uguale al precedente.

Quasi dimenticavo:
NumPy è opera di Travis E. Oliphant, con il cotributo di numerosi collaboratori, vedi Origins of NumPy in Guide to NumPy.
f2py è stato scritto da Pearu Peterson.

Bellissima (anche se l’ho solo scorsa per adesso) la Guide to Numpy, se si ha fretta c’è anche An introduction to Numpy and Scipy. E poi c’è il Cookbook, il sito di SciPy dove ci sono ovviamente NumPy e f2py.

Quindi, conclusione provvisoria: si può usare? Forse, dipende. Leggo il manuale e vi dico.
Prossimamente… forse… ;-)

Sul FORTRAN di EWD

f0A seguito del post su Dijkstra  e di un paio di telefonate mi è venuta voglia di farvi vedere com’era il codice quando ho iniziato a lavorarci io. In realtà sto barando: con il Fortran 77 l’if era ormai multilinea, come if/then/else ma i vecchi (di allora) continuavano a usare il IV e lo imponevano.

Sieve_of_Eratosthenes_animation

L’immagine che illustra il crivello di Eratostene che in questi giorni circola su G+ (credo sia di John Baez Skopp (vedi commento, grazie Annarita!)) è l’ideale per implementare un programmino nel Fortran delle origini. Eccolo

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
C     ERAHOSTHENES

      DIMENSION NUMERI(20), LPRIMI(10)

C     NUMERI CONTIENE I NUMERI CHE POSSONO ESSERE PRIMI  
      DO 10, I = 1, 20
10    NUMERI(I) = I
      NUMERI(1) = 0  

C     NPRIMI = IL NUMERO DI PRIMI TROVATI
      NPRIMI = 0
      NTROV = 2

C     QUANDO NE TROVO UNO LO AGGIUNGO ALLA LISTA DEI PRIMI
20    NPRIMI = NPRIMI + 1
      LPRIMI(NPRIMI) = NTROV

C     ELIMINO I SUOI MULTIPLI NELLA LISTA DEI NUMERI
      J = 1 
30    J = J + 1
      NT = NTROV * J
      IF(NT .GT. 20) GOTO 40
      NUMERI(NT) = 0
      GOTO 30

C     PASSO ALL ESAME DEL CANDIDATO SUCCESSIVO
40    N = NTROV
41    N = N + 1
      IF(N .GT. 20) GOTO 60
      IF(NUMERI(N) .EQ.0) GOTO 41
50    NTROV = N
      GOTO 20

C     FINE DELL ESAME DELLA LISTA
60    DO 65 I = 1, NPRIMI
65    WRITE(*, 90) I, LPRIMI(I)
90    FORMAT(2I4)

      END

ERAT

Nello spirito del tempo ho scritto tutto in maiuscolo ;-) Però non ce l’ho fatta a non lasciare spazi: allora si sarebbe scritto tutto come J = J + 1. E un’altra cosa ancora: la riga a cui l’esecuzine saltava a seguito di un GOTO era tipicamente del tipo:

50    CONTINUE

essendo quest’istruzione segnaposto fatta proprio apposta per questo. Ma il mio programma sarebbe diventato troppo lungo, inutilmente.
Un’altra cosa: è quasi immediata l’ottimizzazione considerando che 2 è l’unico primo pari, dimezzando il numero di numeri da testare ma il codice si allunga, lasciato come esercizio.

Ecco, allora si programmava così. Non ho più il libro ma nel testo che si usava al Poli c’era scritto: “finito di perforare le schede fai una stampa del listato e controllala attentamente prima di provare l’esecuzione“, vero anche perché a me al primo tentativo è successo questo:

0

Avevo sbagliato il numero di una label.
Bon! Nostalgia OFF! :-D

lcm – Fortran e linguaggi normali

Nuova puntata della telenovela lcm, le precedenti sono qui e qui.
Perché ci sono un paio di cose nuove (per me almeno) da dire sull’argomento. Roba che non finisce oggi, forse, credo.

lcm-matFinora ho utilizzato linguaggi che rendono il problema facile, anche se non hanno una core function che lo implementa.

Mi interessava vedere i linguaggi normali, tipo C/C++, Java, Go, Io che sono vecchio mi sono detto: e il Fortran? Recentemente ha subito una cura ringiovanente, è diventato come gli altri, indistinguibile (quasi).
Verificare è facile, c’è Google. Ho trovato questo: FORTRAN PROGRAMS FOR SOLVIMG NUMERICAL PROBLEMS, Designed by T K Rajan.
Bello, ci sono tutte (o parecchie) delle cose che servivano e che ogni programmatore Fortran aveva una volta. c’è anche lui, l’lcm, a p.7. Vediamo:

PROGRAM gcdlcm
integer a,b,l

write(*,*)'To find the gcd and lcm of two numbers.'
write(*,*)'Enter 2 numbers:'
read(*,*)a,b

m=a
n=b
IF(a.GT.b)THEN
    j=a
    a=b
    b=j
ENDIF
10  i=mod(b,a)
    IF(i.EQ.0)THEN
        write(*,20) a
20      format(1x,'GCD= ',I4)
    ELSE
        b=a
        a=i
        GOTO 10
    ENDIF

l=m*n/a
write(*,30)l
30 format(1x,'LCM= ',I4)

END

l0

Prima delusione: è scritto come si usava una volta, diciamo Fortran 77, 1978 e via fino al ’90. Provo a attualizzarla un pochino, ecco la mia versione:

program lcm

implicit none
integer a, b, i, l, m, n

print *, 'Enter 2 numbers:'
read(*,*) a, b

m = a
n = b
if(a > b) then
    i = a
    a = b
    b = i
end if

i = mod(b, a)
do while(i /= 0)
    b = a
    a = i
    i = mod(b, a)
end do
print *, 'GCD = ', a

l = m * n / a
print *, 'lcm = ', l

end program

Però, cosa più grave, siamo sempre limitati a due numeri. E non abbiamo l’equivalente di reduce() o apply(). Forte di quello visto in precedenza mi son detto: “si può fare”. Anzi è risultato semplicissimo, un ciclo for do (in Fortran il for si chiama do), righe 14-17.

program lcm_ftn

implicit none

integer :: i, numero_dati, dati(10)
integer :: lcm, lcm_2

call leggi_input(numero_dati, dati)
if (numero_dati < 2) then
    print *, 'dati insufficienti'
    call exit(1)
end if

lcm = dati(1)
do i = 2, numero_dati
    lcm = lcm_2(lcm, dati(i))
end do

print '(a6, i4)', 'lcm = ', lcm

end program

subroutine leggi_input(narg, arg)
implicit none

character(20) :: argv ! argomenti linea di comando
integer :: narg ! numero argomenti
integer :: i, t, arg(10)

narg = command_argument_count()
do i = 1, narg
    call get_command_argument(i, argv)
    read(argv, *), arg(i)
end do

end subroutine

integer function lcm_2(a, b) result(res)
implicit none

integer :: a, b, i, m, n

m = a
n = b
if(a > b) then
    i = a
    a = b
    b = i
end if

i = mod(b, a)
do while(i /= 0)
    b = a
    a = i
    i = mod(b, a)
end do

res = m * n / a

end function

lcmftn

OK, fatto. Semplice vero? E si capisce perché non ci sia nei linguaggi normali. Avevo pensato di rifarlo anche in awk (come Bit3Lux sono affascinato da awk, e io ho cominciato prima) ma poi ho pensato che era solo una trascrizione per l’algoritmo complicato dalla necessità di usare un linguaggio nato per script one-liner (o poco più). Nel mio caso le righe di codice sarebbero state parecchie.

C’è anche un altro motivo: sapete com’è, con Google, a volte, quando meno te l’aspetti…
Ma non voglio anticipare l’argomento della prossima puntata (anche se JeanMM…) ;-)

Non so se avete notato che ho sempre parlato di righe e mai di schede (almeno credo di averle corrette tutte).

Fortran, common o moduli?

Ecco la seconda puntata della telenovela, iniziata qui.


È possibile continuare a usare i common per le variabili globali o conviene passare ai moduli?

Dal punto di vista prestazionale non ci sono differenze, come vedremo tra breve, ma dal lato logico le cose cambiano. Vediamo un piccolo esempio.

Ecco la versione che utilizza i common (oldc):

program oldc
implicit none

call fillmat
call showval

end

subroutine fillmat
implicit none

integer, parameter :: nmax = 10000
integer :: mat(nmax, nmax)
common /bigmat/ mat

integer :: r, c

do c = 1, nmax
    do r = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

end

subroutine showval
implicit none

integer, parameter :: nmax = 10000
integer :: mat(nmax, nmax)
common /bigmat/ mat

print *, mat(nmax, 1), mat(nmax, nmax), mat(1, nmax)

end

Notare che il parameter nmax non si può mettere nel common, anche se si riferisce a questo. Converrebbe mettere in un file include le righe [12..14], ripetute a [33..35] per evitare errori in aggiornamento.

Lo stesso programma riscritto utilizzando un modulo diventa:

program newm
use bigmatm

implicit none

call fillmat
call showval

end program

 

module bigmatm
implicit none

integer, parameter :: nmax = 10000
integer :: mat(nmax, nmax)

contains
    subroutine fillmat
        integer :: r, c

        do c = 1, nmax
            do r = 1, nmax
                if (r >= c) then
                    mat(r, c) = r * c
                else
                    mat(r, c) = 0
                end if
            end do
        end do

    end subroutine fillmat

    subroutine showval

        print *, mat(nmax, 1), mat(nmax, nmax), mat(1, nmax)

    end subroutine showval

end module bigmatm

Per comodità il modulo è su un file separato (cosa normale per i programmi veri, molto più lunghi di questo esempio minimale). In questo caso le caratteristiche della matrice sono specificate una sola volta, diminuisce la possibilità di errori. Inoltre le due subroutine risultano private al modulo (e a chi lo usa).

Per contro le due versioni sono equivalenti rispetto al tempo di esecuzione. Io sostengo l’uso dei moduli, in qualche misura cambia anche l’impostazione logica del problema. Per contro e possibile mischiare codice vecchio, con common, definizione implicita delle variabili, formato fisso e quant’altro con il nuovo standard. Ma quello che si deve ancora scrivere, dai, facciamolo corrente.

Ci sarebbe poi un’ulteriore questione, alla quale non so rispondere, anzi benvenuti i suggerimenti. Si è deciso di usare il Fortran per un paio di motivi: c’è già parecchia roba scritta in Fortran che gira qui e poi perché il Fortran è abbastanza simile al Basic (una volta si diceva che il Basic era un Fortran semplificato). Poi è gratis. Il problema è che per averlo su Windows occorre installare Cygwin e pare non sia così immediato. Dai test effettuati abbiamo visto che Java è veloce e la JVM c’è dappertutto. Allora non converrebbe fare tutto in Java?

Ordine per cicli nidificati e trasposizione matrici

Credevate che fosse finita la telenovela del Fortran, invece no!
Di ritorno dalle ferie si riprende con due puntate: * 2 puntate 2 *, ma brevissime. Questa è la prima.

La questione è un classico ma siccome salta fuori continuamente lo ripeto ancora una volta. L’ultima! (A meno che…).

Supponiamo di avere una matrice bidimensionale, per esempio punti X, Y del piano cartesiano, come conviene gestirli?
Il Fortran, a differenza degli altri linguaggi (C/C++, Basic, Python, tutti) ordina i dati per colonna. Fintanto che i dati sono pochi la differenza tra la versione corretta e quella usuale è oggi, con le macchine attuali, trascurabile. Diventa però sensibile e poi non è bello.

Ecco un esempio, minimo, didascalico. Talmente didascalico che lo riciclerò tutte le volte che la questione salta fuori. (Capita sovente dopo le ferie).

Le due versioni seguenti producono lo stesso risultato: rempiono la parte sopra la diagonale principale (compresa) con la somma degli indici di riga e colonna. Con una differenza: la versione per colonne è 3.2 volte più veloce!

program righe
implicit none

integer, parameter :: nmax = 10000
integer :: r, c
integer :: mat(nmax, nmax)

do r = 1, nmax
    do c = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

! per controllo
print *, mat(nmax, 1), mat(nmax, nmax), mat(nmax, 1)

end

program colonne
implicit none

integer, parameter :: nmax = 10000
integer :: r, c
integer :: mat(nmax, nmax)

do c = 1, nmax
    do r = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

! per controllo
print *, mat(nmax, 1), mat(nmax, nmax), mat(1, nmax)

end

OK, si tratta di secondi ma voglio solo vedere quella giusta. Avvisati nèh!

Però metti che io la mia matrice la voglio usare in un programma scritto, per esempio, in C++; devo leggere i dati in modo innaturale?
No! il Fortran ha pensato anche a te e a questo tuo caso: esiste la funzione transpose(). Visualizzarne il funzionamento con una matrice così grossa sarebbe problematico, diciamo che la riduciamo a 10 * 10, ecco:

program inverti
implicit none

integer, parameter :: nmax = 10
integer :: r, c
integer :: mat(nmax, nmax), inv(nmax, nmax)

do c = 1, nmax
    do r = 1, nmax
        if (r >= c) then
            mat(r, c) = r * c
        else
            mat(r, c) = 0
        end if
    end do
end do

inv = transpose(mat)

print '(10i4)' , mat
print *
print '(10i4)' , inv

end

Nota: non è possibile sovrascrivere direttamente la matrice, quest’istruzione è errata: mat = transpose(mat).

Prossimamente la seconda (forse ultima) puntata, altrettanto –come si dice– ecco ;-) :-o :-D

Ottimizzare il codice

We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil.

Questa frase è stata attribuita a Donald E. Knuth, e se lo dice lui che ha scritto “The Art of Computer Programming”, ha inventato il TeX, e ha sviluppato la teoria della complessità degli algoritmi, potete stare certi che sia vero.

Il significato della frase è che bisognerebbe evitare di ottimizzare il codice prima di aver finito di scriverlo. Il modo di programmare “corretto” sarebbe:

  • primo, si scrive il codice nella maniera più semplice possibile, possibilmente in maniera elegante e leggibile
  • poi ci si accerta della sua correttezza (testing , debugging, e quello che pare a voi)
  • Infine, pensiamo alle performance, ottimizzando il codice dove necessario

Spesso invece il programmatore comincia fin dall’inizio a complicare il codice e la strrutra del programma, alla ricerca del programma già perfettamente ottimizzato. Nel fare così, spesso viene fuori codice molto complicato e probabilmente pieno di bug. Da qui la sentenza di Knuth: l’ottimizzazione prematura è la radice di tutti i mali!

L’osservazione chiave è che spesso il codice contiene alcuni (pochi) “colli di bottiglia”, ed è lì che l’ottimizzatore deve concentrarsi per far andare veloce il suo programma. Come identificare questi “colli di bottiglia”? Per esempio, funzioni chiamate molto frequentemente, strutture dati usate pesantemente in tutto il codice, ecc. I programmatori esperti ormai lo fanno ad occhio, ma per programmi molto grandi e complessi non è affatto facile. Sarebbe bene ci fosse qualche strumento di ausilio al programmatore.

oggi vediamo come fare analisi delle prestazioni con due tool disponibili su Linux con licenza GPL: valgrind e kcachegrind.

Il primo è una specie di “emulatore” del vostro PC. Prende in ingresso un file eseguibile, e lo esegue su un “processore virtuale”. Nel portare avanti l’esecuzione, però, è in grado di effettuare una serie di analisi sul codice. Valgrind è una piattaforma generica su cui sono stati sviluppati diversi strumenti. Il più famoso è memchek che serve a controllare l’uso della memoria, ed individuare possibili “memory leak“. Quello di cui ci occupiamo oggi, invece, è callgrind che serve a stimare le performance delle varie parti di codice del nostro programma.

Per installare questi tool, su Ubuntu si deve scrivere sul terminale:

sudo apt-get install valgrind
sudo apt-get install kcachegrind

e siamo pronti per i nostri esperimenti.

Proveremo ad utilizzarlo sui programmini che Juhan ci ha proposto per fare confronti fra il Fortran e vari altri linguaggi. Cominciamo dalla versione C/C++ di tale programma, che riporto qui per comodità (anche perché ho fatto due modifiche).

#include <stdio.h>
#include <math.h>

double f[14];

#define POW(x,y) pow(x,y)
//#define POW(x,y) mypow(x,y)

inline double mypow(double x, int y)
{
    double r = 1;
    for (int i=0; i<y; i++) r *= x;
    return r;
}

double fsin(double a) {
    return a - POW(a, 3) / f[3] + POW(a, 5) / f[5]
        - POW(a, 7) / f[7] + POW(a, 9) / f[9]
        - POW(a, 11) / f[11] + POW(a, 13) / f[13];
}

static double fcos(double a) {
    return 1.0 - POW(a, 2) / f[2] + POW(a, 4) / f[4]
        - POW(a, 6) / f[6] + POW(a, 8) / f[8]
        - POW(a, 10) / f[10] + POW(a, 12) / f[12];
}

static double myln(double x) {
    if (x == 0.0) {
        return -1.0e20;
    } else {
        double r = (x - 1.0) / (x + 1.0);
        return 2.0 * (r + POW(r, 3) / 3.0
                      + POW(r, 5) / 5.0
                      + POW(r, 7) / 7.0
                      + POW(r, 9) / 9.0);
    }
}

double mylog(double x) {
    static double ln10 = myln(10.0);

    return x / ln10;
}

int main(int argc, char **argv)
{

    int i, j;
    f[0] = 1.0;
    for (i = 1; i < 14; i++)
        f[i] = i * f[i - 1];
    int deg = 60;// * 60;
    int nsteps = 180 * deg;
    double step = M_PI / nsteps;
    double ssum;
    double csum;
    double a, s, c, t = 0.0;
    double ls, lc;
    for (j = 1; j < 11; j++)
    {
        ssum = 0.0;
        csum = 0.0;
        for (i = 0; i < nsteps; i++)
        {
            a = i * step;
            s = fsin(a);
            ls = mylog(myln(s));
            ssum += s;
            c = fcos(a);
            lc = mylog(myln(c));
            csum += c;
            if ((i % (10 * deg)) == 0)
            {
                if (c != 0.0)
                    t = s / c;
                printf("%3d %11.8f %11.8f %11.8f %15.8e\n", (i / deg), a, s, c, t);
                printf(" %15.8e %15.8e\n", ls, lc);
            }
        }
    }
    printf("%g\n", ssum);
    printf("%g\n", csum);
}

Le due modifiche fatte sono le seguenti: ho sostituito tutte le chiamate a pow() con la macro POW. Quest’ultima si traduce per il momento nella chiamata di libreria pow() (che, come discusso nel post precedente, è più lenta). Inoltre, ho diminuito il numero di cicli da svolgere per accorciare un po la durata dell’esecuzione: adesso il conto viene fatto per deg = 60 invece che per deg = 60*60 (vedi linea 54).

Salvate il programma qui su sul vostro HD con il nome prova.cpp, e poi compilatelo con il seguente comando:

g++ -g prova.cpp -o provacpp

Il flag -g serve per inserire i simboli di debug nel file eseguibile, e la qual cosa ci permetterà di associare il costo delle operazioni alle relative istruzioni nel codice sorgente.

Adesso siamo pronti per l’esecuzione, ecco il comando da dare:

valgrind --tool=callgrind ./provacpp

Questo manda in esecuzione il tool callgrind sul file eseguibile provacpp. Il quale viene eseguito a una velocità molto inferiore a quella normale (parliamo di un rallentamento di quasi 100 volte!) e vengono collezionati i dati di performance. L’output è un file che si chiama

callgrind.out.PID

Dove PID è l’ID del processo che è stato eseguito.
Per visualizzare i dati, date il comando:

kcachegrind callgrind.out.PID

(dove naturalmente avrete sostituito PID con il numero giusto), e si aprirà una finestra come questa:

La finestra principale di KCacheGrind

Questo tool ci visualizza i risultati dell’esecuzione. Faccio zoom sulla finesta a sinistra dove viene mostrato il “call graph” (ovvero il grafo delle chiamate):

Come vedete, il main() chiama le funzioni myln, fsin e fcos, e ciascuna chiama la funzione pow, la quale a sua volta chiama le funzioni interne per fare operazioni floating point. Subito sotto il nome delle funzioni viene riportato il tempo (in percentuale) speso dal programma dentro quella funzione. Questo programma ha speso il 95.73% del suo tempo dentro la funzione pow()! Notare anche una certa asimmetria: myln prende il 50% del tempo, mentre fsin e fcos ciascuna circa il 25%.

Adesso riproviamo utilizzando la funzione mypow() (quella che fa le moltiplicazioni per fare le potenze): basta commentare la linea 6 e scommentare la 7, e ricompilare. Se rifate tutto il percorso otterrete questo callgraph:

Come vedete, le cose sono un po’ cambiate. In particolare, mypow() non chiama alcuna funzione in floating point (come ci si aspetta), e inoltre le percentuali di myln, fsin e fcos sono un po’ cambiate.

È anche possibile mostrare questi numeri direttamente sul sorgente, cliccando sulla tab “Source Code”, ed ecco lo snapshot relativo:

Infine: si può fare una cosa del genera anche per il fortran? Ovviamente sì! Vi riporto semplicemente il comando di compilazione e di esecuzione:

gfortran -g m_c10l.f08 fc10l.f08 -o provaf
valgrind --tool=callgrind ./provaf

E poi tutto come prima. Ecco lo screenshot del callgraph per il fortran:

Qui notiamo che la funzione per la stampa a video (_gfortran_transfer_real_write) prende ben il 13% del tempo, ed è probabilmente lei la causa delle peggiori prestazioni del fortran rispetto al C/C++. Inoltre, notate come da nessuma parte ci sia una chiamata ad alcuna pow(), perché nel fortran l’elevazione a potenza per un numero intero è una operazione del linguaggio.

Ovviamente, questi tool funzionano bene per codice eseguibile, non funzionano affatto bene su Java e su Python. Ci sono però altri tool, più o meno professionali, che fanno cose analoghe. Semmai se ne parla un’altra volta!

Spero che questo post possa esservi utile nel vostro lavoro di “ottimizzatori!”, e ricordate le parole di Knuth!

Nota di Juhan: Mi ricordo il profiler id Borland; questo è molto più bello e anche il prezzo. Vero? Vero dikiyvolk?

Fortran – domande e risposte

Queste sono alcune domande che mi sono stare rivolte riguardo ai post sul Fortran. Di alcune ho la risposta, di altre the answer is blowin’ in the wind, (Dylan 1962).

Contributi benvenuti, come al solito. Non c’è un ordine particolare, non me le ricordo tutte (non ho preso appunti, azt!).

D: Come si fa a passare da stringa a numero e viceversa? L’equivalente di questo codice Python.

R: In fortran esiste, da F77, un modo elegante e atipico che consiste nel leggere dalla variabile stringa il valore numerico o scrivere quest’ultimo in una stringa. Il codice equivalente a quello Python è il seguente:

program intRW
implicit none

integer :: i
real :: r
character(16) st, rst

st = '1234'
read(st, *) i
print *, st, i

r = 42.8
write(rst, *) r
print *, r, rst

end


Anticamente esisteva una coppia di funzioni encode e decode che i vecchi ancora rimpiangono, ma niente da fare: rimosse.

D: A proposito scrivere in modalità ‘*‘ (default) è spesso brutto; quali sono le specifiche per format e fmt?

R: Una volta quando i manuali erano di carta c’erano un paio di tabelle (probabilmente le più fotocopiate) che unite a un righello meccanografico con le dimensioni dei caratteri erano uno strumento indispensabile del programmatore. Poi i numeri si mettevano comunque dentro per successive approssimazioni. Oggi si può trovare qualcosa di simile sul Web, per esempio questo

D: È possibile mescolare codice vecchio (formattato) e nuovo (free-format)?

R: Sì, senza alcun problema. In particolare il compilatore gfortran ipotizza il free-format in base all’estensione:

Source files with ‘.f’, ‘.for’, ‘.fpp’, ‘.ftn’, ‘.F’, ‘.FOR’, ‘.FPP’, and ‘.FTN’ extensions are treated as fixed form. Source files with ‘.f90’, ‘.f95’, ‘.f03’, ‘.f08’, ‘.F90’, ‘.F95’, ‘.F03’ and ‘.F08’ extensions are treated as free form. The capitalized versions of either form are run through preprocessing. Source files with the lower case ‘.fpp’ extension are also run through preprocessing.

(da Using GNU Fortran).
Esistono poi le opzioni -ffree-form e -ffixed-form per casi particolari.
Secondo me se si hanno file vecchi scritti in fixed-form OK ma i nuovi vanno scritti in free-form: cambiano solo le convenzioni per il commento e la continuazione su più righe, dai!

D: Conviene suddividere il programma in più file?

R: Quelli grossi sì. Tra le altre cose è più agevole la compilazione, specie quando si è alla caccia di errori. Anzi non sarebbe una cattiva idea usare make.
Probabilmente conviene avere un file per il main e uno per ogni modulo, oltre a quelli per gli include.

D: I moduli possono essere visti come librerie?

R: I moduli sono fatti per contenere parametri e variabili globali e procedure (funzioni e subroutine) per cui sì, possono essere visti come librerie. Ma sono una cosa un po’ diversa da quelle presenti negli usuali linguaggi; a me ricordano gli unit del Borland Pascal poi diventato Delphi. OK, c’è using del Visual Basic e chissà quanti altri ma io parlavo di cose serie ;-)

D: Negli esempi usi spesso e funzioni e poco le subroutine; c’è un motivo?

R: No, e sono al corrente che di solito i fortranisti preferiscono le subroutine. Però quando ci si aspetta un solo valore di ritorno la funzione è più comoda, secondo me. Poi possono esserci effetti collaterali (side effects) ma questo è un altro discorso e probabilmente un male. Giova ricordare che in Fortran tutti i dati vengono passati per riferimento, a differenza della maggior parte degli altri linguaggi (C e Python p.es.). Per evitare errori comuni nelle ultime versioni del Fortran è stato introdotto intent che può essere in, out o inout.
Ci sono poi linguaggi in cui mancano le subroutine, hanno solo le funzioni, al limite possono non ritornare nessun dato (void) o puoi ignorare il valore ritornato, come il C e il Lisp.

D: Puoi elencare i link ai documenti che hai usato?

R: Certo, eccoli:

il manuale di gfortran: Using GNU Fortran;
TACC: Modern Programming Languages: Fortran90/95/2003/2008;
Tanja van Mourik Fortran 90/95 Programming Manual;

D: Nomini spesso il Lisp: ma davvero ti piace?

R: Sì :-D

Nota: lo so che ho promesso a Walter ulteriori confronti tra linguaggi, compreso il C++ ma questo post è in cantiere da una settimana, sono in ritardo :-(

Iscriviti

Ricevi al tuo indirizzo email tutti i nuovi post del sito.

Unisciti agli altri 71 follower