Les bibliothèques scientifiques sont des ensembles d’algorithmes testés, validés, optimisés. Elles sont indépendantes de l’architecture machine, et supportent différents types de données : réel, complexe, simple ou double précision. Les constructeurs ont optimisé certaines de ces bibliothèques pour leurs machines (ex : Intel MKL), elles sont plus performantes mais dépendantes de l’architecture.

Il est donc recommandé de les utiliser quand c'est possible, cela permet de se consacrer uniquement aux nouveaux développements.

Exemple : BLAS, LAPACK, FFT, …

BLAS ( Basic Linear Algebra Subprograms) est un ensemble de routines effectuant des opérations élémentaires d'algèbre linéaire de façon optimisée

Plusieurs autres bibliothèques y font référence. On distingue trois niveaux de BLAS:

BLAS 1 : opérations de type vecteur-vecteur

                 y = alpha.x + y

BLAS 2 : opérations de type vecteur-matrice

                   y = alpha.A.x + beta.y

BLAS 3 : opérations de type matrice-matrice

                    C + alpha.A.B + beta.C

Exemple BLAS 1 : la routine sdot (produit scalaire de deux vecteurs)

stod.f

program dot_main
real x(10), y(10), sdot, res
integer n, incx, incy, i
external sdot
n = 5
incx = 2
incy = 1
do i = 1, 10
SDOT = 10.000
  x(i) = 2.0e0
  y(i) = 1.0e0
end do
res = sdot (n, x, incx, y, incy)
print*, `SDOT = `, res
end

Compilation :

$ gfortran stod.f -o stod -lblas

demm.c

#include <stdio.h>
#include <stdlib.h>
 
#include <errno.h> 
 
extern void dgemm_(const char*, const char*, int*, int*,int*, double*, double*, int*, double*, int*, double*, double*, int*);
 
 
int main(int argc, char **argv){
 
  int i,j;
  double one = 1.0;
  double zero = 0.0;
 
  int n = 1024;
 
  if(argc==2)
  {
      n=atoi(argv[1]);
  }
 
 
 
  double *A = (double*)malloc(n*n*sizeof(double));
 
  double *B = (double*)malloc(n*n*sizeof(double));
 
  double *C = (double*)malloc(n*n*sizeof(double));
 
  if(C==NULL || B==NULL || A==NULL){
 
      fprintf(stderr, "Error allocating matrix");
 
      if(A!=NULL) free(A);
      if(B!=NULL) free(B);
      if(B!=NULL) free(C);
 
      return -1;
   }
 
  for(i = 0; i < n*n; i++) {
      A[i] = i;
      B[i] = i;
  }
 
  //make BLAS 3  matrix matrix operation
  //C = C + alpha.A.B + beta.C
  dgemm_("N", "N", &n, &n, &n, &one, A, &n, B, &n, &zero, C, &n);
 
 // #elif defined(CBLAS) || defined(MKL) || defined(VECLIB)
 
 
 
for(i=n-2; i<n;i++){
 
   for(j=n-2;j<n;j++){
 
              printf(" %5.1f ", C[j+n*i]);
       }
     printf("\n");  
}
 
 
 free(A);
 free(B);
 free(C);
 
  return(0);
}

LAPACK ( Linear Algebra PACKage ) est une bibliothèques de collection de routines écrites en Fortran 77 permettent de réaliser les opérations suivantes :

• résolution séquentielle des systèmes d'équations linéaires
• problèmes aux valeurs propres
• factorisations de matrice (QR, …)
• …

Exemple : appel de la routine DGESV pour résoudre un système linéaire.

randomsys.f90

program randomsys
    implicit none
    real(kind=8), dimension(:),allocatable :: x,b
    real(kind=8), dimension(:,:),allocatable :: a
    real(kind=8) :: err
    integer :: i, info, lda, ldb, nrhs, n
    integer, dimension(:), allocatable :: ipiv
 
    ! initialize random number generator seed
    ! if you remove this, the same numbers will be generated each
    ! time you run this code.
    call init_random_seed()  
 
    print *, "Input n ... "
    read *, n
 
    allocate(a(n,n))
    allocate(b(n))
    allocate(x(n))
    allocate(ipiv(n))
 
    call random_number(a)
    call random_number(x)
    b = matmul(a,x) ! compute RHS
 
    nrhs = 1 ! number of right hand sides in b
    lda = n  ! leading dimension of a
    ldb = n  ! leading dimension of b
 
    call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
 
    ! Note: the solution is returned in b
    ! and a has been changed.
 
    ! compare computed solution to original x:
    print *, "         x          computed       rel. error"
    do i=1,n
        err = abs(x(i)-b(i))/abs(x(i))
        print '(3d16.6)', x(i),b(i),err
        enddo
 
    deallocate(a,b,ipiv)
 
end program randomsys

Compilation :

$ gfortran randomsys.f90 -o appli -lblas -llapack

Scalapack (SCAlable Linear Algebra PACKage) est la version scalable de routines issues de la bibliothèque lapack pour architecture à mémoire distribuée

La bibliothèque de FFT (incluse dans ESSL) permet de réaliser les opérations suivantes :

• transformées de Fourier 1D, 2D, 3D complexe/complexe, complexe/réel et réel/complexe,
• produit de convolution et corrélation.

• Blas
• Lapack
• Scalapack
• Documentation intel