TD/TP 5 de Calcul num´erique M´ethodes it´eratives de base

R´esolution de l’´equation de la chaleur en 1D stationnaire

L’objectif de ce TD/TP est d’appliquer un ensemble d’algorithmes ´etudi´es en cours et TD pour la r´esolution d’un syst`eme lin´eaire obtenu par discr´etisation par la m´ethode des diff´erences finies de l’´equation de la chaleur 1D stationnaire. Les impl´ementations seront faites en C avec BLAS et LAPACK. Afin de valider vous pouvez vous reposer sur les impl´ementations Scilab r´ealis´ees en TD. Une analyse critique de vos r´esultats est demand´ee `a partir de l’´etude des complexit´es en temps et en espace. Pour chaque algorithme,  le temps d’ex´ecution en fonction de la taille de la matrice devra  ˆetre mesur´e et des courbes de performances devront  ˆetre pr´esent´ees.   Un rapport ainsi qu’un lien vers un d´epˆot git doivent ˆetre fourni et rendu `a la fin de la derni`ere s´eance.

Langage C, Biblioth`eque BLAS et LAPACK, Gnuplot

1    Travail pr´eliminaire : Etablissement d’un cas de test

Soit `a r´esoudre l’´equation de la chaleur dans un milieu immobile lin´eaire homog`ene avec terme

source et isotrope :

, ìk = g, x 冬 ]0, 1[

．  T (1) = T1

O`u g est un terme source, k > 0 le coefficient de conductivit´e thermique, et T0  < T1  les temp´eratures au bord du domaine consid´er´e.

On se propose de r´esoudre cette  ´equation par une m´ethode de diff´erence finie centr´ee d’ordre 2. On discr´etise le domaine 1D selon n + 2 noeuds xi , i = 0, 1, 2, ...n + 1, espac´es d’un pas h constant.

En chaque noeud l’´equation discr`ete s’´ecrit :

ìk ╱ 、i = gi

Pour la suite du TP on notera ce syst`eme :

Au = f, A 冬 Rn xn , u, f 冬 Rn

(2)

(3)

Dans la suite du TP on consid`ere qu’il n’y a pas de source de chaleur, i. e .  g = 0.  La solution analytique de cette ´equation est :

T (x) = T0 + x(T1 ì T0)                                                   (4)

ì Exercice 1. Fait en TD

1. Approximer la d´eriv´ee seconde de T au moyen d’un sch´ema centr´e d’ordre 2.

2. E´ crire le syst`eme lin´eaire de dimension n correspondant au probl`eme 1.

3.  Pr´esentez le probl`eme `a r´esoudre et sa discr´etisation en introduction du rapport.

Afin de calculer efficacement la solution de ce probl`eme nous souhaitons d´evelopper un code en langage C faisant appel aux biblioth`eques de calcul BLAS et LAPACK.

ì Exercice   2.     Mise  en  place  de  l’environnement  de  d´eveloppement   (C  + BLAS/LAPACK)

1.  Cr´eez un projet et le versionner avec git.  Cr´eer un d´epˆot public sur la plateforme de votre choix.

2.  V´erifiez l’installation  et  la pr´esence  des  biblioth`eques  CBLAS  et  CLAPACK  (dans /usr/local/lib  et  /usr/local/include)   (sous   Ubuntu   :    apt-get  install  libblas-dev liblapack-dev )

3. E´ crivez un makefile et un code de test de votre environnement de travail (on pourra utiliser le fichier testenv.c)

Dans la suite de ce TP, vous  ´etudierez et compl`eterez le code C associ´e `a ce TP faisant appel `a BLAS et LAPACK et t´el´echargeable sur moodle.

Dans la suite du TP il est conseill´e de cr´eer un fichier lib poisson1D.c qui contient toutes les fonctions d´evelopp´es pour r´esoudre le probl`eme de Poisson1D. Et de cr´eer des fichiers s´epar´es pour les codes de r´esolution.

2    M´ethode directe et stockage bande

