Sandro Tosi Website > Programming > C per applicazioni numeriche

Contact me sending an e-mail (antispam defense activated)
Title: C per applicazioni numeriche
Author: Sandro Tosi
Last modified: 30/01/2004

Cerchero' qui di fare un  breve vademecum per utilizzare il linguaggio
di  programmazione C  per sviluppare  applicazioni rivolte  al calcolo
numerico. Queste  note si  basano sull'esperienza maturata  durante la
tesi di laurea su sistemi  GNU/Linux, Gentoo e Debian, con compilatore
gcc-3.3; pertanto, al variare del sistema operativo o del compilatore,
alcune  indicazioni  potrebbero risultare  inesatte  e errate:  _state
attenti_!

***Utilizzo di funzioni matematiche

Il  C fornisce  una  nutrita  serie di  routine  matematiche che  sono
dichiarate nel  file `math.h'  ma, mentre per  alcune di  essere (come
`fabs',  che  calcola  il  valore  assoluto di  un  numero  reale)  e'
possibile farvi  riferimento senza  problemi, per altre  (come `sqrt',
che calcola la radice quadrata  di un numero) e' necessario aggiungere
il  parametro  di compilazione  `-lm'  che  indica  al compilatore  di
linkare le librerie matematiche.

***I vettori

Per  chi  proviene da  altri  linguaggi  di  programmazione (come  per
esempio Fortran  o Pascal) questa  e' una novita',  ma in C  i vettori
dichiarati  di  n   elementi,  partono  da  0  ed   arrivano  ad  n-1:
dichiarando, come esempio,

int x[n];

potro' accedere agli elementi da x[0] ad x[n-1];

***I vettori multidimensionali

Per  dichiarare un  array multidimensionale,  si usa  la dichiarazione
classica, accostando tra loro  tutte le dimensioni: per dichiarare una
matrice di n x m elementi di interi si usa

int a[n][m]

Si tenga  presente che in questo modo  si e' definito un  array di `m'
elementi a loro volta array di `n' elementi.

***Utilizzo di vettori molto grandi

Nel caso  si abbia  la necessita' di  utilizzare vettori  molto grandi
(diciamo   con   n=2x10^6  float),   non   e'  possibile   dichiararli
staticamente come

float x[n]

poiche', una volta compilato, l'eseguibile risultera' tanto grande che
il  sistema  si  rifiutera'  di  caricarlo  completamente  in  memoria
inviando all'applicazione il messaggio SEGFAULT.

Dunque,  per poterne  far uso  si devono  utilizzare le  istruzioni di
gestione dinamica della memoria:  `malloc' e `calloc'; vediamo come si
usano:  si deve  definire un  puntatore al  tipo di  dati  del vettore
(continuando l'esempio precedente)

float *x;

e successivamente si chiama una delle due funzioni in questo modo:

x = (float *) malloc( n * sizeof(float) )

oppure

x = (float *) calloc( n, sizeof(float) )

dove  si esegue  il cast  esplicito al  tipo del  puntatore  in quanto
entrambe  restituiscono  un  puntatore  ad  un area  di  memoria  `non
tipata'.  `malloc' richiede come  parametro la  dimensione dell'intera
area da allocare, ottenuta come il numero di elementi moltiplicato per
la  loro   dimensione  (si  e'   utilizzato  `sizeof'  in   quanto  su
architetture differenti  la dimensione dei tipi  puo' variare), mentre
`calloc'  chiede  2  parametri,  il  numero  di  elementi  e  la  loro
dimensione.

Le  differenze principali,  che  fanno preferire  per le  applicazioni
numerice  `calloc'   rispetto  a  `malloc'  (o  in   generale  per  la
dichiarazione di array),  sono il fatto che la  prima azzera l'area di
memora che alloca, mentre la seconda  non lo fa, e che `calloc' alloca
lo spazio di memoria in modo contiguo, importante per l'allocazione di
array e matrici in applicazioni numeriche.

Una volta eseguite queste istruzioni, e' possibile accedere al vettore
come al  solito: quindi, per  accedere al decimo elemento  del vettore
`x' si usa

x[9]

IMPORTANTE:  terminato di  utilizzare il  vettore dinamico,  si _DEVE_
liberare la  memoria allocata  tramite l'istruzione `free'  che prende
come parametro il puntatore all'area di memoria da liberare:

free(x)

nel caso non  si liberi la memoria occupata  da vettori dinamici, essa
rimarra' occupata fino al termine del programma principale.

***Vettori & parametri

Un vettore  e', fondamentalmente, un  puntatore ad un'area  di memoria
(anche se definito implicitamente  come int x[n]); inoltre i parametri
delle funzioni  C sono tutti  passati per valore: quindi,  sebbene non
sia possibile  modificare il valore  del puntatore (cioe'  spostare il
puntatore  nella  memoria), sara'  possibile  modificare gli  elementi
dell'array. Il puntatore all'array e' passato per valore, ma l'"array"
in toto possiamo dire che e' passato per riferimento.

Non facciamoci poi illudere dal modo di definire una funzione:

int f(int *a)

e

int f(int a[])

sono definizioni  equivalenti, ed  e' ancora consentito  di modificare
gli elementi di `a'.

Attenzione  dunque  che  tutte   le  modifiche  apportate  al  vettore
all'interno  delle  proprie   funzioni  avranno  valenza  globale  sul
vettore.

***Utilizzo delle macro del preprocessore

Nelle applicazioni numeriche  si cerca in ogni modo  di ottimizzare le
operazioni  svolte, spesso  a scapito  della leggibilita'  del codice;
come  sappiamo,  pero',  scrivere  un  programma in  modo  ordinato  e
leggibile consente  una sua successiva modifica in  tempi molto rapidi
oltre a consentire la sua comprensione anche ad una persona che non ha
scritto il codice.

La prima  cosa che si  e' soliti fare  e' creare una funzione  per una
parte di  codice utilizzata spesso  e poi richiamare  questa funzione:
questo  approccio   consente  di  scrivere  codice   piu'  compatto  e
maggiormente  modificabile  e  manutenibile  (una  correzione  od  una
modifica a quella  parte di codice resa funzione  deve essere eseguita
una sola  volta, invece  che in  tutti i posti  in cui  se ne  fa uso,
diminuendo  il tempo  necessario alla  modifica e  la  possibilita' di
introdurre errori).

Purtroppo, questa metodologia male  si adatta ai calcoli numerici, nel
caso queste  funzioni debbano essere richiamate  molte volte: infatti,
la  necessita' di  dover salvare  le  informazioni per  la chiamata  a
funzione ed il loro ripristino introducono dei ritardi inaccettabili.

Per   risolvere   questo   inconveniente,   possiamo   sfruttare   una
particolarita' del C:  le macro del preprocessore (da  ora in poi solo
macro).

Per spiegare  cosa sono le  macro, facciamo l'esempio  delle costanti:
nel  caso  il nostro  programma  faccia  uso  frequente di  un  valore
numerico  fissato,  risulta  piu'  comodo definire  una  costante  che
contenga questo valore  e nel codice far riferimento  a detta costante
invece che  al valore  numerico che essa  contiene; in questo  modo si
riducono  gli  errori  nel  codice  (non si  correra'  il  rischio  di
sbagliare a scrivere  questo numero dove viene usato)  e si diminuisce
il  tempo di  modifica del  programma  (in quanto  si dovra'  soltanto
cambiare valore  alla costante, invece  che a tutte le  occorrenze nel
codice).

Le macro  funzionano in maniera  molto simile: si definisce  una macro
tramite la direttiva del preprocessore

#define  

che associa al nome scelto il  valore che lo segue (per convenzione il
nome scelto  viene scritto in  maiuscolo, consentendo di  avere subito
chiaro quali sono  macro e quali sono funzioni);  quando nel codice si
fara'   riferimento  al  nome   della  macro,   durante  la   fase  di
precompilazione (che,  come il nome lascia intendere,  precede la fase
di compilazione  vera e propria) tutte  le occorrenze del  nome di una
macro vengono sostituite dal valore ad esso associato.

Mentre per una  costante non si notano grandi  benefici, per il nostro
problema iniziale, le macro risultano molto utili.

Infatti, con  lo stesso meccanismo di definizione  visto in precedenza
e'  possibile anche  definire  delle `funzioni':  per  esempio, se  si
avesse  bisogno di  calcolare molte  volte il  minimo tra  due numeri,
potrebbe essere utile definire una macro come questa:

#define MIN(x, y)  ((x) < (y) ? (x) : (y))

per essere richiamata come

z=MIN(a,b)

che verra' poi sostituita dal preprocessore in

z=((a) < (b) ? (a) : (b))

E'  anche possibile  definire macro  piu'  complesse e  su piu'  linee
(tramite  l'uso  del  carattere  `\'  a  fine  linea,  che  indica  di
continuare  sulla riga  successiva nella  dichiarazione  della macro),
come per esempio

#define DIV(x, y)                \
	if ( y==0 )              \
	  { *error handling* }   \
	else                     \
	  { printf("%f", x/y); } 

Puo' essere utile, in fase di verifica del codice, controllare come
vengono sostituite le macro dal preprocessore: per fare questo si usa
l'opzione `-E' del `gcc'

gcc -E -o  

che restituisce in   il risultato del passaggio del preprocessore
sulla lista  di file in , senza  eseguire la compilazione
(il risultato e', quindi, un file di testo).

***Metodo di memorizzazione delle matrici

Questa  e' una  nozione importante  per poter  sviluppare applicazioni
numeriche efficenti:  come? Seguendo l'ordine  di memorizzazione delle
matrici nella memoria utilizzato dal linguaggio che si sta usando.

Supponiamo di dover operare su una matrice 3x3 come questa:

+ - + - + - +
| 1 | 2 | 3 |
+ - + - + - +
| 4 | 5 | 6 |
+ - + - + - +
| 7 | 8 | 9 |
+ - + - + - +

Il  C memorizza la  matrice riga  per riga,  in un  vettore: otterremo
quindi che la nostra matrice avra' questa forma in memoria:

+ - + - + - + - + - + - + - + - + - +
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
+ - + - + - + - + - + - + - + - + - +

se quindi dobbiamo operare sulla matrice sara' conveniente organizzare
le  nostre operazioni secondo  le righe  della matrice  (piuttosto che
secondo le colonne), in modo  da accedere linearmente e senza eseguire
`salti' all'interno del vettore.

Come nota  ulteriore, aggiungiamo che il Fortran  memorizza le matrici
per colonne:

+ - + - + - + - + - + - + - + - + - +
| 1 | 4 | 7 | 2 | 5 | 8 | 3 | 6 | 9 |
+ - + - + - + - + - + - + - + - + - +

***Utilizzare le routine LAPACK in C

Lapack  e' un  pacchetto di  librerie di  algebra  lineare liberamente
utilizzabili scritte in Fortran; dal momento che abbiamo intenzione di
utilizzare queste  librerie all'interno  del nostro codice  C dobbiamo
seguire alcune accortezze per farlo.

Prima di tutto vediamo la  corrispondenza dei tipi delle variabili dei
due linguaggi:

C			Fortran

float			real
double			double precision
struct {float r, i;}	complex
struct {double r, i;}	double complex

per i complessi, che in Fortran sono built-in, in C si deve utilizzare
una coppia di reali.

Supponiamo ora di  conoscere il nome della routine  Lapack che abbiamo
intenzione di  usare (per  esempio, la routine  per la  risoluzione di
sistemi lineari  di numeri reali in singola  precisione utilizzando la
fattorizzazione  LU  con pivoting  si  chiama  SGESV): per  richiamare
questa  funzione  all'interno   del  nostro  programma  e'  necessario
aggiungere un  underscore `_' al  nome della routine scritto  tutto in
minuscolo, in  quanto il compilatore Fortran  aggiunge l'underscore al
nome  della funzione quando  genera il  codice compilato.  Inoltre, si
deve  dichiararla  come  `extern'   (qualcosa  del  tipo  extern  void
sgesv_(......) ); io ho fatto senza la dichiarazione e non
ho avuto  problemi, ma un  `Warning' mi avvisava che  la dichiarazione
delle funzioni era implicita, ma poi tutto ha funzionato.

Se,  comunque,  si  vuole  rimuovere questo  fastidioso  messaggio  di
avviso, e'  possibile trovare sul  sito di Michael R.  Haggerty, nella
sezione  riguardante la  programmazione  (il riferimento  si trova  in
calce) un archivio contenente  le intestazioni C delle routine Lapack.
Lo stesso Netlib Repository, fornisce  un file di include contenente i
prototipi  delle funzioni  LAPACK per  il pacchetto  CLAPAK  (il link,
anche questa volta si trova in calce).

Fortran usa esclusivamente le chiamate per riferimento, quindi tutti i
parametri  devono essere  passati per  riferimento, cioe'  deve essere
passato  il  loro puntatore:  questo  implica  che  tutti i  parametri
passati  direttamente  (per intenderci,  se  passiamo valori  numerici
diretti, come  in sqrt(4)) devono  essere prima messi in  variabili di
cui  successivamente  si passera'  il  puntatore;  se  i parametri  da
passare  sono  dei  caratteri  (per  esempio,  la  routine  SPOSV  per
risolvere sistemi lineari con matrici simmetriche e definite positive,
richiede un  parametro per sapere  se la matrice e'  memorizzata nella
parte triangolare inferiore, L, o superiore, U) si possono passare tra
doppi apici  direttamente alla routine Lapack  (nell'esempio di prima,
per usarlo in C, si passera' "U" oppure "L").

Veniamo  ora alla  parte un  po'  piu' complessa:  il passaggio  delle
matrici.  Come detto sopra,  Fortran memorizza  le matrici  in vettori
ponendo una intera colonna di  seguito all'altra: vediamo uno pezzo di
codice per chiarire meglio la conversione

float AC[size][size], AF[size*size];

[...]

for (i=0; i  -llapack -lblas -lg2c

g2c e' la  libreria di interfaccia per il C  al Fortran, necessaria in
quanto le librerie Lapack sono scritte proprio in Fortran.

***Link Utili

Using Lapack Routines in C Programs
		http://www.nacse.org/demos/lapack/cprog.html

How Matrices Are Stored In Memory
		http://www.nacse.org/demos/lapack/storage.html

The C Preprocessor (Richard M. Stallman)
		http://www.lns.cornell.edu/public/COMP/info/egcs/cpp/cpp_toc.html

C, MEX, BLAS, LAPACK, Fortran and Java
		http://www.cs.nyu.edu/cs/faculty/overton/g22_nm/programs/info.html

readme for CLAPACK
		http://www.netlib.org/clapack/readme

Include file of C prototypes for CLAPACK
		http://www.netlib.org/clapack/clapack.h

Using Fortran Library Routines in C Programs
		http://www.hpc.sfu.ca/bugaboo/ftnlibs_for_c.html

Michael R. Haggerty (Programming section)
		http://www-heller.harvard.edu/people/mhagger/download/