Sul FORTRAN di EWD

A 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.

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

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:

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

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.

Finora 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

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

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

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ì

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

Fortran – Ultima puntata

Pointers
OK, adesso anche nel Fortran ci sono i puntatori (pointer). Siccome questi post sono legati a un programma specifico non ne parlerò adesso, forse in futuro.

Preprocessore
Come per il C è disponibile il preprocessore, quindi #ifdef, #ifndef &co. In realtà si usa GCC, sempre lui.

Ricorsività
Adesso una procedura può chiamare sé stessa ricorsivamente, cosa vietata fino al 77 compreso. Ecco per esempio il calcolo del fattoriale.

program nfact ! da Tanja van Mourik
implicit none

integer :: n
integer :: nfactorial

do n = 2, 9, 3
    print *, n, nfactorial(n)
end do

end program nfact

recursive function nfactorial(n) result(fac)
implicit none

integer :: fac
integer, intent(in) :: n

select case(n)
    case (0:1)
        fac = 1
    case default
        fac = n * nfactorial(n-1)
end select

end function nfactorial

Notare il select case e che la funzione dev’essere dichiarata recursive.

Operazioni sugli array
Gli array possono essere inizializzati sia con (/ e /) che con [ e ].
Possono essere copiati con una singola istruzione, evitando un do specifico.

program arr
implicit none

integer :: a(5) = (/1, 2, 3, 4, 5/)
integer :: b(5) = [1, 4, 9, 16, 25]

print *, 'inizio'
print "(5i3)", a, b

a = b ! copia dell'array
print "(/, 'dopo la copia')"
print "(5i3)", a

a(2:5) = a(1:4)
a(1) = 0
print  "(/, 'dopo lo spostamento')"
print "(5i3)", a

end

Il tentativo di invertire ‘al volo’ un array da errore, la copia funziona solo per slice (chissà come si dice da noi?).

program arr
implicit none

integer :: a(5) = (/1, 2, 3, 4, 5/)
integer :: b(5) = [1, 4, 9, 16, 25]

a(1:5) = b(5:1)
print "(5i3)", a, b

end

In questo caso si deve ricorrere al solito ciclo.

program arr
implicit none

integer :: a(5)
integer :: b(5) = [1, 4, 9, 16, 25]
integer :: i

do i = 1, 5
    a(i) = b(6 - i)
end do
print "(5i3)", a, b

end

Array dinamici
Copiando da Tanja: alle volte la dimensione di un array non è nota fino all’esecuzione del programma, ecco allora gli array dinamici.

program dynarr
implicit none

real, allocatable :: a(:,:)
integer :: dim1, dim2
integer :: i, j

write(*, fmt = '("Give dimensions DIM1 and DIM2: ")', &
                advance = "NO")
read(*, *) dim1, dim2

! now that the size of a is know, allocate memory for it
allocate(a(dim1,dim2))
do i = 1, dim2
    do j = 1, dim1
        a(j,i) = i*j
        write(*, fmt = '("a(",i2,",",i2,") = ", f6.2)') &
              j, i, a(j, i)
    end do
end do
deallocate(a)

end

Domanda di G. (ripetuta, ecco la risposta): si può continuare a usare quello che si usava una volta? Certo, ecco:

program oldstyle
implicit none
common /mydat/ numero

integer :: i, numero, molt

i = 5
numero = 3

print *, i, molt(i)

end

integer function molt(n)
implicit none
common /mydat/ mul

integer n, d, mul

d = mul * n

molt = d
return !non serve ma si usava

end

Il common è comunque sconsigliato, usare i moduli. Con questi posso importare selettivamente le variabili che mi servono effettivamente, anche con nomi differenti.

program es_only
use m_only, only: alpha,  b => beta
! b è beta del modulo

implicit none

integer :: gamm, delta
! gamm e delta NON sono quelli del modulo

print *, alpha, b, gamm, delta

end

module m_only
implicit none

integer :: alpha = 1, beta = 2, gamm = 3
integer, private :: delta = 4

end module

