MPI IO Bloque de lectura y escritura Matriz cíclica

Question

Tengo un proyecto escolar para hacer la multiplicación de matrices en un sistema distribuido hpc.

Necesito leer en una matriz desde un sistema paralelo de IO y usar pblacs para realizar la multiplicación de matrices en paralelo en muchos nodos de procesamiento (procesadores). Los datos deben leerse utilizando comandos MPI IO. Sé que PBlacs usa distribuciones cíclicas de bloques para realizar la multiplicación.

El profesor no nos ha dado mucha información sobre MPI IO, y estoy teniendo problemas para encontrar mucha información/recursos sobre el mismo. Específicamente, ¿hay formas de leer en una matriz desde un sistema paralelo en un bloque de forma cíclica y fácilmente conectarlo a pblacs pdgemm?

Cualquier apuntador a los recursos útiles sería muy apreciado. Tengo poco tiempo y estoy frustrado por la falta de dirección en este proyecto.

Answer 1

Esto es relativamente sencillo de hacer (si usted ya sabe algo sobre blacs/scalapack y mpi-io!) Pero incluso entonces la documentación, incluso en línea, es como usted ha descubierto, algo pobre.

Lo primero que debe saber sobre MPI-IO es que le permite usar tipos de datos MPI normales para especificar la "vista" de cada proceso del archivo y luego solo leer los datos que entran en esa vista. En nuestro centro tenemos diapositivas y código fuente para un curso de medio día sobre IO paralelo; el primer tercio aproximadamente es sobre los conceptos básicos de MPI-IO. Hay diapositivas here y muestra el código fuente here.

Lo segundo que debe saber es que MPI tiene una forma integrada de crear tipos de datos de "matriz distribuida", una combinación de los cuales le permite diseñar una distribución de datos cíclica de bloques; eso se discute en términos generales en mi respuesta a esta pregunta: What is the difference between darray and subarray in mpi?.

Esto significa que si tiene un archivo binario que contiene una matriz grande, puede leerlo con mpi-io usando MPI_Type_create_darray y se distribuirá por tareas de forma cíclica. Entonces solo se trata de hacer la llamada de blacs o scalapack. Un programa de ejemplo del uso de scalapack psgemv para la multiplicación de matrices-vectores en lugar de psgemm figura en la lista en mi respuesta a question en el intercambio de pilas de Ciencias computacionales.

Para hacerte una idea de cómo encajan las piezas, el siguiente es un programa simple que lee en un archivo binario que contiene una matriz (primero el tamaño de la matriz cuadrada N y luego los elementos N^2) y luego calcula los autovalores y vectores utilizando la rutina pssyevr de scalapack (nueva). Combina las cosas de MPI-IO, darray y scalapack. Está en Fortran, pero las llamadas de función son las mismas en los lenguajes basados ​​en C.

! 
! Use MPI-IO to read a diagonal matrix distributed block-cyclically, 
! use Scalapack to calculate its eigenvalues, and compare 
! to expected results. 
! 
program darray 
     use mpi 
     implicit none 

     integer :: n, nb ! problem size and block size 
     integer :: myArows, myAcols ! size of local subset of global array 
     real :: p 
     real, dimension(:), allocatable :: myA, myZ 
     real, dimension(:), allocatable :: work 
     integer, dimension(:), allocatable :: iwork 
     real, dimension(:), allocatable :: eigenvals 
     real, dimension(:), allocatable :: expected 
     integer :: worksize, totwork, iworksize 

     integer, external :: numroc ! blacs routine 
     integer :: me, procs, icontxt, prow, pcol, myrow, mycol ! blacs data 
     integer :: info ! scalapack return value 
     integer, dimension(9) :: ides_a, ides_z ! scalapack array desc 
     integer :: clock 
     real :: calctime, iotime 

     character(len=128) :: filename 
     integer :: mpirank 
     integer :: ierr 
     integer, dimension(2) :: pdims, dims, distribs, dargs 
     integer :: infile 
     integer, dimension(MPI_STATUS_SIZE) :: mpistatus 
     integer :: darray 
     integer :: locsize, nelements 
     integer(kind=MPI_ADDRESS_KIND) :: lb, locextent 
     integer(kind=MPI_OFFSET_KIND) :: disp 
     integer :: nargs 
     integer :: m, nz 

! Initialize MPI (for MPI-IO) 

     call MPI_Init(ierr) 
     call MPI_Comm_size(MPI_COMM_WORLD,procs,ierr) 
     call MPI_Comm_rank(MPI_COMM_WORLD,mpirank,ierr) 

! May as well get the process grid from MPI_Dims_create 
     pdims = 0 
     call MPI_Dims_create(procs, 2, pdims, ierr) 
     prow = pdims(1) 
     pcol = pdims(2) 

! get filename 
     nargs = command_argument_count() 
     if (nargs /= 1) then 
      print *,'Usage: darray filename' 
      print *,'  Where filename = name of binary matrix file.' 
      call MPI_Abort(MPI_COMM_WORLD,1,ierr) 
     endif 
     call get_command_argument(1, filename) 

