Matlab : du prototype

1
Matlab du prototype à lapplication

• François Bodin et Stéphane Chauveau
• IRISA/INRIA
• campus de Beaulieu
• 35042 Rennes cédex

2
Sommaire de la présentation

• contexte du projet
• objectifs
• Menhir
• processus de compilation
• description du système cible
• exemples de performances
• perspectives

3
Contexte

• de nombreux prototypes sont écrits en Matlab
• développement rapide mais performances limitées
• ils sont fréquemment recodés
• mises en œuvre des applications numériques
• exploiter les bibliothèques existantes
• exploiter les architectures parallèles
• sabstraire de Fortran
• spécification mathématique

4
MatLab versus Fortran
MatLab
Fortran
for k 1M if ( vv lt toltol ), break, end z
Av t c/(vz) xxtv rr-tz drr
vr(d/c)v cd end
do k1,M if (dot_product(v,v)lttoltol) exit z
matmul(A,v) tc/ dot_product(v,z) xxtv
rr-tz ddot_product(r,r) vr(d/c)v cd
enddo
...
Scalapack
Lapack
bibl. constr.
exécutable
5
Objectifs de Menhir
function x2, error, iter, flag cg(...)
.... r b - Ax error norm( r ) /
bnrm2 .... for iter 1maxit z M \ r
rho (r'z) .... q Ap alpha rho /
(p'q ) .... end ....
6
Objectifs de Menhir (suite)

• génération de code à partir de Matlab 4
• exécution performante
• exploitation du parallélisme
• compilateur paramétrable
• exploiter des bibliothèques optimisées
• produire des implémentations variées (C et
  Fortran)
• contrôler la génération de code
• ensemble de directives

7
Menhir (Matlab ENvironment for High peRformance)

• compilateur  reciblable  pour Matlab 4

script
Menhir
fichier(.f.c)
M-file
système cible

• processus de compilation dun code Matlab
• analyse des types
• génération de code

8
Analyse des types
aucune information sur le type de a,b et p ne
peux être déduite (par défaut mis en œuvre par
des matrices de complexes)

• réduire les tests à lexécution
• peu dinformations dans les programmes
  (opérateurs polymorphes)
• utilisation de directives
• interprocédurale avec clonage

function ppgcd(a,b) while (a b) if (a gtb) a
a-b else b b-a end end p a
Sémantique dynamique des opérateurs I
masque, indirection avec ou sans
redimensionnement
A(I) B
9
Analyse des types

• doit déterminer
• la nature des identificateurs
• la forme (scalaire, matrice, ...),
• la taille des dimensions,
• le type de valeur (complexe, réelle, ...)

script res.m
constante p
res a pi pi 3 res pi
A 1 2 for k 3bsup A(k) 0 end
constante 3
la taille de A augmente à chaque itération
10
Analyse des types

• les directives indiquent les propriétés des
  objets
• VAR var ,var no RESIZE
• VAR var ,var no MASK
• RANGE var int,...
• PRAGMA var ,var shape
• ...
• attention Matlab directives ! F90

11
Génération de code
f1 rend un vecteur dans une structure de type
matrice

• trois tâches principales interdépendantes
• sélection des classes (structures de données)
• sélection des méthodes de calculs
• calcul des conversions

norm(complex v) matriceR f1(..) ... m
f1(...) norm (m)
mat. réelle
...
vecteur ligne réel
conversion de structure de données
vecteur ligne complexe
12
Exemple de code C généré
if (((((_b).dim1)((_A).dim1))(((_b).dim2)
((_A).dim2)))) __tmp94 get_VRi((_x),1) __t
mp126 (_b).dim __tmp127 (_b).dim1 __tmp128
(_b).dim2 alloc_R(_r, __tmp127,
__tmp128) __tmp129 (_b).data __tmp130
(_b).ld - (_b).dim1 __tmp131
(_A).data __tmp132 (_A).ld -
(_A).dim1 __tmp133 _r.data __tmp134 _r.ld
- _r.dim1 for (__tmp136 1 __tmp136 lt
__tmp128 __tmp136) for (__tmp135 1
__tmp135 lt __tmp127 __tmp135)
(__tmp133) ((__tmp129)-((__tmp131)__tm
p94)) __tmp129 __tmp129__tmp130 __tmp131
__tmp131__tmp132 __tmp133
__tmp133__tmp134
13
Description du système cible

• mises en œuvre des variables
• classe (structure de données et leurs propriétés)
• déclarations
• accès aux données (agrégats, élémentaires, locaux
  et globaux)
• conversions entre structures de données
• mises en œuvre des opérateurs
• constructions de base du langage
• fonctions Matlab

14
Description des structures de données
même structure de donnée
propriétés de PosInt

• déclaration
• orientée objet
• comporte les champs
• elem, shape
• minreal, minimag
• maxreal, maximag
• prop

class PosInt elem real minreal 0 prop
int end class MatReal shape matrix elem
real end class ColumnAbsintMatReal shape
column elem PosInt end
inline _at_declare(MatReal MAT) gt StructName MAT
spécialisation en matrice colonne avec des
éléments entiers positifs
15
Description des instructions et opérateurs

• la déclaration de la mise en œuvre des fonctions
  (ici Mat rand(m,n))
• produit le code

inline rand_R2(out MatReal MAT,int DIM1,int
DIM2, local int v1, int
v2) gtalloc_R(MAT,DIM1,DIM2) gtfor (v11
v1ltDIM1 v1) for (v21 v2ltDIM2 v2) gt
MAT.datav1v2 rand()
alloc_R(Mat,m,n) for (v11 v1lt m v1) for
(v21 v2lt n v2) Mat.datav1v2
rand()
16
Exemple de Lapack
class boolean int minreal0 maxreal1
end class RowInt elemint shaperow end
inline int _at_dim0(RowInt MAT) "MAT.dim"
... class PMatComplex elemcomplex
shapematrix end inline _at_declar(PMatComplex
VAR) gt PMatComplex VAR
inline _at_alloc(PMatComplex VAR,int DIM1,int
DIM2) gt alloc_PC(VAR,DIM1,DIM2)
... inline _at_op_matdiv(PMatReal RES, PMatReal
PARAM1, PMatReal PARAM2) gt
div_PR_PRPR(RES,PARAM1,PARAM2)
inline _at_op_matdiv(PMatComplex RES, PMatComplex
PARAM1, PMatComplex PARAM2) gt
div_PC_PCPC(RES,PARAM1,PARAM2)
17
Quelques exemples de performances

• codes séquentiels gauss, cholesky, ...
• codes parallèles gradient conjugué avec et sans
  préconditionnement, Jacobi
• comparaison de performances
• Matlab 4.2, Matlab 5
• MCC
• Menhir avec Lapack et Scalapack
• deux architectures cible
• Sparc station 5, 110 Mhz
• Onyx 4 processeurs MIPS R10000, 194Mhz

18
Codes séquentiels

• gauss élimination de Gauss (matrice 100x100)
• chol factorisation incomplète de Cholesky
  (matrice 300)
• pfl méthode itérative de recherche de
  sous-espace propre par dichotomie spectrale
• INT 1, 4 intégration déquations différentielles

• tailles Matlab C
• gauss 70 200
• chol 180 620
• pfl 200 9000
• INT 1 140 4700
• INT 4 200 3100

pivot 0 for jin temp abs(A(j,i))
if (pivotlttemp) pivottemp
ipivotj end end
19
Codes séquentiels (suite)
20
Gradient conjugué
for k1M if (norm(v)lttol) ,break ,end zAv
tc/(v'z) xxtv rr-tz
dnorm(r)2 vr(d/c)v cd end

• algorithme du gradient conjugué (tailles de
  problème 936 et 1408)

21
Gradient conjugué avec préconditionnement
for iter 1max_it z M \ r
rho (r'z) if ( iter gt 1 ),
beta rho/rho_1 p z betap
else p z end q Ap
alpha rho / (p'q ) ...

• algorithme du gradient conjugué (tailles de
  problème 481 et 936) avec préconditionnement
  (matrice triangulaire)

22
Jacobi

• résolution de systèmes linéaires par la méthode
  de Jacobi
• Jacobi 1 test diagonale dynamique
• Jacobi 2 test statique

... msize(A,1) nsize(A,2) M, N
split( A , b, 1.0, 1 ) for iter 1max_it,
x_1 x y Nx b x M \
(y) error norm(x-x_1
)/ norm( x ) if ( error lt tol ), break,
end end ...
23
Autres projets fondés sur Matlab

• compilateur Matlab
• MCC
• Falcon
• Matcom
• Conlab
• Ramaswamy et al.
• environnements interactifs
• Octave
• Scilab
• ...

• extensions parallèles de Matlab
• Scilab//
• MultiMatlab
• MathServer

24
Perspectives

• gestion des matrices creuses
• exploitation du parallélisme des boucles Matlab
• exploitation des propriétés algébriques
• développement conjoint Matlab/Fortran