OK, con questo credo di aver elencato (quasi) tutti le nuove caratteristiche del linguaggio, rispetto al 77. In realtà ho saltato quelle che non credo saranno utilizzate in un caso reale che partirà a breve. Probabilmente posterò ancora un elenco di domande tipiche, tipo questa:
D: con il formato libero puoi mettere i commenti con una C in prima colonna?
R: prova ma no, ecco:

program c_com

C questi una volta 
c sarebbero stati commenti

i = 42

print *, i

end

Beh, ce ne sono anche di meno banali

Ancora confronti sui linguaggi

Intanto questa dovrebbe essere l’ultima puntata del tormentone sviluppato negli ultimi post: comincio a detestarlo al pari del caldo afoso che perdura.

Poi devo confessare un mio errore, grave: nella versione Fortran del listato ho postato un file non aggiornato, quello del modulo m_conf.f08. Questa è la versione corretta:

module m_conf
implicit none

real(8) :: f(13)

contains
    real(8) function fsin(a) result(r)
        real(8) :: a
        r = a - a ** 3 / f(3) + a ** 5 / f(5) &
            & - a ** 7 / f(7) + a ** 9 / f(9) &
            & - a ** 11 / f(11) + a ** 13 / f(13)
    end function fsin

    real(8) function fcos(a) result(r)
        real(8) :: a
        r = 1 - a ** 2 / f(2) + a ** 4 / f(4) &
            & - a ** 6 / f(6) + a ** 8 / f(8) &
            & - a ** 10 / f(10) + a ** 12 / f(12)
    end function fcos

end module

Conseguenza della mancata dichiarazione di real(8) per le funzioni fsin() e fcos() queste restituivano un valore real(4). Ma attenzione: siccome tutte le variabili erano in doppia precisione (ahemm, real(8), una volta si diceva così) il tempo di esecuzione non cambia!

Allora l’errore –OK, il grave errore– è quasi solo metodologico. E a mia giustificazione potrei citare il fatto di essere costantemente interrotto da infinite incombenze. E dal caldo, non dimentichiamo il caldo. C’è poi un altro motivo, più grave: a lavorare da solo non hai la possibilità di farlo vedere a qualcuno, raccontarglielo e allora queste cose salterebbero fuori subitissimo.

In ogni caso per verificare se il Fortran è davvero così più performante rispetto a Python ho modificato i codici facendo ripetere i calcoli 10 volte e introducendo il calcolo dei logaritmi decimali di seno e coseno.

Io da piccolo invece del liceo ho frequentato l’istituto per geometri (sono figlio di contadini). Uh! anche Piergiorgio Odifreddi e Flavio Briatore hanno fatto i geometri, pensa te!
A quei tempi non esistevano le calcolatrici elettroniche e i conti si dovevano fare a mano. In particolare per il corso di topografia questi assumevano un’importanza preponderante e si usavano i logaritmi, quelli del Bruhns. Chi fosse Karl Christian Bruhns l’ho scoperto solo l’altro giorno, grazie al solito Wiki. Devo confessarvi che adesso il Bruhns mi risulta simpatico; ma solo da adesso

OK, fine dell’intermezzo nostalgico; il calcolo viene fatto usando la serie

presa da qui.

Poi la trasformazione da logaritmo naturale a decimale richiede solo una divisione ma nello spirito del test si è resa come una funzione. Questo perché anticamente, quando i ‘puter erano grossi e lenti, sarebbe stato un errore grave: chiamare una funzione richiede(va) aggiungere istruzioni, in questo caso inutili. Ma basta divagazioni, ecco le versioni aggiornate.

Python

#!/usr/bin/env python

import math

v = 2.0
f = ([1.0, 1.0, v])
for i in range(3, 14):
    v *= i
    f.append(v)

def fsin(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]

def fcos(a):
    return 1 - 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]

def ln(x):
    if x == 0.0:
        return -1.0e100
    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)

def log(x):
    return x / ln10

pi = math.pi
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
ln10 = ln(10.0)
fmt = "{:3d} {:10.8f} {: 10.8f} {: 10.8f} {: 10.8e}"
fml = "{:12s} {: 12.8e} {: 12.8e}"

