module processor_grid 

  implicit none 
 
  integer :: nproc, mype 
  integer :: rowid, colid 
  integer :: idest, isour 
  integer :: nprow, npcol 
  integer :: icontext 
 
contains 
 
  subroutine processor_grid_init( ) 
    integer :: ierr 
    integer :: ndims, dims(2), coor(2) 
    integer :: itest 
    ! 
    !     Initialize BLACS, get processor ID (mype) 
    !     and the number of processors (nproc) 
    !     NOTE: mype = 0, 1, .... nproc - 1 
    ! 
    CALL BLACS_PINFO( mype, nproc ) 
 
    !     find out nprow, npcol so that: 
    !     nprow * npcol = nproc 
 
    CALL gridsetup( nproc, nprow, npcol ) 
 
    !     Get a default BLACS context 
    ! 
    CALL BLACS_GET( -1, 0, icontext ) 
 
    !     Initialize a default BLACS context 
    ! 
    CALL BLACS_GRIDINIT( icontext, 'R', nprow, npcol ) 
 
    !     Get processor coordinates = rowid, colid 
    ! 
    CALL BLACS_GRIDINFO( icontext, nprow, npcol, rowid, colid ) 
 
    WRITE(6,*) 'mype = ', mype, ' row = ', rowid, ' col = ', colid 
    ! 
  end subroutine 
 
  subroutine processor_grid_finalize 
    ! 
    !     Close-up Blacs context 
    ! 
    CALL BLACS_GRIDEXIT( icontext ) 
    ! 
  end subroutine 
 
end module 
 
! 
! 
! 
! 
 
PROGRAM scalapack_example 
  ! 
  USE processor_grid 
  ! 
  implicit none 
  ! 
  integer :: n, nb 
  ! 
  CALL processor_grid_init() 
  ! 
  n  = 2000  !  matrix size 
  nb = 120   !  block size 
  ! 
  CALL sample_pdsyev_call( n, nb ) 
  ! 
  CALL processor_grid_finalize() 
  ! 
END PROGRAM 
 
! 
! 
! 
! 
 
