Exemple d'utilisation: la fonction gaussband

La résolution d'un système linéaire "Ax=b" dont la matrice est au format Bandmatrix et le second membre ""b" au format Vector peut s'effectuer par une seule fonction quelque soit le type de base car l'algorithme de Gauss est algébrique:

   
template <class T,class R> 
  float gaussband (Bandmatrix<T,R>& a, Vector<R>& x, int first)
{
   int i,j,k, n = a.size, bdthl = a.bdth;    
   T s, s1;
   R s2;
   float smin = 1e9, eps = 1/smin;

   if (first)    // factorization      
     for (i=0;i<n;i++) 
      {
         for(j=mmax(i-bdthl,0);j<=i;j++)
         {    s=0; for (k=mmax(i-bdthl,0); k<j;k++)s += a(i,k)*a(k,j) ; 
            a(i,j) -= s ; 
         }
         for(j=i+1;j<=mmin(n-1,i+bdthl);j++)
         {
             s= a(i,j); 
            for (k=mmax(j-bdthl,0);k<i;k++) s -= a(i,k)*a(k,j) ;
            s1 = a(i,i);
            if(s1.modul() < smin) smin=s1.modul();
            if(s1.modul() < eps) s1 = eps;
            a(i,j) = s / s1;
         }
      }
    for (i=0;i<n;i++)                     // resolution
      {    s2 = x[i]; 
         for (k=mmax(i-bdthl,0);k<i;k++) s2 -= a(i,k) * x[k]; 
         x[i] = s2 / a(i,i) ; 
      }
    for (i=n-1;i>=0;i--)
      { s2=0; for (k=i+1; k<=mmin(n-1,i+bdthl);k++) s2 += a(i,k)* x[k];
       x[i] -= s2 ;
      }
   return smin;   
}

Le paramètre "first" permet de réutiliser la fonction pour la résolution du système linéaire dont la matrice aurait été factorisé lors d'un appel précédent.

Prendre le fichier TNFEMX.TAR sur le web dans http://www.ann.jussieu.fr/ pironneau/ sur le site ftp. Remplacer dans le makefile le mot tnstokes.o par exo.o et utiliser le fichier exo.cpp contenant:

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include "config.h"
#include "aarray.h"
#include "tnvect.h"

const int N=1;                    // Scalar base type follows
typedef float Basetype;           // real numbers (not complex)
typedef VectN<Basetype,N> vectNB; // base type is float (vector of size 1)
typedef MatN<Basetype,N> matNB;   // same for matrices

void main()
{     int nv = 1000, bdth = 3;
      Vector<vectNB>  sol(nv);
      Bandmatrix<matNB,vectNB> aa(nv,bdth);
      time_t now1,now2; time(&now1);  // initialize the clock
      cout<< "pivot=" << gaussband(aa,sol,1) << endl;      
	     time(&now2); cout<<"time ellapsed "<<now2-now1<<endl;
}

Prendre comme exemple le cas ou aa(i,i) = 2, aa(i,i+1) = aa(i,i-1) = -1 et sol[i]=1, sol = aa*sol; de sorte que la solution sera identiquement 1.

Ensuite faire un cas bloc 2x2 sur le meme exemple mais avec nv moitié, aa(i,i) = b, aa(i,i+1) = c, aa(i,i-1) = d et sol[i]=(1,1), sol = aa*sol, avec b=(2,-1;-1,2) du type matNB, c=(0,0;-1,0), et d= c transposé. Comparer les temps d'exécution.