Tecnicamente's Weblog

Just another WordPress.com weblog

Archivi Categorie: Fortran

Fortran95, il ciclo “DO” su matrici di grandi dimensioni

Se state usando il Fortran95 i motivi sono semplici, o vi stanno obbligando oppure volete usare un linguaggio di programmazione che non ha eguali nelle applicazioni del calcolo numerico ed in particolar modo nei calcoli dove si fa uso pesante di matrici.

Questo breve articolo ha due scopi:

  1. mostrare che bisogna organizzare le matrici nel modo migliore per permettere ai cicli di lettura di essere molto efficienti;
  2. mostrare le ottimizzazioni intrinseche del linguaggio e del compilatore;

Supponiamo di avere due matrici di grande dimensione la [A] e la [B], il calcolo che si vuole eseguire è c(i,j) = a(i,j) + b(i,j). Voglio trovare la matrice somma [C]=[A]+[B].

Operazione estremamente semplice, supponiamo ancora che decidiamo di usare il linguaggio C e che i dati siano interi, vediamo cosa accade usando il seguente codice:

#include<stdio.h>
#include<time.h>

#define X 10000
#define Y 10000

int main(void){
 int i, j;
 float t0, t1;
 //int a[X][Y], b[X][Y], c[X][Y];
 int **a, **b, **c;

/* allocazione della memoria */
 a = (int **)calloc(X,sizeof(int *));
 b = (int **)calloc(X,sizeof(int *));
 c = (int **)calloc(X,sizeof(int *));

 for(i=0;i<X;i++){
    a[i] = (int *)calloc(Y,sizeof(int));
    b[i] = (int *)calloc(Y,sizeof(int));
    c[i] = (int *)calloc(Y,sizeof(int));
 }

/* inizializzazione delle matrici */
 for(i=0;i<X;i++){
    for(j=0;j<Y;j++){
       a[i][j] = 1;
       b[i][j] = 2;
    }
 }
 
/* somma scorrendo per le colonne */
 t0 = clock();
 for(i=0;i<X;i++){
   for(j=0;j<Y;j++){
     c[i][j] = a[i][j] + b[i][j];
   }
 }
 
 t1 = clock();
 printf("%d\n",c[X-1][Y-1]);
 printf("tempo richiesto somme su colonne = %f\n",(double)(t1 - t0) / CLOCKS_PER_SEC);

/* somma scorrendo le righe */
 t0 = clock();
 for(j=0;j<Y;j++){
   for(i=0;i<X;i++){
     c[i][j] = a[i][j] + b[i][j];
   }
 }
 t1 = clock();
 printf("%d\n",c[X-1][Y-1]);
 printf("tempo richiesto somme su colonne = %f\n",(double)(t1 - t0) / CLOCKS_PER_SEC);

 return 0;
}

Il cui risultato sul mio povero PC è rappresentato nella seguente immagine:

C

Ok ma perchè ho messo un codice in C ed un test del risultato?

Ovviamente per mostrare come in Fortran95 sia comodo descrivere lo stesso algoritmo e per vedere quali sono i risultati in termini di tempo di accesso, quindi qui il codice:

program prova
implicit none
  integer :: i,j
  integer, parameter :: x=10000, y=10000
  real :: t0,t1
  integer, dimension(x,y) :: a,b,c

  a=1 ; b = 2

!SOMMA SCORRENDO LE COLONNE
  call cpu_time(t0)
  do i=1,x
    do j=1,y
      c(i,j) = a(i,j) + b(i,j)
    enddo
  enddo
  call cpu_time(t1)
  print*, c(x,y)
  print*, "tempo richiesto somme su colonne = ",t1-t0
  print*,""

!SOMMA SCORRENDO LE RIGHE
  do j=1,y
    do i=1,x
      c(i,j) = a(i,j) + b(i,j)
    enddo
  enddo
  call cpu_time(t1)
  print*, c(x,y)
  print*, "tempo richiesto somme su righe = ", t1-t0
  print*,""

!SOMMA IN MODO NATIVO DEL FORTRAN
  call cpu_time(t0)
  c=a+b
  call cpu_time(t1)
  print*, c(x,y)
  print*, "tempo richiesto somme tra matrici = ", t1-t0

