ScaLAPACK

De Wiki de Calcul Québec
Aller à : Navigation, rechercher
Autres langues :anglais 100% • ‎français 100%

Sommaire

Description

ScaLAPACK (Scalable Linear Algebra Package) est une bibliothèque mathématique regroupant un ensemble de fonctions pour résoudre des systèmes d'équations linéaires ou des problèmes aux valeurs propres. Ces fonctions peuvent être utilisées avec des matrices pleines, symétriques ou non, ou des matrices en bande. On peut considérer cette bibliothèque comme la version parallèle de LAPACK. En fait, les noms des fonctions sont similaires afin de passer facilement de l'une à l'autre. Plusieurs fonctions existent sous trois formes : QR, LU et « Divide and Conquer ». De base, la version QR est la plus lente, mais c'est celle qui offre le meilleur efficacité parallèle (presque linéaire). Donc, avec un très grand nombre de cœurs, sa performance se rapproche de celle des autres. Par expérience, les fonctions de la catégorie « Divide and Conquer » sont les plus rapides.

La documentation et la version de base de la bibliothèque sont disponibles en ligne sur le site de Netlib. ScaLAPACK requiert l'installation préalable de LAPACK, BLAS et BLACS. Toutefois, le programmeur qui veut utiliser ScaLAPACK a généralement seulement besoin d'appeler des fonctions de la bibliothèque ScaLAPACK et quelques fonctions de BLACS.

Implémentations disponibles

Intel Math Kernel Library

Math Kernel Library contient une version complète de ScaLAPACK. Il est facile de l'utiliser conjointement avec les bibliothèques LAPACK, BLAS et BLACS de ce même paquet.

Exemple

Voici un exemple écrit en langage Fortran90 pour le calcul des cinq plus petites valeurs propres d'une matrice et les vecteurs propres correspondants à l'aide de Scalapack.

Fichier : ex_scalapack.f90
!  Steve Allen, CCS, mars 2006
!
!  Programme pour le calcul de cinq valeurs propres d'une matrice de taille
!  2000 x 2000. Ce programme utilise la bibliothèque Scalapack.

 
  Program ex_scalapack
  implicit none
  character,parameter              :: porder='R', jobZ='V', prange='I'
  character,parameter              :: uplo='U'
  integer,parameter                :: NTOT = 2000
  integer                          :: mypnum, nproc, nprow, npcol
  integer                          :: ictxt, locrow, loccol, m
  integer                          :: i, n, neig, ia, ja, il, iu
  integer                          :: lwork, liwork, nz, iz, jz
  integer,dimension(9)             :: desc
  integer,dimension(NTOT)          :: ifail
  integer,dimension(:),allocatable :: iwork, icluster
  real                             :: vl, vu, abstol, orfac
  real,dimension(NTOT)             :: eigval
  real,dimension(:),allocatable    :: work, gap
  real,dimension(:,:),allocatable  :: Z
  real,dimension(:,:),pointer      :: mat
!
!   Liste des fonctions intrinsèques et externes
!
  integer,intrinsic                :: max, nint, mod, size
  real,intrinsic                   :: sqrt, real
  logical,intrinsic                :: btest
  interface
    function genmat( ictxt, NTOT, mat, desc )
      implicit none
      integer                          :: genmat
      integer,intent(in)               :: ictxt, NTOT
      real,dimension(:,:),pointer      :: mat
      integer,dimension(:),intent(out) :: desc
    end function genmat
  end interface
 
!   setup BLACS grid
  call blacs_pinfo(mypnum, nproc)
  call blacs_get(-1, 0, ictxt)
 
  nprow = nint( sqrt( real( nproc ) ) )
  npcol = nproc / nprow
  if (mypnum == 0) then
    print*," BLACS will setup ",nprow," X ",npcol," processor grid"
    if (npcol*nprow /= nproc) then
      print*," BLACS grid does not cover all processors ",nproc
    endif
  endif
 
  call blacs_gridinit(ictxt, porder, nprow, npcol)
 
!   generate the matrix
  i = genmat(ictxt, NTOT, mat, desc)
  if (i /= 0) then
    call blacs_abort(ictxt, i)
    stop
  endif
 
!    parameters needed by pssyevx
  locrow = size( mat, 1 )
  loccol = size( mat, 2 ) 
  n = NTOT
  ia = 1
  ja = 1
  iz = 1
  jz = 1
  il = 1
  iu = 5
  neig = -1
  abstol = 0.0
  orfac = 10.0
  vl = 0.0
  vu = 0.0
  liwork = 200 * locrow
  lwork = locrow * loccol * 10
  allocate( icluster( 2*nprow*npcol ) )
  allocate( gap( nprow*npcol ) )
  allocate( Z( locrow, loccol ) )
  allocate( iwork( liwork ) )
  allocate( work( lwork ) )
  if(size(iwork)==0.OR.size(icluster)==0.OR.size(gap)==0.OR.size(work)==0)then
    if( mypnum == 0 ) then
      print*,"Unable to allocate workspace needed for pssyevx"
    endif
    call blacs_abort(ictxt, i)
    stop
  endif
  i = 0
 