! find the size of the array we'll be using 

     call tick(clock) 
     call MPI_File_open(MPI_COMM_WORLD, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, infile, ierr) 
     call MPI_File_read_all(infile,n,1,MPI_INTEGER,mpistatus,ierr) 
     call MPI_File_close(infile,ierr) 

! create the darray that will read in the data. 

     nb = 64 
     if (nb > (N/prow)) nb = N/prow 
     if (nb > (N/pcol)) nb = N/pcol 
     dims = [n,n] 
     distribs = [MPI_DISTRIBUTE_CYCLIC, MPI_DISTRIBUTE_CYCLIC] 
     dargs = [nb, nb] 

     call MPI_Type_create_darray(procs, mpirank, 2, dims, distribs, dargs, & 
            pdims, MPI_ORDER_FORTRAN, MPI_REAL, darray, ierr) 
     call MPI_Type_commit(darray,ierr) 

     call MPI_Type_size(darray, locsize, ierr) 
     nelements = locsize/4 
     call MPI_Type_get_extent(darray, lb, locextent, ierr) 

! Initialize local arrays  

     allocate(myA(nelements)) 
     allocate(myZ(nelements)) 
     allocate(eigenvals(n), expected(n)) 

! read in the data 
     call MPI_File_open(MPI_COMM_WORLD, trim(filename), MPI_MODE_RDONLY, MPI_INFO_NULL, infile, ierr) 
     disp = 4 ! skip N = 4 bytes 
     call MPI_File_set_view(infile, disp, MPI_REAL, darray, "native", MPI_INFO_NULL, ierr) 
     call MPI_File_read_all(infile, myA, nelements, MPI_REAL, mpistatus, ierr) 
     call MPI_File_close(infile,ierr) 

     iotime = tock(clock) 
     if (mpirank == 0) then 
      print *,'I/O time  = ', iotime 
     endif 

! Initialize blacs processor grid 

     call tick(clock) 
     call blacs_pinfo (me,procs) 

     call blacs_get  (-1, 0, icontxt) 
     call blacs_gridinit(icontxt, 'R', prow, pcol) 
     call blacs_gridinfo(icontxt, prow, pcol, myrow, mycol) 

     myArows = numroc(n, nb, myrow, 0, prow) 
     myAcols = numroc(n, nb, mycol, 0, pcol) 

! Construct local arrays 
! Global structure: matrix A of n rows and n columns 

! Prepare array descriptors for ScaLAPACK 
     call descinit(ides_a, n, n, nb, nb, 0, 0, icontxt, myArows, info) 
     call descinit(ides_z, n, n, nb, nb, 0, 0, icontxt, myArows, info) 

! Call ScaLAPACK library routine 

     allocate(work(1), iwork(1)) 
     iwork(1) = -1 
     work(1) = -1. 
     call pssyevr('V', 'A', 'U', n, myA, 1, 1, ides_a, -1.e20, +1.e20, 1, n, & 
        m, nz, eigenvals, myZ, 1, 1, ides_z, work, -1,   & 
        iwork, -1, info) 
     worksize = int(work(1))/2*3 
     iworksize = iwork(1)/2*3 
     print *, 'Local workspace ', worksize 
     call MPI_Reduce(worksize, totwork, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) 
     if (mpirank == 0) print *, ' total work space ', totwork 
     call MPI_Reduce(iworksize, totwork, 1, MPI_INTEGER, MPI_SUM, 0, MPI_COMM_WORLD, ierr) 
     if (mpirank == 0) print *, ' total iwork space ', totwork 
     deallocate(work,iwork) 
     allocate(work(worksize),iwork(iworksize)) 
     call pssyev('N','U',n,myA,1,1,ides_a,eigenvals,myZ,1,1,ides_z,work,worksize,info) 
     if (info /= 0) then 
     print *, 'Error: info = ', info 
     else if (mpirank == 0) then 
     print *, 'Calculated ', m, ' eigenvalues and ', nz, ' eigenvectors.' 
     endif 

! Deallocate the local arrays 

     deallocate(myA, myZ) 
     deallocate(work, iwork) 

! End blacs for processors that are used 

     call blacs_gridexit(icontxt) 
     calctime = tock(clock) 

! calculated the expected eigenvalues for a particular test matrix 

     p = 3. + sqrt((4. * n - 3.) * (n - 1.)*3./(n+1.)) 
     expected(1) = p/(n*(5.-2.*n)) 
     expected(2) = 6./(p*(n + 1.)) 
     expected(3:n) = 1. 

! Print results 

     if (me == 0) then 
     if (info /= 0) then 
      print *, 'Error -- info = ', info 
     endif 
     print *,'Eigenvalues L_infty err = ', & 
      maxval(abs(eigenvals-expected)) 
     print *,'Compute time = ', calctime 
     endif 

     deallocate(eigenvals, expected) 

     call MPI_Finalize(ierr) 


contains 
    subroutine tick(t) 
     integer, intent(OUT) :: t 

     call system_clock(t) 
    end subroutine tick 

    ! returns time in seconds from now to time described by t 
    real function tock(t) 
     integer, intent(in) :: t 
     integer :: now, clock_rate 

     call system_clock(now,clock_rate) 

     tock = real(now - t)/real(clock_rate) 
    end function tock 

end program darray