end program prova

Con il seguente risultato:F

Come vedete la principale differenza tra i due linguaggi è che in Fortran le MATRICI ed i VETTORI sono tipi nativi, quindi è possibile fare delle operazioni direttamente su di esse.

Occorre tenere presente che ho usato i seguenti compilatori gcc e gfortran, quindi il backend è lo stesso, il perchè della differenza di risultato lo lascio agli esperti, io non so spiegarlo.

E’ interessante notare come una implementazione non oculata degli accessi agli array in C può essere disastrosa.

Tutti i modi usati in Fortran sono migliori, in particolare l’uso delle operazioni a livello “nativo” sulle matrici permette di ridurre notevolmente le righe di codice oltre al non trascurabile fatto avere tempi di calcolo migliori che in C.

Una interessante possibilità dei compilatori fortran è quella di parallelizzare il codice, cosa che purtroppo non riesco a testare :(.

 

Buon Coding in Fortran.

Annunci

Come installare Intel(R) Fortran Compiler XE 12.1 su Ubuntu 12.04 LTS

Intel offre da anni una versione gratis del suo compilatore fortran per scopi non commerciali e non c’è dubbio che tale compilatore sia tra i migliori se non il migliore attualmente disponibile.

Quindi se state realizzando un software non commerciale, una bella compilazione con il compilatore Intel e le sue ottimizzazioni per i processori Intel potrebbe darvi grandi soddisfazioni.

Bene cominciamo, per prima cosa occorre il compilatore, lo trovate qui: http://software.intel.com/en-us/articles/non-commercial-software-download/

Scaticate e scompattate il file e lanciate come root il file “install.sh”, vi apparirà il seguente menù:

date invio e andate avanti.

Ovviamente questo tutorial non avrebbe senso se l’installazione fosse stata semplice :), infatti vi trovere subito dei problemi:

Selezionate 2 per analizzare i dettagli.

Ed ecco l’origine dei nostri mali:

Provate a forzare l’installazione selezionando 1.

Legge la licenza ed alla fine digitate “accept”.

Alla schermata per l’attivazione scegliete il menù 1:

A questo punto occore inserire il codice seriale che vi è stato spedito via email durante la registrazione per il download del compilatore.

Alla seguente schermata indicate dove installare i files, /opt va benissino:

Bene siamo quasi alla fine seguite le indicazione dell’ultima schermata:

date quindi un bel “$ sudo /opt/intel/bin/compilervars.sh ia32”

Nel mio caso tale comando non ha dato esito, ho quindi provveduto per via manuale alla modifica della variabile $PATH aggiungendo il percorso dei file binari del compilatore.

date i seguenti comandi:

$ cd
$ gedit .profile

e modificate il file come segue:

tutto fatto 🙂

Buon lavoro

Fortran parte 3: Fortran Vs. C/C++ o compilatore Vs. compilatore?

Basta fare una ricerca dove si vogliono recuperare informazioni circa le performances offerte dal Fortran95 contro C/C++ e si trovano centianai di benchmark a favore dell’uno e dell’altro linguaggio.

Se poi si va su gruppi di discussione si assiste a vere e proprie guerre di religione per partito preso e che quasi sempre si concludono con insulti.

Non è mia intenzione schierarmi dall’una o dall’altra parte, credo semplicemente che i software, come qualsiasi altra cosa, sono dei prodotti che devono essere ingegnerizzati.

Quindi quando c’è da metter giù due righe di codice si può partire alla “garibaldina”, ma quando le cose sono serie ed impegnative occorre spendere del tempo per svolgere delle analisi e successivamente impostare l’architettura del software e la scelta degli strumenti da usare.

Tra gli strumenti c’è sicuramente il linguaggio di programmazione, difficilmente un solo linguaggio è adatto a coprire tutti gli scopi.

Attenzione, con un linguaggio come ad esempio il C o il C++ si può fare praticamente tutto, non a caso sono “general purpose“, e devo essere sincero, sono un grande ammiratore di questi linguaggi e delle loro potenzialità.

Ma in un ambiente di produzione il fattore tempo è fondamentale, il ricorso a strumenti che facilitano certe parti del progetto è, a mio parere, obbligatorio.

