HPC-Wiki (temp)

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
======  ScaLAPACK  ======

=====  Parallel Linear Algebra with ScaLAPACK  =====

ScaLAPACK is a library of high-performance linear algebra routines for parallel distributed memory machines. On GWDG's scientific computing cluster ScaLAPACK is available from the Intel® Math Kernel Library, currently installed in version 11.1. ScaLAPACK provides often used algorithms for dense and banded matrices: the solution of systems of equations, of east squares problems, eigenvalue problems, and singular value problems.

The distributed memory parallelisme of ScaLAPACK is based on the Basic Linear Algebra Communication Subprograms (BLACS) supporting the 2-dimensional data- and process-patterns used in ScaLAPACK for the efficient parallelization of the matrix algorithms of linear algebra. The BLACS are built on top of a message passing library, which is MPI for the Intel MKL implementation of BLACS. The parallel basic matrix operations used in ScaLAPACK are provided by the Parallel Basic Linear Algebra Subprograms (PBLAS), which in turn at the single process level call the Basic Linear Algebra Subprograms (BLAS) for matrix computations. This software hierarchy of ScaLAPACK is displayed in the following figure: 

{{:wiki:hpc:1200px-scalapack.png|}}

The performance of the ScaLAPACK routines depends on the performance of the underlying implementations of BLAS and the message passing library. With Intel MPI and Intel MKL optimal performance will be achieved on GWDG's clusters.

