BLACS solution

C

#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#include "math.h"
#include "wrapper.h"
#define N 10
#define NB 2

#ifndef max
#define max( a, b ) ( ((a) > (b)) ? (a) : (b) )
#endif

int main( int argc, char** argv)
{
  int i, j, zero = 0, uno = 1, sei = 6;
  // MPI
  int myrank, nprocs, ndims = 2, dims[2] = {0,0};
  MPI_Status status;
  // BLACS/SCALAPACK
  int nprow, npcol, info, ic = -1, context, myrow, mycol, desca[9];
  char order = 'R';
  // MATRIX
  int n = N, nb = NB, nrhs = 1, nr, nc, lld;
  double *A, *B, *work, alpha;
  int *ipiv;

  // Initialize MPI environment
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

  // Initialize a default BLACS context and the processes grid
  MPI_Dims_create(nprocs, ndims, dims);
  nprow = dims[0];
  npcol = dims[1];
  Cblacs_get( ic, zero, &context );
  Cblacs_gridinit( &context, &order, nprow, npcol);
  Cblacs_gridinfo( context, &nprow, &npcol, &myrow, &mycol);

  // Computation of local matrix size
  nr = numroc_( &n, &nb, &myrow, &zero, &nprow );
  nc = numroc_( &n, &nb, &mycol, &zero, &npcol );
  lld = max( 1 , nr );
  A = malloc(nr*nc*sizeof(double));
  B = malloc(nr*nc*sizeof(double));
  work = malloc(nb*sizeof(double));

  // Descriptors
  descinit_( desca, &n, &n, &nb, &nb, &zero, &zero, &context, &lld, &info );

  // Matrix elements generation
  alpha = (double)(0);
  for( i = 1; i <= n; i++ ){
    for( j = 1; j <= n; j++ ){
      alpha += (double)(1);
      pdelset_( A, &i, &j, desca, &alpha );
    }
  }

  // Print matrix to file
  Cpdlawrite( "matrix.dat", &n, &n, A, &uno, &uno, desca, &zero, &zero, work );

  // Read matrix from file
  Cpdlaread( "matrix.dat", B, desca, &zero, &zero, work );

  // Print matrix to stdout
  Cpdlaprnt( &n, &n, B, &uno, &uno, desca, &zero, &zero, "A", &sei, work );

  free(A);
  free(B);
  free(work);

  // Close BLACS environment
  Cblacs_gridexit( context );
  Cblacs_exit( zero );

  // Close MPI environment if blacs_exit paramater is not equal zero
  //MPI_Finalize();
  return 0;
}
 

FORTRAN

program es_blacs
use mpi
implicit none

integer :: i, j, stdout
parameter ( stdout=6 )
! MPI
integer :: ierr, nprocs, myrank
! BLACS/SCALAPACK
integer :: context, nprow, npcol, myrow, mycol, info, ndims, dims(2)
parameter ( ndims=2 )
external :: blacs_exit, blacs_gridexit, blacs_gridinfo, blacs_get
external :: blacs_gridinit, descinit, numroc, pdelset
external :: pdlaprnt, pdlaread, pdlawrite
integer :: numroc
! Matrix
integer :: M, N, MB, NB, Nloc, Mloc, aval
parameter ( M=10, N=10, MB=2, NB=2 )
real*8 , allocatable :: A(:,:), B(:,:), work(:)
real*8 :: alpha
integer :: descA(9)

! Initialize MPI environment
call MPI_Init(ierr)
call MPI_Comm_Size(mpi_comm_world,nprocs,ierr)
call MPI_Comm_Rank(mpi_comm_world,myrank,ierr)

! Initialize a default BLACS context and the processes grid
dims = 0
call MPI_Dims_Create( nprocs, ndims, dims, ierr)
nprow = dims(1)
npcol = dims(2)
call blacs_get( -1, 0, context )
call blacs_gridinit( context, 'Row-major', nprow, npcol)
call blacs_gridinfo( context, nprow, npcol, myrow, mycol )

! Computation of local matrix size
Mloc = numroc( M, MB, myrow, 0, nprow )
Nloc = numroc( N, NB, mycol, 0, npcol )
write(*,*) "myrank= ",myrank,"     Mloc= ",Mloc,"     Nloc= ",Nloc 
allocate( A( Mloc, Nloc ) )
allocate( B( Mloc, Nloc ) )
allocate( work( M ) )

! Descriptor
call descinit(descA, M, N, MB, NB, 0, 0, context, max(1,Mloc), info)

! Matrix elements generation
do j=1,N
  do i=1,M
    alpha = dble((i-1)*N+j)
    call pdelset( A, i, j, descA, alpha )
  enddo
enddo

! Print the matrix to file
call pdlawrite('matrix.dat', M, N, A, 1, 1, descA, 0, 0, work)

! Print the matrix to file
call pdlaread('matrix.dat', B, descA, 0, 0, work)

! Print the matrix to stdout
call pdlaprnt(M, N, B, 1, 1, descA, 0, 0, 'A', stdout, work)

deallocate( A )
deallocate( B )
deallocate( work )

! Close BLACS environment
call blacs_gridexit( context )
call blacs_exit( 0 )

! Close MPI environment if blacs_exit paramater is not equal zero
!call MPI_Finalize(ierr)

end program es_blacs