Intel Math Kernel Library (Intel MKL) : Ensemble de bibliothèques mathématiques fortement optimisées et thread-safe pour les applications scientifiques, techniques et financières nécessitant des performances de haut niveau.

La bibliothèque Intel Math Kernel (MKL) installée au mésocentre prend en compte les routines suivantes :

BLAS (Basic Linear Algebra Subprograms) : Level 1,2 et 3
LAPACK : toutes les routines LAPACK version 3
Sparse Solver (Solveur Creux) : Routines pour la résolution de Systèmes linéaires creux
Vector Math Library (VML) : Permet des opérations de type tan, sqrt,exp, sur des vecteurs
Vector Statistical Library (VSL) : Génération de nombre aléatoires
FFT
ScaLAPACK
BLACS

MKL possède une interface C pour BLAS

Options de compilation des applications MKL : http://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/

Depuis la version 13, la compilation des applications MKL est rendu plus simple :

-mkl=sequential pour les applications MKL non threadés exemple (Blas1,2,3)
-mkl=parallel (valeur par defaut) pour activer les support multithreads exemple (Lapack)
-mkl=cluster pour activer le mode cluster i.e MPI, BLACS, exemple (scalapack)

Les variables d'environnements relatives à MKL sont chargés dans votre environnement en chargeant le module intel/14.0.2

$ module load intel/14.0.2

Exemple : multiplication matricielle en utilisant la routine BLAS DGEMM

matmul.f90

 program blas_matmul
c declare variables
      integer dim1, dim2, dim3
      parameter (dim1 = 1000, dim2 = 1000, dim3 = 1000)
      double precision A(dim1, dim2), B(dim2, dim3), C(dim1, dim3)
      double precision alpha, beta
      real start, finish, tarray(2)
      alpha = 1.0d0
c init data      
      beta = 0.0d0
      call srand(86456)
      do i = 1, dim1
        do j = 1, dim2
          A(i, j) = rand()
        enddo
      enddo
      do i = 1, dim2
        do j = 1, dim3
          B(i, j) = rand()
        enddo
      enddo
c etime is not a part of the standard 
      call etime(tarray, start)
      call dgemm('N','N',dim1,dim3,dim2,alpha,A,dim1,B,dim2,beta,C,dim1)
      call etime(tarray, finish)
      print *,'time for C(',dim1,',',dim3,') = A(',dim1,',',dim2,') B(',
     1dim2,',',dim3,') is',finish - start,' s'
      end program blas_matmul

Compilation

%ifort matmul.f90 -o matmul_mkl -mkl=sequential

Exemple : utilisation de la routine DGETRI pour calculer l'inverse d'une matrice, à partir de la décomposition LU de celle-ci.

Paramètres :

A : Matrice (décomposition LU de) d'ordre N

WORK : Tableau de travail (de taille LWORK, LWORK doit être au moins égal à N)

IPIV : Tableau d'entiers (liste des permutations, N éléments)

LDA : Entier (Nombre effectif de lignes de A)