On consid`ere dans notre ´etude des matrices tridiagonales. Une m´ethode de stockage efficace pour ce type de matrice est le stockage par bande.

Lors de l’utilisation de LAPACK et BLAS, diff´erents stockages bande sont propos´es. Nous ´etudions dans un premier temps le stockage General Band.

Une matrice de dimension m K n avec kl sous-diagonales et ku sur-diagonales peut ˆetre stock´ee de mani`ere compacte dans un tableau `a deux dimensions avec kl + ku + 1 lignes et n colonnes.  Les colonnes de la matrice sont stock´ees dans les colonnes correspondantes du tableau, et les diagonales de la matrice sont stock´ees dans les lignes du tableau. Ce stockage ne doit ˆetre utilis´e que si kl, ku ← min(m, n). Dans LAPACK les matrices ayant un stockage de ce type ont un nom se teminant par B. Il suit que l’´el´ement aij de la matrice A est stock´ee dans AB(ku + 1 + i ì j, j). Par exemple, pour m = 5, n = 5, kl = 2, ku = 1 : Soit,

sont stockage GB dans le tablea AB est tel que :

╱  、

(    ． .

Les ´el´ements ì doivent ˆetre affect´es. On choisira g´en´eralement 0.

Dans la suite de ce TP, vous concevez et validez un code C faisant appel `a BLAS et LAPACK.

R´ef´erences :

Matrix storage scheme: http://www .netlib .org/lapack/lug/node121 .html Band Storage: http://www .netlib .org/lapack/lug/node124 .html

BLAS Documentation:https://netlib.org/blas/

LAPACK Documentation: http://www .netlib .org/lapack

The LAPACKE C Interface to LAPACK: http://www .netlib .org/lapack/lapacke CLAPACK The Fortran to C version of LAPACK:https://netlib.org/clapack/

ì Exercice 3. R´ef´erence et utilisation de BLAS/LAPACK

A l’aide de la documentation de BLAS et LAPACK r´epondez aux questions suivantes :

1. En C, comment doit on d´eclarer et allouer une matrice pour utiliser BLAS et LAPACK ?

2.  Quelle est la signification de la constante LAPACK COL MAJOR ?

3. A` quoi correspond la dimension principale (leading dimension) g´en´eralement not´ee ld ?

4.  Que fait la fonction dgbmv ? Quelle m´ethode impl´emente-t-elle ?

5.  Que fait la fonction dgbtrf ? Quelle m´ethode impl´emente-t-elle ?

6.  Que fait la fonction dgbtrs  ? Quelle m´ethode impl´emente-t-elle ?

7.  Que fait la fonction dgbsv ? Quelle m´ethode impl´emente-t-elle ?

8.  Comment calculer la norme du r´esidu relatif avec des appels BLAS ?

ì Exercice 4. Stockage GB et appel  `a DGBMV

1. E´ crire le stockage GB en priorit´e colonne pour la matrice de Poisson 1D.

2.  Utilisez la fonction BLAS dgbmv avec cette matrice.

3.  Proposez une m´ethode de validation.

ì Exercice 5. DGBTRF, DGBTRS, DGBSV

1.  R´esoudre le syst`eme lin´eaire par une m´ethode directe en faisant appel `a LAPACK.

2. E´ valuer les performances.  Que dire de la complexit´e des m´ethodes appel´ees ? ì Exercice 6. LU pour les matrices tridiagonales

1. Impl´ementez la m´ethode de factorisation LU pour les matrices tridiagonales avec le format GB.

2.  Proposez une m´ethode de validation.

3    M´ethode de r´esolution it´erative

Nous souhaitons maintenant r´esoudre l’´equation de la chaleur par une m´ethode it´erative. Dans les exercices suivant nous  ´etudions diff´erents algorithmes et leur conception pour une matricie tridiagonale et plus particuli`erement pour ce probl`eme.   Dans un premier temps on effectue une  ´etude th´eorique puis un maquettage avec Scilab.  Dans un second temps on d´eveloppe le code C reposant sur l’appel des BLAS.