Il Fortran non si pone come un linguaggio per risolvere qualsiasi problema, il Fortran si pone come un linguaggio per la soluzione di problemi numerici.

In tale ottica offre una serie di strumenti che altri linguaggi non hanno, per capirci dovendo spiegare quale è il linguaggio più somigliante al fortran, io dire Matlab.

La sostanziale differenza è che Fortran necessita di compilatori, Matlab mi pare sia in pseudocodice.

Ma torniamo alla questione Fortran Vs. C/C++, voglio solo fare degli esempi per dare una idea di quello che mi piace del Fortran ma anche del C/C++.

Si prenda in considerazione l’esempio seguente, una matrice B contiene dei valori numerici, ed una seconda matrice A è nulla.

Si vuole che A assuma gli stessi elementi di B quando questi sono pari.

Il problema proposto è solo uno spunto per mostrare le capacità del linguaggio.

In C si può scrivere il seguente programma:

#define X 15000
#define Y 15000
#define ESECUZIONE 5

#include<stdio.h>
#include<stdlib.h>
#include<malloc.h>

int main(void)
{
  int i,j,k;
  int **a;
  int **b;

  for(k=0;k<ESECUZIONE;k++)
  {
    a = malloc(sizeof(int)*X);
    b = malloc(sizeof(int)*X);

    for(i=0;i<Y;i++)
    {
      a[i] = malloc(sizeof(int)*Y);
      b[i] = malloc(sizeof(int)*Y);

      for(j=0;j<Y;j++)
      {
        b[i][j] = j+1;
      }
    }

    for(i=0;i<X;i++)
    {
      for(j=0;j<Y;j++)
      {
        if((b[i][j]%2)==0)
        {
          a[i][j] = b[i][j];
        }
      }
    }
    printf("%d\n",a[10-1][10-1]);

    for(i=0;i<X;i++)
    {
      free(a[i]);
      free(b[i]);
    }
    free(a);
    free(b);
    return 0;
}

return 0;
}

Come si potrebbe scrivere lo stesso programma in Fortran95?
Ecco la risposta:

program esempio
implicit none

integer, parameter :: X = 15000
integer, parameter :: Y = 15000
integer, parameter :: ESECUZIONE = 5

integer :: i,j,k

integer, allocatable :: A(:,:), B(:,:)

do k=1,ESECUZIONE
 allocate(A(X,Y))
 allocate(B(X,Y))

 A = 0

 do j=1,Y
     B(:,j) = (/(i,i=1,X)/)
 enddo

 where(mod(B,2)==0) A=B

 write(*,*) B(10,10)

 deallocate (A)
 deallocate (B)

enddo
end program esempio

o più sinteticamente:

program esempio_2
implicit none

integer, parameter :: X = 15000
integer, parameter :: Y = 15000
integer, parameter :: ESECUZIONE = 5

integer :: i,j,k

integer, allocatable :: A(:,:), B(:,:)

do k=1,ESECUZIONE
	allocate(A(X,Y))
	allocate(B(X,Y))

	A = 0

	forall(j=1:Y)
		B(:,j)=j
	end forall

	where(mod(B,2)==0) A=B

	write(*,*) B(10,10)

        deallocate (A)
	deallocate (B)
enddo
end program esempio_2

La cosa evidente, a prescindere da fatto che abbiate capito qualcosa nel codice Fortran :), è che in Fortran le linee di codice richieste sono effettivamente poche, e vi assicuro che se si conosce il linguaggio sono anche estremamente leggibili.

Vi ho mostrato due versioni in Fortran95, la differenza tra i due sta nel fatto che nel primo codice uso un ciclo DO…ENDDO per riempire la matrice B lungo le colonne, ciò è necessario per velocizzare gli accessi in memoria visto che il Fortran organizza i dati in memoria seguendo l’ordine delle collone, a differenza del C/C++ che seguno le righe.

Nel secondo codice il ciclo DO…ENDDO precedende è sostituito dal un ciclo FORALL…END FORALL, tale ciclo è molto particolare, e dovrebbe dare il meglio di se su architetture multiprocessore.