INFO : Entier (Retour d'information sur le déroulement des opérations)

En sortie, A contient l'inverse de la matrice originelle et INFO = 0 si le processus s'est correctement déroulé.

dgetri.f

	PROGRAM TEST_DGETRI
C...Inversion d'une matrice
	IMPLICIT NONE
	DOUBLE PRECISION MAT,MATINV,M	
	DIMENSION MAT(4,4),MATINV(4,4),M(4,4)
C...
	INTEGER INFO
	INTEGER IPIV
	DIMENSION IPIV(4)
C...Espace de travail pour DGETRF
	DOUBLE PRECISION WORK
	DIMENSION WORK(80)	
 
	INTEGER I,J
 
C...Elements de la matrice MAT 
C...Premiere ligne
	MAT(1,1)=3
	MAT(1,2)=-1
	MAT(1,3)=2
	MAT(1,4)=-4
C...Deuxieme ligne	
	MAT(2,1)=5
	MAT(2,2)=-8
	MAT(2,3)=4
	MAT(2,4)=-1
C...Toisieme ligne
	MAT(3,1)=5
	MAT(3,2)=6
	MAT(3,3)=1
	MAT(3,4)=7
C...Quatrieme ligne
	MAT(4,1)=1
	MAT(4,2)=2
	MAT(4,3)=-8
	MAT(4,4)=7
 
	WRITE(6,*)'MAT:'
	DO I=1,4
	WRITE(6,42)(MAT(I,J),J=1,4)
	ENDDO
 
C...Copie de MAT dans MATINV, via DCOPY (routine BLAS)
	CALL DCOPY(4*4,MAT,1,MATINV,1)
 
C...Decomposition LU de MATINV (==MAT pour l instant)
C...DGETRF( M, N, A, LDA, IPIV, INFO )
	CALL DGETRF(4,4,MATINV,4,IPIV,INFO)
C.. en sortie, MATINV contient la decomposition LU
C.. en sortie, INFO=0 si tout va bien
C.. en sortie, IPIV contient la liste des permutations
	WRITE(6,*)'EN SORTIE DE DGETRF: INFO=',INFO
	DO I=1,4
	WRITE(6,*)'IPIV(',I,')=',IPIV(I)
	ENDDO
 
C...Inversion
C...DGETRI(N,A,LDA,IPIV,WORK,LWORK,INFO)
C...Rem:: il faut au minimum LWORK >= N
	CALL DGETRI(4,MATINV,4,IPIV,WORK,80,INFO)
C.. en sortie, MATINV contient l inverse de MAT
C.. en sortie, INFO=0 si tout va bien
	WRITE(6,*)'EN SORTIE DE DGETRI: INFO=',INFO
	WRITE(6,*)'MATINV:'
	DO I=1,4
	WRITE(6,42)(MATINV(I,J),J=1,4)
	ENDDO
 
C..Verification; produit MAT*MATINV, via DGEMM (routine BLAS)
C...Produit matrice matrice C:=ALPHA*(A*B)+BETA*C
C...DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
	CALL DGEMM('N','N',4,4,4,1.D0,MAT,4,MATINV,4,0.D0,M,4)
 
	WRITE(6,*)'M=MAT*MATINV:'
	DO I=1,4
	WRITE(6,42)(M(I,J),J=1,4)
	ENDDO
 
 42	FORMAT(4(1X,1PE22.14))
 
	END PROGRAM TEST_DGETRI

Compilation

%ifort dgetri.f -o dgetri -mkl

Pour activer le support multithreading il faut rajouter l'option -mkl=parallel à la compilation

A l'exéxution, il faut spécifier le nombre de threads à utiliser par la commande : export MKL_NUM_THREADS=N

Si aucune valeur n'est donnée, le système utilisera la totalité des coeurs
voir la section : Contrôle des Threads

Exemple : utilisation de la routine PDGEMV(), qui calcul en parallèle (mémoire distribuée) le produit matrice-vecteur y =alphaAx+beta*y.

Avec :

A= [1 4 7 10 13 ; 3 6 9 12 15 ; 5 8 11 14 17 ; 7 10 13 16 19; 9 12 15 18 21]

et

x=[1 ; 1; 0; 0; 1], x est un vecteur colonne

Appel la routine PDGEMV() pour calculer y=Ax

Le résultat est y=[18; 24; 30; 36; 42]

scala.c

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "mpi.h"
#define AA(i,j) AA[(i)*M+(j)]
 
int main(int argc, char **argv) {
 
   int i, j, k;
 
/************  MPI ***************************/
 
   int myrank_mpi, nprocs_mpi;
   MPI_Init( &argc, &argv);
   MPI_Comm_rank(MPI_COMM_WORLD, &myrank_mpi);
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs_mpi);
 
/************  BLACS ***************************/
 
   int ictxt, nprow, npcol, myrow, mycol,nb;
   int info,itemp;
   int ZERO=0,ONE=1;
   nprow = 2; npcol = 2; nb =2;
   Cblacs_pinfo( &myrank_mpi, &nprocs_mpi ) ;
   Cblacs_get( -1, 0, &ictxt );
   Cblacs_gridinit( &ictxt, "Row", nprow, npcol );
   Cblacs_gridinfo( ictxt, &nprow, &npcol, &myrow, &mycol );
 
   int M=5;
   double *AA = (double*) malloc(M*M*sizeof(double));
 
   for(i=0;i<M;i++ )
     for(j=0;j<M;j++)
        AA[i*M+j]=(2*i+3*j+1);
 
   double *X = (double*) malloc(M*sizeof(double));
 
   X[0]=1;X[1]=1;X[2]=0;X[3]=0;X[4]=1;
 
   int descA[9],descx[9],descy[9];
 
   int mA = numroc_( &M, &nb, &myrow, &ZERO, &nprow );
 
   int nA = numroc_( &M, &nb, &mycol, &ZERO, &npcol );
 
   int nx = numroc_( &M, &nb, &myrow, &ZERO, &nprow );
 
   int my = numroc_( &M, &nb, &myrow, &ZERO, &nprow );
 
   descinit_(descA, &M,   &M,   &nb,  &nb,  &ZERO, &ZERO, &ictxt, &mA,  &info);
   descinit_(descx, &M, &ONE,   &nb, &ONE,  &ZERO, &ZERO, &ictxt, &nx, &info);
   descinit_(descy, &M, &ONE,   &nb, &ONE,  &ZERO, &ZERO, &ictxt, &my, &info);
 
   double *x = (double*) malloc(nx*sizeof(double));
   double *y = (double*) calloc(my,sizeof(double));
   double *A = (double*) malloc(mA*nA*sizeof(double));
 
   int sat,sut;
 
   for(i=0;i<mA;i++)
     for(j=0;j<nA;j++){
        sat= (myrow*nb)+i+(i/nb)*nb;
        sut= (mycol*nb)+j+(j/nb)*nb;
        A[j*mA+i]=AA(sat,sut);
        }
 
   for(i=0;i<nx;i++){
    sut= (myrow*nb)+i+(i/nb)*nb;
        x[i]=X[sut];
     }
 
   double alpha = 1.0; double beta = 0.0;
 
   pdgemv_("N",&M,&M,&alpha,A,&ONE,&ONE,descA,x,&ONE,&ONE,descx,&ONE,&beta,y,&ONE,&ONE,descy,&ONE);
 
   Cblacs_barrier(ictxt,"A");
 
   for(i=0;i<my;i++)
   printf("rank=%d %.2f \n", myrank_mpi,y[i]);
 
   Cblacs_gridexit( 0 );
   MPI_Finalize();
   return 0;
}