ssum = 0.0
csum = 0.0
for i in range (0, nsteps + 1):
    a = i * step
    s = fsin(a)
    ls = log(ln(s))
    ssum += s
    c = fcos(a)
    lc = log(ln(c))
    csum += c
    if (i % (10 * deg)) == 0:
        if c != 0.0:
            t = s / c
        print fmt.format(int(i / deg), a, s, c, t)
        print fml.format('', ls, lc)
print ssum, csum

Fortran

program fc10l
use m_c10l
implicit none

real(8) :: v, pi, step, ssum, csum
real(8) :: a, s, c, t, ls, lc
integer(4) :: deg, nsteps, i, j
character(30) fmt_, fml_

f(1) = 1.0
do i = 2, 13
    f(i) = i * f(i-1)
end do

pi = 4 * atan(1.0)
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
ln10 = myln(dble(10.0))
fmt_ = "(i3, 3(1x, f11.8), 1x, e12.4)"
fml_ = "(12x, 2(e12.4))"

do j = 1, 10
    ssum = 0.0
    csum = 0.0

    do i = 0, nsteps
        a = i * step
        s = fsin(a)
        ls = mylog(myln(s))
        ssum = ssum + s
        c = fcos(a)
        lc = mylog(myln(c))
        csum = csum + c
        if (mod(i,10 * deg) == 0) then
            if (c /= 0.0) t = s / c
            print fmt_, int(i / deg), a, s, c, t
            print fml_, ls, lc
        end if
    end do

    print *, ssum, csum
end do
end program

 

module m_c10l
implicit none

real(8) :: f(13)
real(8) ln10

contains
    real(8) function fsin(a) result(r)
        real(8) :: a
        r = a - a ** 3 / f(3) + a ** 5 / f(5) &
            & - a ** 7 / f(7) + a ** 9 / f(9) &
            & - a ** 11 / f(11) + a ** 13 / f(13)
    end function fsin

    real(8) function fcos(a) result(r)
        real(8) :: a
        r = 1 - a ** 2 / f(2) + a ** 4 / f(4) &
            & - a ** 6 / f(6) + a ** 8 / f(8) &
            & - a ** 10 / f(10) + a ** 12 / f(12)
    end function fcos

    real(8) function myln(x) result(res)
        real(8) :: x, r
        if(x == 0.0) then
            res = -1.0e20
        else
            r = (x - 1.0) / (x + 1.0)
            res = 2.0 *(r ** 3 / 3.0 + r ** 5 / 5.0 &
                & + r ** 7 / 7.0 + r ** 9 / 9.0)
        end if
    end function

    real(8) function mylog(x) result(res)
        real(8) :: x
        res = x / ln10
    end function

end module

Java

import static java.lang.Math.pow;
import static java.lang.System.out;

class jc10l {
    static double f[];
    //static double ln10;
    static double ln10 = myln(10.0);

    public static double fsin(double a) {
        return a - pow(a, 3.0) / f[3] + pow(a, 5.0) / f[5]
                 - pow(a, 7.0) / f[7] + pow(a, 9.0) / f[9]
                 - pow(a, 11.0) / f[11] + pow(a, 13.0) / f[13];
    }

    public static double fcos(double a) {
        return 1.0 - pow(a, 2.0) / f[2] + pow(a, 4.0) / f[4]
                   - pow(a, 6.0) / f[6] + pow(a, 8.0) / f[8]
                   - pow(a, 10.0) / f[10] + pow(a, 12.0) / f[12];
    }

    public 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.0) / 3.0
                            + pow(r, 5.0) / 5.0
                            + pow(r, 7.0) / 7.0
                            + pow(r, 9.0) / 9.0);
        }
    }

    public static double mylog(double x) {
        return x / ln10;
    }

    public static void main(String args[]) {
        int i, j;
        f = new double[14];
        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 = Math.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                 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;
                    }
                    out.printf("%3d %11.8f %11.8f %11.8f %15.8e\n",
                          (i / deg), a, s, c, t);
                    out.printf("              %15.8e %15.8e\n",
                               ls, lc);
                }
            }
            out.println(ssum);
            out.println(csum);
        }
    }

}