Infatti FORALL implementa il calcolo parallelo nativamente e direttamente dal compilatore, cosa non da poco, purtroppo io dispongo di un sistema non proprio performante e per di più monoprocessore infatti eccovi la scheda del mio sistema:

Adesso passiamo alla copilazione ed al “benchmark” dei codici.

Per la precisione userò una versione per usi non commerciali del compilatore Fortran di Intel e l’immancabile gcc per il codice in C.

Compilo i primi due codici ed ottengo:

Quale è la conclusione?

La velocità dipende:

– dal tempo impiegato a scrivere tutte le righe di codice

– da come il programmatore implementa il codice

– dal tempo impiegato a debuggare

– dalle ottimizzazioni realizzate dal compilatore

Circa l’ultimo punto ho una domanda, come mai gfortran che poi è gcc impiega più tempo rispetto alla versione in C, mentre il compilatoreo ifort impiega praticamente lo stesso tempo di gcc?

Ovviamente questo test non è un test valido per qualificare le prestazioni dei compilatori.

Se qualcuno prova ed eseguire i test su una macchina dual core mi farebbe piacere avere i risultati.

In tal caso il codice fortran da usare è il secondo dove è presente l’istruzione FORALL.

Appena recupero i sorgenti della mia tesi, vedo di eseguire un test sui diversi compilatori, a suo tempo impiegava 60 ore, compilatore Compaq Visual Fortran e un bel PC da 3200 Mhz.

Fortran parte 2: Passare i parametri dalla riga di comando

Mi è sempre piacuta la possibiltà di passare dei parametri ad un programma nel momento in cui lo si lancia dalla riga di comando.

Ciò è estremamente semplice in C/C++, basta infatti:

#include<stdio.h>
int main(int argc, char* argv[])
{
    int i;
    printf("Passati %d parametri\n",argc);
    printf("I parametri sono:\n");
    for(i=1;i<argc;i++) { printf("%s\t",argv[i]);}
    printf("\n");  
    return 0;
}

Sono rimasto sorpreso e stupito che lo si possa fare con il fortran9x, occore richiamare due funzioni native ed ecco:

program parametri
 implicit none
 integer :: argc, i
 character(len=256) :: argv

 argc = command_argument_count()
 argnum = command_argument_count()
 write(*,*)  "Numero argomenti: ",argc 
 do i = 1, argc
  call get_command_argument(i,argv)
  write(*,*) argv
 end do
end program parametri


Un modo diverso per ottenere lo stesso risultato  è il seguente:

program parametri
    implicit none
    
    integer :: i, n_of_arg
    character(10) :: argv
    
    
    write(*,*) "Numero argomenti: ", iargc()

    do i=1,iargc()
        call getarg(i,argv)
        write(*,*) argv
    enddo
end program parametri

Fortran Parte I: Gli array in Fortran95

Molti ridono quando sentono parlare del Fortran, ma di fatto chi programma in tale linguaggio non è mai un informatico di professione.

In genere chi usa tale linguaggio è una persona che necessita di scrivere del software che si deve interfacciare direttamente a persone o sistemi, ma per la precisione, scrive  software per il calcolo numerio.

Tutti i linguaggi sono pensati per eseguire calcoli, ma nessuno permette di ottenere le prestazioni e la semplicità di base che si può ottenere con il Fortran95.

In questa breve nota mostrerò l’aspetto che più apprezzo di questo linguaggio, e che alla fine è quello che mi impedisce anche di passare al C o similari, parliamo di ARRAY!

ARRAY in FORTRAN

Un Array è un insieme di dati omogenei, e non differesiche in nulla da quelli che sono gli array in altri linguaggi.

Possiamo quindi avere:

  • Array di numeri interi
  • Array di numeri reali
  • Array di numeri complessi
  • Array di caratteri
  • Array di tipi personali, cioè di tipi creati per uno scopo preciso

Ancora una volta gli array possono essere statici o dinamici, e vedremo più avanti cosa vuol dire.

Per un ingegnere, un fisico, un matematico, un chimico, una persona che vuol fare dei calcoli, l’array non è altro che una MATRICE.