Compilation : Tout d'abord il faut charger le module ompi

%module load mpi/openmpi

% mpicc scala.c -o scala -mkl=cluster

Exécution

%mpirun -np 4 ./scala
 
rank=1 0.00 
rank=2 30.00 
rank=2 36.00 
rank=0 18.00 
rank=0 24.00 
rank=0 42.00 
rank=3 0.00 
rank=3 0.00 
rank=1 0.00 
rank=1 0.00

Contrôle des Threads

Pour forcer MKL à s'exécuter en mode séquentiel :

export MKL_SERIAL=YES

Il est possible de spécifier le nombre de threads à utiliser par une application. Pour cela, il est nécessaire d'exporter certaines variables d'environnements.

Environment	Variable Service	Function Comment	Equivalent OpenMP*
MKL_NUM_THREADS	mkl_set_num_threads	Suggests the number of threads to use.	OMP_NUM_THREADS
MKL_DOMAIN_NUM_THREADS	mkl_domain_set_num_threads	Suggests the number of threads for a particular function domain.
MKL_DYNAMIC	mkl_set_dynamic	Enables Intel MKL to dynamically change the number of threads.	OMP_DYNAMIC

Valeurs possibles pour MKL_DOMAIN_NUM_THREADS

MKL_DOMAIN_NUM_THREADS accepts a string value <MKL-env-string>, which must have the following format:
 
<MKL-env-string> ::= <MKL-domain-env-string> { <delimiter> <MKL-domain-env-string> }
 
<delimiter> ::= [ <space-symbol>* ] ( <space-symbol> | <comma-symbol> | <semicolon-symbol> | <colon-symbol> ) [ <space-symbol>* ]
 
<MKL-domain-env-string> ::= <MKL-domain-env-name> <uses> <number-of-threads>
 
<MKL-domain-env-name> ::= MKL_ALL | MKL_BLAS | MKL_FFT | MKL_VML
 
<uses> ::= [ <space-symbol>* ] ( <space-symbol> | <equality-sign> | <comma-symbol>) [ <space-symbol>* ]
 
<number-of-threads> ::= <positive-number>

Exemple : Application utilisant VML et BLAS

Objectif : Forcer le nombre de threads à utiliser à 1

export MKL_NUM_THREADS=1
 
export MKL_DOMAIN_NUM_THREADS="MKL_ALL=1, MKL_BLAS=1, MKL_VML=1"
 
export MKL_DYNAMIC=FALSE

Pour plus d'information MKL Threading

Introduction à MKL

Environnement MKL

Utilisation de BLAS MKL avec le compilateur Intel

Utilisation de LAPACK MKL avec le compilateur Intel

Utilisation de ScaLAPACK MKL avec le compilateur Intel

Contrôle des Threads

Liens