ScaLAPACK

De Wiki de Calcul Québec
Aller à : Navigation, rechercher
Cette page est une traduction de la page ScaLAPACK et la traduction est complétée à 100 % et à jour.

Autres langues :anglais 100% • ‎français 100%

Sommaire

Description

ScaLAPACK (Scalable Linear Algebra Package) is a mathematical library that contains a set of functions for solving systems of linear equations and eigenvalue problems. These functions can be used with full matrices, symmetric or not, or with banded matrices. This library can be considered as a parallel version of LAPACK. In fact, the function names are similar to be able to easily go from one to the other. Multiple functions exist in three forms: QR, LU and "Divide and Conquer". The basic (single-core) QR function is slower, however it is the algorithm with the highest parallel efficiency (almost linear). Hence, if you have many cores, its performance approaches that of the others. From experience, the divide-and-conquer functions are the fastest.

Documentation and the reference version of the library are available online at Netlib. ScaLAPACK requires an existing installation of LAPACK, BLAS and BLACS. Usually a programmer who wants to use ScaLAPACK only needs to call ScaLAPACK functions and a few BLACS functions.

Available implemenations

Intel Math Kernel Library

Math Kernel Library contains a complete version of ScaLAPACK. It is easy to use jointly with the LAPACK, BLAS, and BLACS libraries within the same package.

Example

Here is an example written in Fortran90 to calculate the five smallest eigenvalues of a matrix and their corresponding eigenvectors using Scalapack's help.

File : ex_scalapack.f90
!  Steve Allen, CCS, March 2006
!
!  Program to calculate five eigenvalues of a matrix of size
!  2000 x 2000. This program uses the Scalapack library.

 
  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
!
! List of intrinsic and external functions
!
  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
 
! Function to generate the matrix
  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
!
! Intrinsic and external functions
  integer,external                 :: numroc
  integer,intrinsic                :: size, mod
  real,intrinsic                   :: real, sqrt
 
!  Collect info on the dimensions and evaluate the size of the
!  local matrix
  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
 
!   Put this information in the desc array
!
  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
 
!  Generate the local components of the matrix
!
  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


To compile this example using the Math Kernel Library, you need to use a compilation command like this:

[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