Possiamo avere array monodimensionali, quindi avremo un VETTORE, o array multimensionali quindi MATRICI.

In molte materie scientifiche si realizza un modello matematico basato su matrici ed il calcolo matriciale per semplificare il problema.

Gli esempi sono innumerevoli:

  • L’analisi del moto dei corpi celesti
  • La scienza delle costruzioni
  • L’analisi agli elementi finiti
  • La fluido dinamica numerica
  • L’aerodinamica applicata
  • ecc.

Per cui possedere un linguaggio che gestisca in modo diretto le matrici non è cosa da poco.

Passiamo agli esempi 🙂

Dichiarazione degli ARRAY

integer :: matrice(10,10)
integer, dimension(10,10) :: matrice

con le due dichiarazioni precedenti si crea una matrice di interi di rango 10×10

analogamente si può creare un vettore di numeri reali o complessi:

double precision :: vettore(10)
double precision, dimension(10) :: vettore

Le dichiarazioni precedenti sono usate quando si conosce a priori la dimensione dell’array da usare, in molti casi si può avere a che fare con array dalle dimensioni sconosciute, in tal caso occore usare gli array dinamici ed allocare di volta in volta lo spazio i memori necessario.

La dichiarazione è simile alla precedende, ma al posto della dimensione si inserisce il simbolo “:”.

integer, allocatable :: matrice(:,:)
integer, allocatable, dimension(:) :: vettore
ALLOCATE(matrice(10,10))
ALLOCATE(vettore(10))

Per liberare la memori quando le matrici non sono più necessarie basta:

DEALLOCATE(matrice)
DEALLOCATE(vettore)

ASSEGNAZIONI DI VALORI AGLI ARRAY

Gli array possono essere usati come semplici variabi sia con la metodologia di accesso in lettura e scrittura basata su indici, quindi:

WRITE(*,*) A(i,j) -> stampa l'elemento di coordinate (i,j)
A(i,j) = 5 -> assegnail valore 5 all'elemento di coordinate (i,j)

Ma fin qui tutto è in linea con le possibilità offerte da C/C++, vediamo quindi le peculiarità circa l’uso degli array.

INTEGER :: A(10,20;30)  ! creo un array di interi 10x20x30
A = 1 ! tutti gli elementi dell'array sono pari ad 1
A(1,:,:) = 2 ! tutti gli elementi che hanno il primo indicie pari a 1 assumono valore 2, quindi A(1,1,1) = 2, A(1,1,3) = 2, ..., A(1,20,30) = 2
A(1,1,:) = 3 ! tutti gli elementi di coordinate inizilia (1,1) hanno valore = 3, quindi A(1,1,1)=3, A(1,1,2)=3, ...; A(1,1,30) = 3
C = A+B !somma indice ad indice gli elementi di un array ovviamente si può avere anche la sottrazione
 
C=A*B ! moltiplica indice ad indice 

Tali modi di accesso possono essere usati per ridurre le righe di codice, semplificare la lettura e permette al compilatore di implementare l’ottimizzazione migliore di accesso, specialmente quando si ha a che fare con migliaia o milioni di elementi.

FUNZIONI NATIVE PER GLI ARRAY

Il linguaggio offre delle funzioni native per particolari operazioni su array quali:

DOT_PRODUCT(A, B) ! prodotto scalare di A e della B
MAXVAL(A)  !  restituisce il valore maggiore dell'array A
MINVAL(A)  !restituisce  il valore minore di A
MAXLOC(A)  !  restituisce le coordinate le valore maggiore di A
MINLOC(A)  ! restituisce le coordinate del valore minore di A
PRODUCT(A)  ! prodotto degli elementi di A
SUM(A) ! somma gli elementi di A
TRANSPOSE(A) ! crea la matrice trasposta

e molte altre ancora per fare lo scorrimento degli elementi, la scelta di elementi da leggere o scrive con l’applicazioni di filtri (chiamati maschera) ecc.

Il mio intento non è quello di scrivere una guida, se il linguaggio vi interessa potete chiedere qui o consultare i diversi manuali che si trovano in rete.

vi lascio con un piccolo esempio:


Qui il risultato:

Ci sarebbe ancora molto da scrivere, ma sarà per la prossima volta.