!   call ScaLAPACK to compute the five smallest eigenvalues
  call pssyevx(jobZ, prange, uplo, n, mat, ia, ja, desc, vl, vu, &
       &  il, iu, abstol, neig, nz, eigval, orfac, Z, iz, jz, &
       &  desc, work, lwork, iwork, liwork, ifail, icluster, gap, i)
 
  if (i == -21 .OR. btest(i,3)) then ! need more workspace
    lwork = -1 * nint( work(0) )
    deallocate(work)
    allocate( work( lwork ) )
    if (size(work) == 0) then
      if( mypnum == 0 ) then
        print*,"Unable to reallocate work[", lwork,"]"
      endif
      call blacs_abort(ictxt, i)
      stop
    endif
    call pssyevx(jobZ, prange, uplo, n, mat, ia, ja, desc, vl, vu, &
         &  il, iu, abstol, neig, nz, eigval, orfac, Z, iz, jz, &
         &  desc, work, lwork, iwork, liwork, ifail, icluster, gap, i)
  endif
 
!   output somethings from the first processor
  if (mypnum == 0) then
    print*," PSSYEV found ",neig," eigenvalues, return code ",i
    if (neig > 0) then
      do i=1, neig
        print*,"E[",i,"] = ",eigval(i),", error=",eigval(i)-i
      end do
    endif
  endif
 
!   check the eigenvectors
  call blacs_gridinfo(ictxt, nprow, npcol, ia, ja);
  if (mypnum == 0) then
    print*," PSSYEV found ",nz," eigenvectors (5 request)"
  endif
  if (nz > 0) then
    do i=1, nz
      iz = mod( i, nprow )
      jz = mod( i, npcol )
      if (iz==ia .AND. jz==ja) then
        vl = Z(mod(i,locrow),i/locrow)
        print*,"Z(",i,",",i,") = ",vl,", error = ",abs(vl-1.)
      endif
    end do
  endif
 
!   cleaning up
  deallocate(work)
  deallocate(Z)
  deallocate(gap)
  deallocate(iwork)
  deallocate(icluster)
  deallocate(mat)
  call blacs_gridexit( ictxt )
  call blacs_exit(0)
  end
 
!   Fonction pour générer la matrice
  function genmat(ictxt, NTOT, mat, desc)
  implicit none
  integer                          :: genmat
  integer,intent(in)               :: ictxt, NTOT
  real,dimension(:,:),pointer      :: mat
  integer,dimension(:),intent(out) :: desc
  integer,parameter                :: brow = 64, bcol = 64
  integer                          :: locrow, loccol
  integer                          :: nprow, npcol, myprow, mypcol, rsrc, csrc
  integer                          :: i, j, ib, jb, ii, jj, ig, jg
!
!  Fonctions intrinsèques et externes
  integer,external                 :: numroc
  integer,intrinsic                :: size, mod
  real,intrinsic                   :: real, sqrt
 
!   Recueille l'info sur la grille et évalue la taille de la
!   matrice locale
  call blacs_gridinfo(ictxt, nprow, npcol, myprow, mypcol)
 
  rsrc = 0
  csrc = 0
  locrow = numroc( NTOT, brow, myprow, rsrc, nprow )
  loccol = numroc( NTOT, bcol, mypcol, csrc, npcol )
 
  allocate( mat(locrow, loccol) )
  if (size(mat) == 0) then
    write(6,*) "Unable to allocate matrix on processor ",myprow," ", mypcol
    genmat = -1
    return
  endif
 
!   Inscrit l'information dans le tableau desc
!
  call descinit(desc, NTOT, NTOT, brow, bcol, rsrc, csrc, ictxt, locrow, i)
  if (i /= 0) then
    print*,'Problème i =', i
    genmat = -2
    return
  endif
 
!   Génère les composantes locales de la matrice
!
  do j=1, loccol
    jb = ((j - 1) / bcol) * npcol + mypcol
    jj = mod(j - 1, bcol) + 1
    jg = jb * bcol + jj
    do i=1, locrow
      ib = ((i - 1) / brow) * nprow + myprow
      ii = mod(i - 1, brow) + 1
      ig = ib * brow + ii
      if (ig /= jg) then ! off-diagonal
        mat(i,j) = 1.0e-16 / sqrt(real(ig*jg))
      else
        mat(i,j) = real(ig)
      endif
    end do
  end do
  genmat = 0
  end function genmat


Pour compiler l'exemple en utilisant le paquet Math Kernel Library, vous aurez besoin d'une commande de compilation qui ressemble à ce qui suit :

[nom@serveur $] mpif90 ex_scalapack.f90 -lmkl_scalapack_lp64 -lmkl_blacs_openmpi_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm


Outils personnels
Espaces de noms

Variantes
Actions
Navigation
Ressources de Calcul Québec
Outils
Partager