		program solvesystemMPI
		implicit none
!!!!		use mpi
                include "mpif.h"
		integer i, j
		integer ierr, nprocs, me
! BLACS/SCALAPACK
		integer  context, nprow, npcol, myrow, mycol, 
     %                                   info, ndims, dims(2)
		parameter ( ndims=2 )
		integer  numroc
		real*8 myrandom
	
		integer  M, N, MB, NB, Nloc, Mloc, aval
!		parameter ( M=8192, N=8192, MB=MYBLOCK, NB=MYBLOCK )
		parameter ( M=8192, N=8192)
                character*128 cnb
		real*8 , allocatable :: A(:,:), B(:), work(:)
		integer, allocatable :: ipiv(:)
		integer  descA(9), descB(9)
		real*8  t1, t2, dt
                integer iseed
                integer iac,iargc
                iac=iargc()
                call getarg(1,cnb)
!                read(cnb,'(I8)')NB
                NB=256
                MB=NB
	
! Initialize MPI environment
		call MPI_Init(ierr)
		call MPI_Comm_Size(mpi_comm_world,nprocs,ierr)
		call MPI_Comm_Rank(mpi_comm_world,me,ierr)
		call init_randoms_generator(iseed)	
! 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(*,*) "me= ",me,"     Mloc= ",Mloc,"     Nloc= ",Nloc 
		allocate( A( Mloc, Nloc ) )
		allocate( B( Mloc ) )
		allocate( ipiv( Mloc+MB ) )
		allocate( work( M ) )
	
		call descinit(descA,M,N,MB,NB,0,0,context,max(1,Mloc),info)
		call descinit(descB,M,1,MB, 1,0,0,context,max(1,Mloc),info)
	
		do j=1,N
		  do i=1,M
		      call pdelset( A, i, j, descA, myrandom(iseed) )
		  enddo
		enddo
	do i=1,M
	  call pdelset( B, i, 1, descB, myrandom(iseed) )
	enddo
	
	call MPI_Barrier( mpi_comm_world, ierr)
	t1 = MPI_wtime()
	
! Linear system equations solver
	call pdgesv( N, 1, A, 1, 1, descA, ipiv, B, 1, 1, descB, info )
	
!	call MPI_Barrier( mpi_comm_world, ierr)
	t2 = MPI_wtime()
	dt = t2-t1
	if(me.eq.0) write(*,*) 'with ', nprocs, ' procs, Walltime ',dt
	
! Print the results vector
!	call pdlaprnt(M, 1, B, 1, 1, descB, 0, 0, 'X', 6, work)
	
	deallocate( A )
	deallocate( B )
	deallocate( ipiv )
	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 
	