SUBROUTINE SAMPLE_PDSYEV_CALL( n, nb ) 
 
   USE processor_grid 
 
   IMPLICIT NONE 
 
   integer             :: n, nb 
 
   REAL*8              :: t1, t2 
   REAL*8              :: cclock 
   EXTERNAL            :: cclock 
 
   INTEGER             :: lwork 
   INTEGER             :: i, j, ii, jj, info, ml, nl 
   INTEGER             :: descb( 50 ), desca( 50 ), descz( 50 ), idum 
   REAL*8, ALLOCATABLE :: a( :, : ), w( : ), z( :, : ), b( :, : ) 
   REAL*8, ALLOCATABLE :: btmp( :, : ) 
   REAL*8, ALLOCATABLE :: work( : ) 
   REAL*8              :: workp( n ) 
   REAL*8              :: wsize, aval 
   REAL*8              :: rtmp( nproc ) 
   ! 
   INTEGER             :: numroc, INDXL2G 
   EXTERNAL            :: numroc, INDXL2G 
 

      CALL BLACS_BARRIER (icontext, 'A') 
      if ( mype == 0 ) then 
        write (6, fmt="(' Matrix and Block size ',2I5)") N, NB 
      end if 
 
      !     Calculate the LOCAL dimensions of A 
      ! 
      nl = numroc( n, nb, rowid, 0, nprow ) 
      ml  = numroc( n, nb, colid, 0, npcol ) 
 
      ALLOCATE( a( nl, ml), w( n ), z(nl, ml) )  
      ! 
      !     These are basic array descriptors 
      ! 
      CALL DESCINIT( desca, n, n, nb, nb, 0, 0, icontext, nl, info ) 
      CALL DESCINIT( descz, n, n, nb, nb, 0, 0, icontext, nl, info ) 
      ! 
      !     Build a matrix that you can create with 
      !     a one line matlab command:  hilb(n) + diag([1:-1/n:1/n]) 
      ! 
      DO j = 1, n 
         DO i = 1, n 
            IF( i == j ) THEN 
               aval = ( DBLE( n-i+1 ) ) / DBLE( n ) + 1.0d0 / ( DBLE( i+j ) - 1.0d0 ) 
            ELSE 
               aval = 1.0d0 / ( DBLE( i+j ) - 1.0d0 ) 
            END IF 
            CALL PDELSET( a, i, j, desca, aval ) 
         END DO 
      END DO 
 
      if( n <= 16 ) THEN 
         !  
         !   print matrix on standard output 
         ! 
         CALL PDLAPRNT( n, n, a, 1, 1, desca, 0, 0, 'A', 6, workp )  
         ! 
      END IF 
 
      CALL BLACS_BARRIER (icontext, 'A') 
      ! 
      T1 = cclock() 
      ! 
      !  diagonalization driver subroutine 
      ! 
      lwork = -1 
      CALL PDSYEV( 'V', 'U', n, a, 1, 1, desca, w, z, 1, 1, descz, wsize, lwork, info ) 
 
      lwork = INT(wsize) + 1 
      allocate( work( lwork ) ) 
      CALL PDSYEV( 'V', 'U', n, a, 1, 1, desca, w, z, 1, 1, descz, work, lwork, info ) 
 
      deallocate( work ) 
 
      IF( info /= 0 ) THEN 
        PRINT *, ' *** DIAGONALIZATION FAILED ! *** '  
      END IF 
 
      T2 = cclock() 
 
      rtmp = 0.0d0 
      rtmp( mype + 1 ) = t2 - t1 
      ! 
      IF( mype == 0 ) THEN 
        WRITE( 6, * ) 'Seconds for PDSYEV = ', rtmp( mype+1 ) 
      END IF 
 
      if( n <= 16 ) THEN 
        CALL PDLAPRNT( n, n, z, 1, 1, descz, 0, 0, 'Z', 6, workp ) 
      END IF 
 
      ! 
      ! Inserire nel programma scal_esercizio.f90 codice necessario a verificare  
      ! l'ortogonalita' degli autovettori memorizzati nella  
      ! matrice z. 
      ! 
      ! traccia: 
      ! 1) allocare una matrice B opportunamente dimensionata 
      ! 2) definire il descrittore per B 
      ! 3) utilizzare la moltiplicazione di matrici parallela  
      !    SCALAPACK PDGEMM per calcolare B = Z^t Z 
      ! 4) i processori verificano che il proprio sottoblocco di 
      !    B corrisponde a quello della matrice identita',  
      !    utilizzando le funzioni SCALAPACK NUMROC e INDXL2G 
      ! 5) utilizzando cclock monitorare il tempo di esecuzione 
      !    di PDGEMM 
      !     
      ! 
      ALLOCATE( b( nl, ml ) ) 

      CALL DESCINIT( descb, n, n, nb, nb, 0, 0, icontext, nl, info ) 
 
      CALL BLACS_BARRIER ( icontext, 'A') 
      T1 = cclock() 
 
      CALL PDGEMM('t', 'n', n, n, n, 1.0d0, z(1,1), 1, 1, descz, z(1,1), 1, 1, descz, 0.0d0, b(1,1), 1, 1, descb) 
 
      T2 = cclock() 
      rtmp = 0.0d0 
      rtmp( mype + 1 ) = t2 - t1 
 
      IF( mype == 0 ) THEN 
        WRITE( 6, * ) 'GFlop/sec for PDGEMM = ', DBLE(n)**3 / rtmp( mype+1 ) * 1.d-9 
      END IF 
 
      DO j = 1, ml 
         DO i = 1, nl 
            ii = INDXL2G( i, nb, rowid, 0, nprow ) 
            jj = INDXL2G( j, nb, colid, 0, npcol ) 
            IF( ii == jj ) THEN 
               IF( ABS( b( i, j ) - 1.0d0 ) >= 1.0d-8 ) WRITE( 6, * ) 'eigenvector not orthogonal' 
            ELSE 
               IF( ABS( b( i, j ) ) >= 1.0d-8 ) WRITE( 6, * ) 'eigenvector not orthogonal' 
            END IF 
         END DO 
      END DO 
 
      IF( n <= 16 ) THEN 
        CALL PDLAPRNT( n, n, b, 1, 1, descb, 0, 0, 'B', 6, workp ) 
      ELSE 
        WRITE(6,*) ' SUM of B on pe ', mype, ' = ', SUM( b ) 
      END IF 
      ! 
      !     Print out diagonalization time and eigenvalues  
      ! 
      IF( mype == 0 ) THEN 
         DO i = 1, MIN(16,n) 
           PRINT *, ' W(', I, ')=', w( i ), ';' 
         END DO 
      END IF 
      ! 
      DEALLOCATE( a, w, z, b ) 
      ! 
   RETURN 
END SUBROUTINE 
 
! 
! 
! 
 
SUBROUTINE gridsetup( nproc, nprow, npcol ) 
 
! 
! This subroutine factorizes the number of processors (NPROC) 
! into NPROW and NPCOL,  that are the sizes of the 2D processors mesh. 
! 
! Written by Carlo Cavazzoni 
! 
 
      IMPLICIT NONE 
 
      integer nproc,nprow,npcol 
 
      integer sqrtnp,i 
 
      sqrtnp = int( sqrt( dble(nproc) ) + 1 ) 
      do i=1,sqrtnp 
        if(mod(nproc,i).eq.0) nprow = i 
      end do 
 
      npcol = nproc/nprow 
 
      return 
END SUBROUTINE 