Documentation for ScaLAPACK and related software can be obtained from the web:
[[http://www.netlib.org/scalapack/|ScaLAPACK]],
[[http://www.netlib.org/scalapack/pblas_qref.html|PBLAS]],
[[http://www.netlib.org/lapack/|LAPACK]],
[[http://www.netlib.org/blas/|BLAS]],
[[http://www.netlib.org/blacs/|BLACS]],
[[http://www.mcs.anl.gov/research/projects/mpi/|MPI]].
Furthermore the  [[https://software.intel.com/en-us/mkl_11.2_ref|reference manual for Intel MKL]] gives detailed descriptions for all routines.


Speed up from parallel computing for linear algebra algorithms can easily be realized by the ScaLAPACK procedures. But even for users familiar with parallel applications based on MPI the calling of ScaLAPACK procedures needs some explanation. This concerns the use of the additional communication interface of the BLACS and the special block-cyclically distributed data layout for matrices, which is the basis for a good load balance in the parallelized matrix multiplication and matrix factorization operations of linear algebra algorithms.

=====  First Steps for Using ScaLAPACK on GWDG Compute Clusters  =====

====  A Parallel HELLO_WORLD Program with BLACS  ====
The following example code shows how to set up the BLACS environment needed for calling ScaLAPACK procedures:

<code>
      program hello_from_BLACS
      implicit none
!  local variables
      integer           info, nproc, nprow, npcol,
     :                  myid, myrow, mycol,
     :                  ctxt, ctxt_sys, ctxt_all

! determine rank of calling process and the size of processor set
      call BLACS_PINFO(myid,nproc)
! get the internal default context
      call BLACS_GET( 0, 0, ctxt_sys )
! Set up a process grid for the process set
      ctxt_all = ctxt_sys
      call BLACS_GRIDINIT( ctxt_all, 'C', nproc, 1)
      call BLACS_BARRIER(ctxt_all,'A')
! Set up a process grid of size 3*2
      ctxt = ctxt_sys
      call BLACS_GRIDINIT( ctxt, 'C', 3, 2)
! All processes not belonging to ctxt jump to the end of the program
      if (ctxt.lt.0) go to 1000
! Get the process coordinates in the grid
      call BLACS_GRIDINFO( ctxt, nprow, npcol, myrow, mycol )
      write(6,*) 'hello from process',myid,myrow,mycol, nprow, npcol
! ...
!  now ScaLAPACK or PBLAS procedures can be used
! ...
 1000  continue
! return all BLACS contexts
      call BLACS_EXIT(0)
      stop
      end
</code>

For the MPI based implementation of the BLACS, the initial call to ''BLACS_PINFO'' is equivalent to calling the MPI routines ''MPI_INIT'', ''MPI_COMM_SIZE'' and ''MPI_COMM_RANK''. Therefore subsequent calls to all MPI routines are possible (don't forget to include the MPI file '''mpif.h'''), except to ''MPI_INIT'', which would stop the program with the MPI error message ''Cannot call MPI_INIT or MPI_INIT_THREAD more than once''. Calling ''BLACS_PINFO(myid,nproc)'' returns in ''myid'' the task id of the calling process and in ''nproc'' the number of processes, which have been created by starting the program with the ''mpirun'' command.

The next step is setting up a two-dimensional process grid by calling ''BLACS_GRIDINIT(context,order,nprow,npcol)''. This call maps nprow*npcol processes from the process set created by starting the program to a nprow x npcol grid. The mapping is defined by assigning to a process with linear rank ''myid, 0<= myid< nprow*npcol'' a unique rank tuple ''(myrow,mycol), 0<=myrow<nprow, 0<=mycol<npcol'' in the grid. The processes are placed in the grid column after column, if ''order = 'C''' and row after row if ''order = 'R'''. Processes with linear ranks ''myid>= nprow*npcol'' will be mapped to the rank tuple ''(-1,-1)''. If the number of available processes is smaller than the size of the process grid to be created, the call to ''BLACS_GRIDINIT'' will cause a program abortion.

The ''context'' parameter in the call to ''BLACS_GRIDINIT'' on input has to be a valid system context, which is retrieved by the previous call to  ''BLACS_GET(0,0,ctxt_sys)'', on output it gives an integer handle to the internal representation of the grid context. The context handle is set to -1 in all processes not belonging to the grid.

As the example shows, it is possible to create different grids with their own contexts form the available processor set, such that some processes will have multiple grid coordinates, depending on the context. The first call of ''BLACS_GRIDINIT'' creates a grid of size nproc x 1 with context ''ctxt_all'', which contains all started processes. It is convenient to have all processes belonging to at least one context, because the call to ''BLACS_EXIT'' at the end of the program, which frees all contexts and releases all memory allocated internally for using the BLACS, is effective only for processes belonging to a BLACS context. 

Next a grid of size 3x2 is created with context ''ctxt'' and all processes not belonging to this context jump to the end of the program. The processes of the grid then inquire the grid size and their ranks by calling the ''BLACS_GRIDINFO'' routine and write this information to the standard output device.

====  Compiling and Running the Program on GWDG Compute Clusters ====
In order to compile this program for running on the GWDG compute cluster, the necessary modules for using Intels compiler, MPI and MKL must be loaded by the command

<code>
module load intel/compiler intel-mpi intel/mkl
</code>

The compilation then must include linking the needed libraries, as described in the Intel MKL User's Guide. An example for using the sequential MKL library (one thread per process) is 

<code>
mpiifort -o hello hello_from_BLACS  -L$(MKLROOT)/lib/intel64/ 
                                    -lmkl_scalapack_lp64 
                                    -lmkl_blacs_intelmpi_lp64 
                                    -lmkl_intel_lp64 
                                    -lmkl_sequential 
                                    -lmkl_core
</code>

These libraries are also selected for linking by using the option ''-mkl=cluster'', for example

<code>
mpiifort -o hello hello_from_BLACS -mkl=cluster
</code>

Finally, the program can be started with at least 6 processes interactively with the mpirun command, for example 

<code>
mpirun -n 7 hello
</code>

producing the output

<code>
hello from process           0           0           0           3           2
hello from process           1           1           0           3           2
hello from process           2           2           0           3           2
hello from process           3           0           1           3           2
hello from process           4           1           1           3           2
hello from process           5           2           1           3           2
</code>

With the ''bsub'' command the program can submitted to the queueing system for running on the compute cluster. 

<code>
bsub -q mpi-short -o out.%J -n 7 -a intelmpi mpirun.lsf  ./hello
</code>

More information on using the queueing system on the GWDG comte cluster is provided on the GWDG Wiki under the header [[https://info.gwdg.de/docs/doku.php?id=en:services:application_services:high_performance_computing:running_jobs]].

=====  Block-Cyclic Data Distribution  =====
In ScaLAPACK, all arrays are block-cyclically distributed over the available process grid. This distribution will be explained in detail for vector- and matrix-data in the next two sections.
====  Block-Cyclic Vector Distribution  ====

The block-cyclic distribution of a vector over a linearly arranged processor set is characterized by 4 numbers

<code>
n       size of vector
nb      size of blocks
nproc   size of linear processs set
psrc    process containing the first element of the vector
</code>

In the block-cyclic distribution of a vector the blocks are placed cyclically to the available processes in the process set, starting on process psrc. The blocks are stored contiguously in the local memory of a process.  If n is not divisible by nb, the last block will contain less than nb values. The folowing figures illustrate the block-cyclic distribution for

<code>
n=9, nb=2, nproc=2, psrc=1:
</code>
  
global vector:   {{ :wiki:hpc:vector.png |}}

\\ distributed vector: {{ :wiki:hpc:local_vector.png |}}  



In order to determine the individual sizes of local parts of the vector stored in process ip, the vector size n is split up as

<code>
n = (nblock-1)*nb + nlast,
</code>

where ''nblock=⌈n/nb⌉'' is the total number of blocks, of which ''nblock-1'' are of full size ''nb''  and one of size ''0<nlast≤nb'', the number of remaining elements of the vector. 

The ''nblock'' blocks are distributed in a round-robin order to the ''nproc'' available processes: 

<code>
nblock = (lblock-1)*nproc + lfirst
</code>

such that the first ''0<lfirst≤nproc'' processes will store ''lblock=⌈nblock/nproc⌉'' blocks, the remaining ''nproc-lfirst'' processes storing ''lblock-1'' blocks. The last block stored in the last of the ''lfirst'' processes is of size ''0<nlast≤nb'', all others are of full size ''nb''.   

ScaLAPACK provides the utility function NUMROC (number of rows or columns), which implements this calculation to determine the local number of elements on process ip, given the four numbers n, nb, prsc, nproc. A demo version providing the functionaly of NUMROC is 

<code>
      integer function NUMROC_demo( n, nb, ip, psrc, nproc )
      implicit none
      integer n, nb, ip, psrc, nproc
      integer nblock, nlast, lblock, lfirst, ip0

!  ip0 is the distance of ip to psrc
      ip0 = mod(nproc+ip-psrc,nproc) 
      nblock = (n+nb-1)/nb; nlast = n - (nblock-1)*nb
      lblock = (nblock+nproc-1)/nproc; lfirst = nblock - (lblock-1)*nproc
      if ( ip0.lt.lfirst ) then
         numrock = lblock*nb
      else if ( ip0.eq.lfirst ) then
         numrock = (lblock-1)*nb +nlast
      else 
         numrock = (lblock-1)*nb
      end if
      return
      end
</code>

The maximal number  of local elements in any process is given by 
<code>
if (lfirst.eq.0) NUMROC ≤ (lblock-1)*nb + nlast < lblock*nb
if (lrest.gt.0) NUMROC ≤ lblock*nb.
</code>
<div id="LELM"></div>
Therefore ''LELM=lblock*nb '' is an upper bound  for the  ''NUMROC'' function:

<code>
LELM(n,nb,nproc) = lblock*nb = nb*⌈⌈n/nb⌉/nproc⌉
                 = nb*⌈n/(nb*nproc)⌉ = nb*(⌊(n-1)/(nb*nproc)⌋+1);
                 ≤ ⌊(n-1)/nproc⌋ + nb
</code>

====  Block-Cyclic Matrix Distribution  ====

The block-cyclic distribution of a two dimensioinal array on an two dimensional process grid is a direct product of block-cyclic distribution in each dimension. The relevant parameters for the matrix distribution are

<code>
m,n            number of rows and columns of the matrix
mb, nb         sizes of blocks in columns and in rows
nprow, npcol   number of rows and columns of the two dimensional process grid
rsrc, csrc     row and column index of the process containing the first element of the matrix
</code>

The folowing figures illustrate the block-cyclic matrix distribution with parameters

<code>
m=5,n=7, mb=nb=2, nprow=npcol=2, rsrc=csrc=0:
</code>
  
global matrix: {{ :wiki:hpc:matrix.png |}}

\\ distributed matrix:{{ :wiki:hpc:local_matrix.png |}}

\\ ====  Mapping between Global and Local Indices  ====
In order to transform a global mxn matrix A to a block-cyclically distributed data set, an  element A(i,j) of the global matrix must be stored as an element al(il,jl) of the local matrix al in a process with task tupel (ip,jp). This mapping ''i↔(il,ip)'' depends on block size nb of rows or columns, the number proceses np in a row or column in the process grid and on the value sp for row or column of the process grid containing the first element of the matrix. For given values of ''nb, np, sp'', the mapping ''i→(il,ip)'' is defined by the steps
 
<code>
im = i-1
i1 = im modulo nb
i2 = im - nb*i1
k = i1 modulo np  
rp = i1 - np*k
ip = sp + rp modulo np 
il = k*nb + i2 + 1
</code>

and the reverse mapping ''i←(il,ip)'' by

<code>
rp = ip - sp + np modulo np
ilm = il-1
k = ilm modulo nb
i2 = ilm - nb*k
i1 = k*np + rp
i = i1*nb + i2+1
</code>

These mappings are provided by ScaLAPACK tool functions

<code>
i = INDXL2G(il,nb,ip,sp,np)
il = INDXG2L(i,nb,id1,id2,np)
ip = INDXG2P(i,np,id1,sp,np)
</code>
 
where ''id1,id2'' in ''INDXG2L'' and ''id1'' in ''INDXG2P'' are dummy integer variables, included to unify the calling sequences of these functions.

Supposing, that the values for the element A(i,j)  of the global mxn matrix is given by the function ''matval(i,j)'' the mapping function can be used to place these values sequentially to their storage locations in the local arrays al, as shown in the following code fragment. The leading dimension for the local matrix in this code is given by ''llda =  max(1,NUMROC( m, mb, myrow, 0, nprow ))'' and the source process containing the first element A(1,1) is assumed to be (0,0).

<code>
! copy the values of the global matrix
! to the local parts of the distributed matrix
      llda =  max(1,NUMROC( m, mb, myrow, 0, nprow ))
      do j = 1 , n
         jl = INDXG2L(j,nb,0,0,npcol)
         ipc = INDXG2P(j,nb,0,0,npcol)
         do i = 1 , m
            il = INDXG2L(i,mb,0,0,nprow)
            ipr = INDXG2P(i,mb,0,0,nprow)
            if (ipr.eq.myrow.and.ipc.eq.mycol)
     :             al(il+llda*(jl-1))= matval(i,j)
         end do
      end do
</code>
 
Alternatively, the local parts al of the distributed matrix can be initialized in parallel using the inverse mapping  
<code>
! initialize in parallel the local parts of the distributed matrix
      mll = NUMROC( m, mb, myrow, 0, nprow )
      nll = NUMROC( n, nb, mycol, 0, npcol )
      do jl = 1 , nll
         j = INDXL2G(jl,nb,mycol,0,npcol)
         do il = 1 , mll
            i = INDXL2G(il,mb,myrow,0,nprow)
            al(il+mll*(jl-1))= matval(i,j)
         end do
      end do
</code>

A complete example program for using these mapping routines for initializing the distributed matrix can be copied from [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/use_index_mapping.f|here]].
<div id="array_descriptor"></div>

====  Array Descriptors in ScaLAPACK  ====
In ScaLAPACK the information about the block-cyclic data distribution for an array is encoded in an integer vector with 9 elements, called descriptor. Let ''desc_A'' denote the descriptor for the array A. It's elements have the following meaning:

<code>
desc_A(1) = 1            only possible value for the descriptor type 
desc_A(2) = ctxt_A       BLACS context for the distribution of the global matrix A
desc_A(3) = m_A          number of rows of global matrix A 
desc_A(4) = n_A          number of columns of global matrix A 
desc_A(5) = mb_A         column block size for the distribution of A
desc_A(6) = nb_A         row block size for the distribution of A 
desc_A(7) = rsrc_A       process row in the ctxt_A grid containing the first row of A 
desc_A(8) = csrc_A       process column in the ctxt_A grid containing the first column of A 
desc_A(9) = llda         leading dimension of the local matrix containing elements of A
</code>

The first 8 elements of the descriptor describe global properties of the distributed matrix, they must be equal for all participating processes in the grid context ctxt_A. Only the 9th element describes a local property of the distribution. ''llda'' must be greater or equal than ''NUMROC( m_A, mb_A, myrow, rsrc_A , nprow )'', the number of rows of the global matrix assigned to the calling process. Although there are no descriptor elements containing the sizes of the grid dimensions, this information is available from the descriptor indirectly via the element containing the context handle.

All calls to ScaLAPACK and PBLAS routines involve the descriptors for the distributed matrices participating in operations performed by the subroutines. The descriptors have to be initialized accordingly before calling the routines. The initialization can be done explicitely by assigning the appropriate values to the descriptor elements or by calling the special initialization routine DESCINIT. The calling sequence for DESCINIT is:

<code>
      call DESCINIT( desc_A, m_A, n_A, mb_A, nb_A, rsrc_A, csrc_A, ctxt_A, llda, info )
</code>

Attention must be payed to the order of parameters in the calling sequence, which does not follow completely the numbering order of these values in the descriptor array.
 
If the descriptor ''desc_A'' was initialized successfully, the integer output paramter ''info'' returns 0,. Otherwise, if parameter number ''pos'' is the first to have an illegal value, the return value of ''info'' is ''-pos'' and in addition an error message will be produced.

====  The Matrix Redistribution Routine PDGEMR2D  ====

A further possibility to initialize the local parts of a distributed matrix is provided within ScaLAPACK by the matrix redistribution subroutine ''PDGEMR2D''. This routine copies a submatrix from a given block-cyclically distributed global matrix A, characterized by a grid context ''ctxt_A'' and a matrix descriptor ''desc_A'' to a submatrix of the same size of another  block-cyclically distributed global matrix B characterized by ''ctxt_B'' and ''desc_B''. The two matrices A and B can have different sizes, block-sizes and can be distributed on different sized process grids defined in the set of participating processes. Of course the sizes of A and B must be large enough to contain the submatrices to be used as input and output for the copy operation. 

The calling sequence for the general case is

<code>
call PDGEMR2D(m, n, al, i_A, j_A, desc_A, bl, i_B, j_B, desc_B, ctxt)
m, n       number of rows and columns of the submatrix to be copied
al         local matrix of size llda * nla, llda = desc_A(9),
           nla is the number of columns of the global submatrix of A stored in the process,
           equal to NUMROC( n_A, nb_A, mycol, csrc_A , npcol ) 
           on input, al must contain the local elements of the submatrix of the global matrix A
i_A, j_A   indices in the global matrix A, where the submatrix to be copied starts
desc_A     descriptor of global matrix A
           for processes not belonging to ctxt_A:
           desc_A(2) must be set to -1, other values of desc_A are ignored
bl         local matrix of size lldb * nlb, lldb = desc_B(9),
           nlb is the number of columns of the global submatrix of B stored in the process,
           equal to NUMROC( n_B, nb_B, mycol, csrc_B , npcol ) 
           on output, bl will contain the local elements of the submatrix of the global matrix B
i_B, j_B   indices in the global matrix B, where the submatrix to be stored starts
desc_B     descriptor of global matrix B
           for processes not belonging to ctxt_B:
           desc_B(2) must be set to -1, other values of desc_B are ignored
ctxt       valid context handel pointing to a process grid, 
           which contains all processes used in ctxt_A and ctxt_B
</code>

 
For generating a block-cyclically distributed copy of a global matrix A, the global matrix is generated as an ordinary two-dimensional local array ''a_0'' on process 0, which then is copied using ''PDGEMR2D'' to the local arrays ''al'' in the processes of the two dimensional process grid. The following lines of code illustrate this:

<code>
      call BLACS_PINFO(myid,nproc)
      call BLACS_GET( 0, 0, ctxt_sys )  
! Set up a process grid of size 3*2 with context ctxt
      ctxt = ctxt_sys
      call BLACS_GRIDINIT( ctxt, 'R', 3, 2)
! Set up a process grid of size 1 with context ctxt_0
      ctxt_0 = ctxt_sys
      call BLACS_GRIDINIT( ctxt_0, 'R', 1, 1)
! Get the process coordinates in the grid context ctxt
      call BLACS_GRIDINFO( ctxt, nprow, npcol, myrow, mycol )
! initialize the descriptor for the distributed mxn matrix:
      llda =  max(1,NUMROC( m, mb, myrow, 0, nprow ))
      call DESCINIT( desca, m, n, mb, nb, 0, 0, ctxt, llda, info )
! initialize the descriptor desca_0 for the global mxn matrix:
      desca_0(2) = -1 
      if (myid.eq.0) 
     :     call DESCINIT( desca_0, m, n, m, n, 0, 0, ctxt_0, m,info )
!  on proc 0 initialize the global matrix a_0
      if (myid.eq.0) then
         do  j = 1, n
            do  i = 1 , m
               a_0(i + m * (j - 1)) = matval(i,j)
            end do
         end do
      end if
! copy a_0 on process 0 block-cyclically to the local arrays al in the process grid
      call PDGEMR2D(m, n, a_0, 1, 1, desca_0, al, 1, 1, desca, ctxt)
</code>

A complete example program for using ''PDGEMR2D'' can be copied from [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/use_PDGEMR2D.f|here]].

=====  Example Programs with PBLAS Routines  =====

PBLAS provides basic linear algebra operations on vectors and matrices for parallel execution in a message passing programming environment.
A matrix operand in PBLAS is a submatrix sub(A) with m rows and n columns of a global matrix A with m_A rows and n_A columns, which starts at the element of A with rowindex i_A and column_index j_A.

<code>
sub(A)(1:m,1:n) = A(i_A:i_A+m-1,j_A:j_A+n-1)
</code>
<div id="submatrix_restrictions"></div>
For a valid submatrix the following restrictions have to be met:

<code>
1<=i_A<=m_A, 1<=j_A<=n_A, i_A+m-1<=m_A, j_A+n-1 <=n_A
</code>

The local array al of a matrix operand in a call to a PBLAS routine is assumed to contain the elements of the global matrix A, which are mapped to the calling process according to the parameters in the descriptor ''desc_A'' describing the block-cyclic distribution of A in the given process grid. This mapping is independent of i_A and j_A defining the location of the global submatrix sub(A) within the global matrix A. Therefore the leading dimension  for the local array al, which has to be provided as ''desc_A(9)'' in the array descriptor for matrix A has to be larger or equal to ''NUMROC( m_A, mb_A, myrow, rsrc, nprow )'', independent of the number of rows m of sub(A) and of the location i_A of the first row of sub(A) within the global Matrix A. The number of columns in the local array al actually needed for the operation with the operand sub(A) is  ''NUMROC(j_A+n-1, nb_A, mycol, csrc, npcol )'' and thus depends on the size and location of sub(A).

A vector operand in PBLAS is a special submatrix of a global matrix X with only one column or one row:

<code>
subcol(X)(1:m) = X(i_X,i_X+m-1,j_X:j_X) 
subrow(X)(1:n) = X(i_X:i_X,j_X:j_X+n-1)
</code>

====  Parallel Matrix Vector Multiplication with PDGEMV  ====


The parallel matrix vector multiplication provided by the PBLAS routine ''PDGEMV'' works with three global matrices A, X, Y, which are distributed according to [[#array_descriptor|descriptors]] desc_A, desc_X, desc_Y. ''PDGEMV'' performs the operation

<code>
subvec(Y) = alpha sub(A) subvec(X) + beta subvec(Y)
</code>

where  ''alpha'' and ''beta'' are scalars, ''sub(A)'' is a submatrix of A, ''subvec(X),subvec(Y)'' are single subrows or subcolumns of X and Y. The detailed specification of the submatrices is dictated by the parameters ''trans,m,n,incx,incy'' in the calling sequence for ''PDGEMV''.

For ''trans='N'''>

<code>
              sub(A)(1:m,1:n) = A(i_A:i_A+m-1,j_A:j_A+n-1)
if (incx=1)   subcol(X)(1:n)  = X(i_X:i_X+n-1,j_X:j_X)
if (incx=m_X) subrow(X)(1:n)  = X(i_X:i_X,j_X:j_X+n-1)
if (incy=1)   subcol(Y)(1:m)  = Y(i_Y:i_Y+m-1,j_Y:j_Y)
if (incy=m_Y) subrow(Y)(1:m)  = Y(i_Y:i_Y,j_Y:j_Y+m-1)
</code>

For ''trans='T'''>

<code>
              sub(A)(1:n,1:m) = transpose of A(i_A:i_A+m-1,j_A:j_A+n-1)
if (incx=1)   subcol(X)(1:m)  = X(i_X:i_X+m-1,j_X:j_X)
if (incx=m_X) subrow(X)(1:m)  = X(i_X:i_X,j_X:j_X+m-1)
if (incy=1)   subcol(Y)(1:n)  = Y(i_Y:i_Y+n-1,j_Y:j_Y)
if (incy=m_Y) subrow(Y)(1:n)  = Y(i_Y:i_Y,j_Y:j_Y+n-1)
</code>

The calling sequence for ''PDGEMV'' is

<code>
 call PDGEMV(trans, m, n,
:            alpha, al, i_A, j_A, desc_A,
:                   xl, i_X, j_X, desc_X, incx
:            beta,  yl, i_Y, j_Y, desc_Y, incy )
</code>

The arguments have the following meaning:

<code>
trans: 
       if trans = 'N' the submatrix is taken directly from A  
       if trans = 'T' the submatrix is taken from the transposed of A
m : 
       if trans = 'N', number of rows in sub(A), number of elements in subvec(Y)
       if trans = 'T', number of columns in sub(A), number of elements in subvec(X)
n :    
       if trans = 'N', number of columns in sub(A), number of elements in subvec(X)
       if trans = 'T', number of rows in sub(A), number of elements in subvec(Y)
alpha:
       scalar of type real*8 multiplying sub(A) 
al:
       local array containing the local values of A 
i_A, j_A: 
       row and column number in A, where the submatrix sub(A) starts
desc_A:
       descriptor for the distribution of A
xl:
       local array containing the local values of X 
i_X, j_X: 
       row and column number in X where the subvector subvec(X) starts
desc_X:
       descriptor for the distribution of X
incx:
       if incx=1:    subvec(X) is a subcolumn of X
       if incx=m_X:  subvec(X) is a subrow of X
beta:
       scalar of type real*8 multiplying subvec(Y)
yl:
       local array containing the local values of Y 
i_Y, j_Y: 
       row and column number in Y where the subvector subvec(Y) starts
desc_Y:
       descriptor for the distribution of Y
incy:
       if incy=1:    subvec(Y) is a subcolumn of Y
       if incy=m_Y:  subvec(Y) is a subrow of Y
</code>

An example code for executing a matrix vector multiplication using ''PDGEMV'' with general values for sizes and distributions of global matrices and sizes and location of submatrices can be found [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/use_PDGEMV.f|here]]. Of course the parameters for sizes of the matrices and the positions of submatrices within the matrices have to be consistent [[#submatrix_restrictions|(submatrix_restrictions)]], otherwise the program will abort with an error message.

In the simplest case of the multiplication of an m x n matrix by an n vector yielding an m vector, the global matrix A has size m x n, the global matrix X has size n x 1, the global matrix Y has size m x 1, the respective  blocksizes are chosen as mb x nb, nb x 1, mb x 1 for A, X, Y resp. and the root process for all three global matrices A, X and Y has grid coordinates (0,0). The parameters for the sizes and distributions of global matrices and sizes and locations of submatrices in this case are:

<code>
m_A = m, n_A = n, mb_A = mb, nb_A = nb, i_A = 1, j_A = 1
m_X = n, n_X = 1, mb_X = nb, nb_X = 1 , i_X = 1, j_X = 1
m_Y = m, n_Y = 1, mb_Y = mb, nb_Y = 1 , i_Y = 1, j_Y = 1
incx = 1, incy = 1, rsrc = 0, csrc = 0
</code>

On a column major process grid of size nprow x npcol the matrix A will be distributed block-cyclically over all processes of the grid, whereas X and Y will be distributed over the processes of the first column in the grid with blocksizes nb resp. mb. The relevant parts of a code for this case are

<code>! Parameters for distributed global matrices A, X, Y
      m = 69; n = 77; mb=12; nb = 13
! Parameters for processor grid
      nprow=2; npcol=3
! set data for A and X; A(i,j) = mati*i + matj*j, X(i,1) = veci*i+vecj
      mati = 1; matj = 1; veci=1; vecj = 1
! Set up a process grid of size nprow*npcol
      ctxt = ctxt_sys; CALL BLACS_GRIDINIT( ctxt, 'C', nprow, npcol)
! Get the process coordinates in the grid
      CALL BLACS_GRIDINFO( ctxt, nprow, npcol, myrow, mycol )
! number of rows and columns of local parts al, xl, yl for A, X, Y
      m_al = NUMROC( m, mb, myrow, 0, nprow )
      n_al = NUMROC( n, nb, mycol, 0, npcol )
      m_xl = NUMROC( n, nb, myrow, 0, nprow ); n_xl = 1
      m_yl = m_al; n_yl = 1
! initializing descriptors for the distributed matrices A, X, Y:
      llda =  max(1,m_al); lldx =  max(1,m_xl); lldy = max(1,m_yl)
      call DESCINIT( desc_A, m, n, mb, nb, 0, 0, ctxt, llda, info )
      call DESCINIT( desc_X, n, 1, nb, 1, 0, 0, ctxt, lldx, info )
      call DESCINIT( desc_Y, m, 1, mb, 1, 0, 0, ctxt, lldy, info )
! initialize in parallel the local parts al of A
      do jl = 1 , n_al
         call ind_lc_gl(nb,npcol,0,jl,mycol,j)
         do il = 1 , m_al
            call ind_lc_gl(mb,nprow,0,il,myrow,i)
            al(il+m_al*(jl-1))= mati*i+matj*j
         end do
      end do
! initialize in parallel the local parts xl of X
      do il = 1 , m_xl
         call ind_lc_gl(nb,nprow,0,il,myrow,i)
            xl(il)= veci*i+vecj
      end do
! calculate matrix vector product
      call PDGEMV( 'N', m, n,
     :             1.d0, al, 1, 1, desc_A,
     :                   xl, 1, 1, desc_X, 1,
     :             0.d0, yl, 1, 1, desc_Y, 1 )
</code>

The complete program for the use of ''PDGEMV'' with this simplest case is [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/demo_PDGEMV.f|here]].


=====  Example Programs with ScaLAPACK Driver Routines  =====

====  Parallel Equation Solution with PDGESV  ====

The ScaLAPACK driver routine ''PDGESV'' is a parallel solver of systems of linear equations

<code>
a x = b
</code>

where ''a'' is the coefficient matrix of size n x n, ''b'' is a given right hand side vector of size n and ''x'' is the vector of size n which solves this system of equations. 

Similar to the PBLAS routines, ''PDGESV'' works with global block-cyclically distributed matrices, A and B, and their descriptors desc_A and desc_B. As input matrix  ''a'' with the equation coefficients a submatrix sub(A) af A is used

<code>
sub(A)(1:n,1:n) = A(i_A:i_A+n-1,j_A:j_A+n-1)
</code>

and several (nrhs) right hand side vectors are given as a submatrix sub(B) of B

<code>
sub(B)(1:n,1:nrhs) = B(i_B:i_B+n-1,j_B:j_B+nrhs-1)
</code>

The nrhs solutions of these equations are stored in sub(B), overwriting the input of the right hand sides. 

A pivot vector indicating the exchange of rows in the solution process is stored as a subvector of size n of a global block-cyclically distributed vector V of size m_A and block siza mb_A:

<code>
subvec(V)(1:n) = V(i_A:i_A+n-1)
</code>

The calling sequence for ''PDGESV'' is

<code>
      call PDGESV (n, nrhs, 
     :             al, i_A, j_A, desc_A, ipvtl, 
     :             bl, i_B, j_B, desc_B, info)
</code>

The arguments have the following meaning:

<code>
n:
       number of columns and rows of the global coefficient matrix sub(A) 
       and number of rows of the global matrix of right hand sides sub(B) 
nrhs:
       number of right hand sides, ie. number of columns in the global matrix sub(B)
al:
       local array containing the local values of A 
i_A, j_A: 
       row and column number in A, where the submatrix sub(A) starts
desc_A:
       descriptor for the distribution of A
ipvtl:
       local integer array of size llda+mb_a.
       This array contains the pivoting information. 
       If i is the index of the global row corresponding to the local row index il
       then ipvtl(il) =  index of the global row which was swapped with global row i.
bl:
       local array containing the local values of B, 
       on entry the values of bl corresponding to sub(B) contain 
       the local parts of the right hand sides,
       on return they contain the local parts of the solution vectors
i_B, j_B: 
       row and column number in B where the submatrix submat(B) starts
desc_B:
       descriptor for the distribution of B
info:
       on return, if info=0 the coeeficient matrix sub(A) is nonsingular
       if info > 0 , the coefficient matrix sub(A) is singular 
           and the solutions are not determined.
</code>

As in the case of the PBLAS routines ''PDGEMV'' and ''PDGEMM'', the parameters for sizes of the matrices and the positions of submatrices within the matrices have to be consistent [[#submatrix_restrictions|(submatrix_restrictions)]], otherwise the program will abort with an error message. For ''PDGESV'' there are further restrictions for the layout of the distributed matrices:

<code>
Blocksizes:
       mb_A = nb_A
       mb_B = nb_A
Positions of submatrices:
       i_A-1 modulo mb_A = 0
       j_A-1 modulo nb_A = 0
       i_B-1 modulo mb_B = 0
       the process row containing the row i_A of A
            must also contain the row i_B of B
</code>

If one of these conditions is violated, the program will abort with an error message.

An example code for solving a system of equations using ''PDGESV'' with general values for sizes and distributions of global matrices and sizes and positions of submatrices can be found [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/use_PDGESV.f|here]].

In the simplest case the global matrix A of size nxn is the coefficient matrix and the global matrix B of size nx1 contains one right hand side on entry and the solution on return, and the root process for A and B has grid coordinates (0,0). The parameters for the sizes and distributions of global matrices and sizes and positions of submatrices in this case are:

<code>
m_A = n, n_A = n, mb_A = nb, nb_A = nb, i_A = 1, j_A = 1
m_B = n, n_B = 1, mb_B = nb, nb_B = 1 , i_B = 1, j_B = 1
rsrc_A = 0, csrc_A = 0, rsrc_B = 0, csrc_B = 0
</code>

On a column major process grid of size nprow x npcol the matrix A will be distributed block-cyclically over all processes of the grid, whereas B will be distributed over the processes of the first column in the grid with blocksizes nb. The relevant parts of a code for this case are

<code>
! Parameters for the distributed global matrices
      n = 9; nb = 2
! Parameters for processor grid
      nprow=2; npcol=3
! Set up a process grid of size nprow*npcol
      call BLACS_GET( 0, 0, ctxt ); call BLACS_GRIDINIT( ctxt, 'C', nprow, npcol)
! Get the process coordinates in the grid
      call BLACS_GRIDINFO( ctxt, nprow, npcol, myrow, mycol )
! number of rows and columns of loacal parts al, bl for A, B
      m_al = NUMROC( n, nb, myrow, 0, nprow )
      n_al = NUMROC( n, nb, mycol, 0, npcol )
      m_bl = NUMROC( n, nb, myrow, 0, nprow )
      n_bl = NUMROC( 1, 1, mycol,0, npcol )
! initializing descriptors for the distributed matrices A, B:
      llda = max(1,m_al); lldb = max(1,m_bl)
      call DESCINIT( desc_A, n, n, nb, nb, 0, 0, ctxt, llda, info )
      call DESCINIT( desc_B, n, 1, nb, 1, 0, 0, ctxt, lldb, info )

! initialize in parallel the local parts of A in al and in al_s
      do jl = 1 , n_al
         j = INDXL2G(jl,nb,mycol,0,npcol)
         do il = 1 , m_al
            i = INDXL2G(il,nb,myrow,0,nprow)
            al(il+m_al*(jl-1))= 1.d0/(1.+5.*abs(i-j))
         end do
      end do
! initialize in parallel the local parts of B in bl
      do il = 1 , m_bl
         i = INDXL2G(il,nb,myrow,0,nprow)
         bl(il)= 1.*i + 1.
      end do

! solve the systems of equations
      call PDGESV( n, 1, al, 1, 1, desc_A, vl, bl, 1, 1, desc_B, info)
</code>

The complete program for the use of ''PDGESV'' with this simplest case is [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/demo_PDGESV.f|here]]

====  Parallel Eigenproblem Solution with PDSYEV  ====

The ScaLAPACK driver routine ''PDSYEV'' determines eigenvalues λ<sub>i</sub> and eigenvectors v<sub>i</sub> of a real symmetric nxn matrix a:

''av<sub>i</sub> = λ<sub>i</sub>v<sub>i</sub> ,    i=1,…,n''

The symmetric matrix a must be provided as a nxn submatrix sub(A) of a global block-cyclically distributed matrix A

<code>
sub(A)(1:n,1:n) = A(i_A:i_A+n-1,j_A:j_A+n-1)
</code>

the ortonormal eigenvectors  v<sub>i</sub> will be returned as columns of a nxn submatrix sub(Z)  of a global block-cyclically distributed matrix Z

<code>
sub(Z)(1:n,1:n) = Z(i_Z:i_Z+n-1,j_Z:j_Z+n-1)
</code>

and all n eigenvalues λ<sub>i</sub>  are returned in ascending order on every process in a vector w.

The complete calling sequence for ''PDSYEV'' is

<code>
      call PDSYEV (jobz, uplow, n, al, i_A, j_A, desc_A, w,
     :                             zl, i_Z, j_Z, desc_Z, 
     :                             work, lwork, info )
</code>

where the meaning of the arguments is

<code>
jobz:
       if jobz = 'N', only eigenvalues are computed
       if jobz = 'V', eigenvalues and eigenvectors are computed
uplo:
       if uplo = 'U', the upper triangular part of the symmetric matrix sub(A) is used
       if uplo = 'L', the lower triangular part of the symmetric matrix sub(A) is used
n:
       order of sub(A) used in the computation 
al:
       local array containing the local values of A, 
       on return, the values of al corresponding to the upper or lower triangular part of sub(A) 
       are overwritten, depending on the value of uplow
i_A, j_A: 
       row and column number in A, where the submatrix sub(A) starts
desc_A:
       descriptor for the distribution of A
w:
       global array of size >=n, containing on return all n eigenvalues in ascending order  
zl:
       local array for the local values of Z, 
       containing on return the local values of the eigenvectors, if jobz = 'V'
i_Z, j_Z: 
       row and column number in Z where the submatrix submat(Z) starts
desc_Z:
       descriptor for the distribution of Z
work:
       local workspace array of size ≥ max(1,lwork)
       if lwork = -1, work(1) returns the size of the needed workspace
lwork:
       size of workspace reserved locally. If lwork >0 on input is less than this size, 
          the program aborts.
       If on input globally lwork=-1, PDSYEV returns the needed workspace in work(1) 
          without solving the eigenproblem
info:
       on return,
       if info = 0, the execution is successful,
       if info < 0: 
          if the i-th argument is an array and the j-entry had an illegal value, 
          then info = -(i*100+j), 
          if the i-th argument is a scalar and had an illegal value, 
          then info = -i,
       if info= n+1, then PDSYEV has detected heterogeneity by finding 
          that eigenvalues were not identical across the process grid;
          in this case, the accuracy of the results from PDSYEV cannot be guaranteed
</code>

The [[http://www.netlib.org/scalapack/explore-html/d0/d1a/pdsyev_8f_source.html|source code]] for ''PDSYEV'' includes a warning concerning the validity of its results: "In its present form, PDSYEV assumes a homogeneous system and makes no checks for consistency of the eigenvalues or eigenvectors across the different processes.  Because of this, it is possible that a heterogeneous system may return incorrect results without any error messages." 
The reason for possible inconsitencies is, that some operations in ''PDGEMV'' are performed redundently in all participating processes. The parallel algorithm guaranties correctness only, if these redundant calculations give identical results on all processes. If the participating processes run on heterogeneous cpu's with different implementations for handling the rounding or overflow of floating point operations, the redundant operatioons might give different sesults, leding to inconsistecies. Since there is no complete testing for such inconsistencies, the results may be incorrect even if ''info=0'' on return. More information on the problem of heterogeneity in parallel algorithms can be fond in [[http://www.netlib.org/utk/people/JackDongarra/pdf/prac-het.pdf|this article by Blackford et al.]]. 


As in the case of the PBLAS routines ''PDGEMV'' and ''PDGEMM'', the parameters for sizes of the matrices and the positions of submatrices within the matrices have to be consistent [[#submatrix_restrictions|(submatrix_restrictions)]], otherwise the program will abort with an error message. For ''PDSYEV'' there are further restrictions for the layout of the distributed matrices:

<code>
Global matrices:
       m_A = m_Z
       n_A = n_Z
       mb_A = nb_A = mb_Z = nb_Z = nb
       rsrc_A = rsrc_Z
       csrc_A = csrc_Z
Positions of submatrices:
       i_A-1 modulo nb = 0
       j_A-1 modulo nb = 0
       i_Z-1 modulo nb = 0
       the process row containing the row i_A of A
            must also contain the row i_Z of Z
</code>

If one of these conditions is violated, the program will abort with an error message.

Furthermore, tests with ''PDSYEV'' have shown, that wrong results are produced without error message, if the offsets for the submatrices ''sub(A)'' and ''sub(Z)'' are different. This restriction is not described in the documentation contained in the [[http://www.netlib.org/scalapack/explore-html/d0/d1a/pdsyev_8f_source.html|source code]] for ''PDSYEV''. 

Each of the local processes has to provide memory space for the local parts of the global matrices A and Z and for all n eigenvalues. In addition memory space for the local workarray ''work'' is needed.  The amount of memory locally allocated for ''work'' can be prescribed with the input parameter ''lwork''.  If on input ''lwork>0'', it must be larger or equal than ''size_w'', the workspace needed locally in a process. Two main components contribute to ''size_w'', one is proportional to the global order  ''n'' of the eigenproblem to be solved, one is proportional to the number of elements of the problem matrix stored locally. 

If ''jobz='N''' (no eigenvectors are requested)

<code>
size_w =  n*5 + lel + 1
     lel = workspace requirement for ScaLAPACK routine PDSYTRD
         ≤ nb * max(LELM(n,nb,nprow)+1,3)
</code>

where [[#LELM|''LELM(n,nb,nprow)=nb⌈n/(nb*nprow)⌉'']] is the maximum of local elements over all processes in a column of the process grid. 
For blocksizes larger than 2 therefore the following value for ''lwork'' will be sufficient on all processes:

<code>
lwork = 5n + nb( ⌊(n-1)/nprow⌋ + nb) + 1
</code>


If ''jobz='V''' (eigenvalues and eigenvectors are requested)

<code>
size_w =  n*5 +  max(2*(n-2), lel1) + n*lel2 +  1
     lel1 = workspace requirement for ScaLAPACK routine PDORMTR when it's SIDE argument is 'L'
          ≤ nb*(LELM(n+1,nb,nprow)+LELM(n+1,nb,npcol) +nb)
     lel2 = number of rows of a block-cyclically (blocksize nb) 
             global matrix with n rows 
             stored locally in a process of a process grid of size (nprow*npcol) * 1
           ≤ LELM(n,nb,nprow*npcol)
</code>

A sufficient large value for  ''lwork'' in this case is (using ''A = ⌊n/(nb*nprow)⌋+⌊n/(nb*npcol)⌋'')

<code>
lwork = 5n + max(2n,(3+A)nb²) +nb(⌊(n-1)/(nb*nprow*npcol)⌋+1)n + 1
</code>

An example code for the determination of eigenvalues and eigenvectors using ''PDSYEV'' can be found [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/use_PDSYEV.f|here]].

A demo code for simple values for all parameters is [[http://wwwuser.gwdg.de/~ohaan/ScaLAPACK_examples/demo_PDSYEV.f|here]]. The relevant patrs for this code are
<code>
! input parameters
      n = 1000; nb = 24
! Parameters for processor grid
      nprow=2; npcol=2
! Set up a process grid of size nprow*npcol
      call BLACS_GET( 0, 0, ctxt ); call BLACS_GRIDINIT( ctxt, 'C', nprow, npcol)
! Get the process coordinates in the grid
      call BLACS_GRIDINFO( ctxt, nprow, npcol, myrow, mycol )
! number of rows and columns of local parts al, zl for A, Z
      m_al = NUMROC( n, nb, myrow, 0, nprow )
      n_al = NUMROC( n, nb, mycol, 0, npcol )
! size of workarray
      lel1 = nb**2*(n/(nb*nprow)+n/(nb*npcol)+3)
      lel2 = nb*((n-1)/(nb*nprow*npcol)+1)
      lwork = n*5 + max(2*n,lel1) + n*lel2 + 1
! initializing descriptors for the distributed matrices:
      llda = max(1,m_al)
      call DESCINIT( desc_A, n, n, nb, nb,
     :               0, 0, ctxt, llda, info )
! initialize in parallel the local parts of A in al 
      do jl = 1 , n_al
         j = INDXL2G(jl,nb,mycol,0,npcol)
         do il = 1 , m_al
            i = INDXL2G(il,nb,myrow,0,nprow)
            al(il+m_al*(jl-1))= 1./(1.+5.*abs(i-j))
         end do
      end do
! calculate eigenvalues and eigenvectors
      call PDSYEV( 'V', 'U', n, al, 1, 1, desc_A, w,
     :                          zl, 1, 1, desc_A,
     :                          work, lwork, info )
</code>
 
[[Kategorie: Scientific Computing]]