Infatti non abbiamo parlato di compilatori, lo faremo prossimamente, per ora vi lascio qualche link per curiosare e provare se queste poche note vi hanno incuriosito.

http://hpff.rice.edu/index.htm

http://www.dmi.units.it/~chiarutt/didattica/parallela/fortran/index.html

http://www.nsc.liu.se/~boein/f77to90/f77to90.html#index

Numeri casuali in Fortran 90, 95 e 2003

Tempo fa, quando su linux, freeBSD ed altri sistemi simili non esisteva un valido compilatore fortran 90 ho usato il compilatore fornito dall intel trovandomi benissimo, ed abituandomi inevitabilmente alla sua implementazione dello standard del linguaggio.

Dopo anni  ho rimesso mano al linguaggio, e non nego con qualche ruggine mentale :), niente di piu’ bello che installare gfortran, scrivere due righe di codice e constatare i primi problemi 🙂 infatti, dopo aver compilato il seguente codice:

program prova
implicit none
integer :: i
real     :: r
call random_seed()
do i=1,10
    call random_number(r)
    write(*,*) r
enddo
end program prova

Mi trovo che ad ogni esecuzione i numeri “pseudo-casuali” generati sono sempre gli stessi.

Installo anche g95 ricompilo e qui le cose cambiano, i numeri generati sono differenti ad ogni esecuzione del codie!

Ma allora che succede?

Penso ad un bug a qualcosa di anomalo, e mi informo… risultato?

Lo stardard è ambiguo circa la generazione di numeri casuali, per cui il comportamento di gfortran è corretto quanto quello di g95, ifort ed altri.

Bene, occorre rimediare, comincio a leggere un po’ di documentazione ed ecco questa piccola subroutine, che consiglio di usare come liberia all’interno di un “module”.

Usandola si otterrà del codice che genera differenti numeri casuali in ogni esecuzione del programma, che quello che ci sia aspetta 🙂

Eccovi il codice e buon lavoro, attenzione non funziona con compilatori fortran77 come il g77 :

SUBROUTINE init_random_seed()
 INTEGER :: i, n, clock
 INTEGER, DIMENSION(:), ALLOCATABLE :: seed

 CALL RANDOM_SEED(size = n)
 ALLOCATE(seed(n))

 CALL SYSTEM_CLOCK(COUNT=clock)

 seed = clock + 37 * (/ (i - 1, i = 1, n) /)
 CALL RANDOM_SEED(PUT = seed)

 DEALLOCATE(seed)
END SUBROUTINE

Programmare in FORTRAN?

Il fortran come molti sanno è uno dei primi linguaggi di programmazione nati per facilitare la scrittura dei programmi.

In passato ha avuto il momento di successo, specialmente in ambito schientifico-ingegneristico.

Poi a poco a poco il suo utilizzo e la sua fama sono decaduti a favore di linguaggi come C e C++.

Ma veramente il FORTRAN è un linguaggio vetusto?

Personalmente penso di no.

Per la tesi di laurea e per passione ho programmato usando diversi linguaggi, tra cui C.

E’ ciò mi permette di dire con assoluta neutralità che il FORTRAN è un linguaggio comodo e potente per le applicazioni di calcolo puramente numerico.

Sia chiaro non parlo delle versioni fortran66 o fortran77, ma delle più recenti Fortran90 e 95.

Infatti con tali nuovi standard il linguaggio è statto ammodernato e svecchiato.

E’ stata introdotta la gestione nativa delle matrici e di svariate funzioni di gestione delle stesse ed altro.

Se siete avvezzi a usare matlab o scilab allora nel nuovo fortran95 troverete un amico potente, facile da usare e veloce.

Altra cosa importante e che nel mio caso ha fatto la differenza nello scrivere un software di calcolo complicato è stata l’attenzione che il fortran95 da alla limitazione di errori per sforamenti della memori, usi errati degli indici degli array filtro sul tipo di dati passati alle funzioni.

Insomma dovete provare per capire la bellezza che questo linguaggio ancora porta con se.

quindi?

quindi eccovi un piccolo link dove cominciare l’avventura: http://www.thefreecountry.com/compilers/fortran.shtml