newLISP

#!/usr/bin/newlisp

(set-locale "C") ;;; sono abituato così, perso un'ora!

;         0 1 2  3  4 5 6  7  8 9 10  11  12 13
(set 'p '(1 1 -1 -1 1 1 -1 -1 1 1 -1  -1  1  1))
(set 'f '(1))
(for (i 1 13)
    (set 't (* (p i) i (abs (f (- i 1)))))
    (push t f i)
)

(define (fsin a)
    (add a (div (pow a 3) (f 3)) (div (pow a 5) (f 5))
           (div (pow a 7) (f 7)) (div (pow a 9) (f 9))
           (div (pow a 11) (f 11)) (div (pow a 13) (f 13))
    )
)

(define (fcos a)
    (add 1 (div (pow a 2) (f 2)) (div (pow a 4) (f 4))
           (div (pow a 6) (f 6)) (div (pow a 8) (f 8))
           (div (pow a 10) (f 10)) (div (pow a 12) (f 12))
    )
)

(define (myln x)
    (if (= x 0.0)
        (set 'r -1.0e100)
        (begin
            (set 'r (div (dec x 1.0) (add x 1.0)))
            (mul 2.0 (add r (div (pow r 3) 3.0)
                            (div (pow r 5) 5.0)
                            (div (pow r 7) 7.0)
                            (div (pow r 9) 9.0)
                      )
            )
        )
    )
)

(define (mylog x)
    (div x ln10)
)

(set 'pi (mul 4 (atan 1.0))
     'deg (* 60 60)
     'nsteps (* 180 deg)
     'step (div pi nsteps)
     'ln10 (myln 10)
)

(for (j 1 10)
    (set 'ssum 0
         'csum 0
    )

    (for (i 0 nsteps)
        (set 'a (mul i step)
             's (fsin a)
             'ls (mylog (myln s))
             'c (fcos a)
             'lc (mylog (myln c))
        )
        (inc ssum s)
        (inc csum c)
        (if (!= c 0)
            (set 't (div s c)))
        (if (= (mod i (* deg 10)) 0)
            (begin
            (println (format "%3d %11.8f %11.8f %11.8f %15.8e"
                     (/ i deg) a s c t))
            (println (format "                %15.8e %15.8e"
                     ls lc))
            )
        )
    )
    (println ssum " " csum)
)

(exit)

Con newLISP ho perso un’ora a inseguire un errore mysteryuoso assay! Continuava a dirmi che le parentesi non erano bilanciate anche se l’editor mi diceva che erano OK. Dimenticare qualche parentesi con il Lisp è la cosa peggiore che un niubbo possa fare, ma io tanto niubbo non sono e poi ormai una certa pratica ce l’ho. Altro errore è quando apri e chiudi dove non devi ma neanche questo mi sembrava il caso. Poi l’illuminazione, la virgola al posto del punto. Atz! dimenticato (set-locale "C"), sì io i numeri sono abituato a scriverli alla ‘mericana.
A dire il vero ho dimenticato Nmila punti-e-virgola nella versione Java ma lì il compilatore te lo dice.

Nell’esecuzione dei confronti è tornato mooolto utile il fatto che time scrive su stderr; quindi ridirigendo l’output (stdout) sul teminale si può avere un risultato leggibile

Il risultato finale conferma sostanzialmente quello della prima versione. I rapporti rispetto al Fortran dei tempi di esecuzione diminuiscono leggermente ma rimangono elevati.

Per motivi non solo informatici pare che si userà il Fortran per i calcoli e Python per cose complementari, come la costruzione dei grafici riassuntivi delle elaborazioni. Mi sa che ne riparlerò, forse…

Un confronto Fortran – Python

OK, ancora Fortran. Ma un attimo, domanda: dovendo decidere che linguaggio usare prossimamente quale conviene?

Nel caso specifico si tratta di una piccola struttura in cui recentemente si è usato esclusivamente il Visual Basic, a parte cose vecchie in Fortran IV e 77. I candidati possibili devono essere semplici, costare poco, possibilmente multipiattaforma (gli applicativi Fortran adesso sono su Ubuntu anche se, forse, chissà, si potrebbero portare su Windows) e usati frequentemente, niente newLISP.
Tra i candidati c’è, ovviamente, Python. Oggi facciamo un confronto su un caso finto ma verosimile.

Proviamo a calcolare i valori di seno e coseno nell’intervallo [0 .. 180] gradi, per ogni secondo d’arco. Usiamo le formule dello sviluppo in serie di Taylor, limitandoci ai primi termini, quelli che consentano di avere una discreta precisione e costringano il ‘puter a fare qualche calcolo.

Questa è la versione Python

#!/usr/bin/env python

import math

v = 2.0
f = ([1.0, 1.0, v])
for i in range(3, 14):
    v *= i
    f.append(v)

def fsin(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]

def fcos(a):
    return 1 - 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]

pi = math.pi
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
fmt = "{:3d} {:10.8f} {: 10.8f} {: 10.8f} {: 10.8e}"

ssum = 0.0
csum = 0.0
for i in range (0, nsteps + 1):
    a = i * step
    s = fsin(a)
    ssum += s
    c = fcos(a)
    csum += c
    if (i % (10 * deg)) == 0:
        if c != 0.0:
            t = s / c
        print fmt.format(int(i / deg), a, s, c, t)

print ssum, csum

Notare l’array f, cosa inusuale in Python, di solito si usano le liste.
Poi per il resto normale amministrazione.

E questa è la versione Fortran. Ormai sono abituato a usare il free format, dichiarare tutte le variabili e inoltre ho scoperto che l’uso dei moduli è comodo, anzi fondamentale.

program fconf !confronto con python
use m_conf
implicit none

real(8) :: v, pi, step, ssum, csum
real(8) :: a, s, c, t
integer(4) :: deg, nsteps, i
character(30) fmt_

f(1) = 1.0
do i = 2, 13
    f(i) = i * f(i-1)
end do

pi = 4 * atan(1.0)
deg = 60 * 60
nsteps = 180 * deg
step = pi / nsteps
fmt_ = "(i3, 3(1x, f11.8), 1x, e12.4)"

ssum = 0.0
csum = 0.0

do i = 0, nsteps
    a = i * step
    s = fsin(a)
    ssum = ssum + s
    c = fcos(a)
    csum = csum + c
    if (mod(i,10 * deg) == 0) then
        if (c /= 0.0) t = s / c
        print fmt_, int(i / deg), a, s, c, t
    end if
end do

print *, ssum, csum

end program

(Guardando le righe 27 e 29 si nota una cosa che ancora manca al Fortran: gli operatori del tipo +=).

module m_conf
<pre>implicit none

real(8) :: f(13)

contains
    function fsin(a) result(r)
        real(8) :: a, r
        r = a - a ** 3 / f(3) + a ** 5 / f(5) &
            & - a ** 7 / f(7) + a ** 9 / f(9) &
            & - a ** 11 / f(11) + a ** 13 / f(13)
    end function fsin

    function fcos(a) result(r)
        real(8) :: a, r
        r = 1 - a ** 2 / f(2) + a ** 4 / f(4) &
            & - a ** 6 / f(6) + a ** 8 / f(8) &
            & - a ** 10 / f(10) + a ** 12 / f(12)
    end function fcos

end module

Piccola osservazione OT: Gedit si confonde quando visualizza file di Fortran, considera commento l’istruzione csum = 0.0 per via della C in prima colonna.

I risultati? L’output è visibile nelle immagini sopra, sostanzialmente coincidenti. I tempi di elaborazione sono interessanti, ecco la tabella di 5 run

Il Fortran è nettamente più veloce, per un fattore di 3.876 / 0.129 = 30. Però stiamo parlando di secondi. E allora questo non è più così importante e saranno da prendere in considerazione altre caratteristiche quali la semplicità e i gusti del programmatore. Anzi principalmente questi ultimi
Poi ci sarebbero, volendo, altri linguaggi

Fortran – i moduli

I moduli sono il punto d’arrivo di una linea di sviluppo antica. Probabilmente conviene esaminarne i vari componenti per comprenderne l’utilità e la bellezza.
Le procedure, funzioni e subroutine, hanno diversi scopi:

• riusare blocchi di codice;
• ripetere operazioni con set diversi di dati;
• evitare conflitti tra le variabili
• semplificare lo sviluppo (divide et impera).

Il loro posto naturale sono i moduli; inoltre questi possono contenere costanti (parameter), variabili, arrays, structures. In questo caso il modulo supera il common, in modo più flessibile. Eccone un esempio.

program e_com !esempio d'uso di common

implicit none
common /my_common/ intv, realv, complv

integer :: intv
real :: realv
complex :: complv

intv = 21
realv = 24.6
complv = (3, 4)

print *, 'inizio'
print *, intv, realv
print *,  complv

call doppio

print *, 'dopo'
print *, intv, realv
print *,  complv

end

subroutine doppio

implicit none
common /my_common/ intv, realv, complv

integer :: intv
real :: realv, re, im
complex :: complv
character :: nl =  new_line(' ')

print *, nl, 'doppio', nl

intv = 2 * intv
realv = 2 * realv
re = 2.0 * real(complv)
im = 2.0 * aimag(complv)
complv = complex(re, im)

end

Lo stresso programma riscritto usando un modulo

program e_mod !esempio d'uso di modulo

use m_mod
implicit none

intv = 21
realv = 24.6
complv = (3, 4)

print *, 'inizio'
print *, intv, realv
print *,  complv

call doppio

print *, 'dopo'
print *, intv, realv
print *,  complv

end

e

module m_mod

implicit none
save
integer :: intv
real :: realv
complex :: complv

contains
    subroutine doppio

    real :: re, im
    character :: nl =  new_line(' ')

    print *, nl, 'doppio', nl

    intv = 2 * intv
    realv = 2 * realv
    re = 2.0 * real(complv)
    im = 2.0 * aimag(complv)
    complv = complex(re, im)

    end subroutine doppio

end module m_mod

Con i moduli si possono scrivere procedure polimorfe, simili a quelle predefinite. Fin dall’antichità (quasi) si può usare, per esempio, max() al posto di amax0(), amax1(), max0() e max1(). Adesso si possono costruire! Ecco un esempio, preso da TACC, la funzione swap().

module mod_swap
public  swap
private swap_real, swap_integer

interface swap
    module procedure swap_real, swap_integer
end interface

contains
    subroutine swap_real(x, y)
    real :: x, y, t
    t = x
    x = y
    y = t
    end subroutine swap_real

    subroutine swap_integer(i, j)
    integer :: i, j, k
    k = i
    i = j
    j = k
    end subroutine swap_integer

end module mod_swap

 

program p_swap
use mod_swap

real :: a = 3., b = 4.
integer :: i = 1, j = 2

call swap(a, b)
print *, a, b

call swap(i, j)
print *, i, j

end

Posso definire operatori, ecco .plus. che fa la somma di due interi

module mod_op

public  :: operator(.plus.)
private :: f_plus
    interface operator(.plus.)
        module procedure f_plus
    end interface
contains
    function f_plus(x, y) result(res)
    integer, intent(in) :: x, y
    integer :: res
    res = x + y
    end function
end module

 

program p_op
use mod_op

i = 3
j = 5
print *, i, j, i .plus. j

end

Esistono inoltre le funzioni elemental, questa è una parola magica che estende la funzione agli array

module m_quad
contains
    elemental function quad(x) result (q)
    real, intent(in) :: x
    real :: q
    q = x * x
    end function
end module

program p_quad
use m_quad

real :: t = 12.
real, dimension(3) :: v = [1., 2., 3.]

print *, t, quad(t)
print *, v
print *, quad(v)

end

Durante la compilazione di un modulo viene creato un file .mod. Devo ancora capire quando viene utilizzato, se si rimuove l’eseguibile funziona ugualmente. Sono file di testo, che iniziano in questo modo:

GFORTRAN module version '6' created from m_mod.f08 on Mon Jun 25 09:10:02 2012
MD5:09d3ef60d0503441cdee588ae159073e -- If you edit this, you'll get what you deserve.

Prossimamente, forse