MPI IO Czytanie i pisanie Matryca cykliczna

Question

Mam projekt szkolny do mnożenia macierzy w rozproszonym systemie hpc.

Muszę odczytywać w macierzy z równoległego systemu IO i używać pblacs do równoległego wykonywania mnożenia macierzy na wielu węzłach obliczeniowych (procesorach). Dane muszą być odczytywane przy użyciu poleceń MPI IO. Wiem, że PBlacs używa blokowych rozkładów cyklicznych do wykonania mnożenia.

Profesor nie przekazał nam zbyt wiele informacji na temat MPI IO i mam problem ze znalezieniem wielu informacji/zasobów na ten temat. Czy istnieją sposoby odczytywania w macierzy z systemu równoległego w cykliczny sposób i łatwego podłączania do pblacs pdgemm?

Wszelkie wskazówki dotyczące przydatnych zasobów byłyby mile widziane. Mam mało czasu i denerwuje mnie brak ukierunkowania na ten projekt.

Answer 1

Jest to względnie proste do zrobienia (jeśli już wiesz coś o blacs/scalapack i mpi-io!), Ale nawet wtedy dokumentacja - nawet online - jest taka, jaką odkryłeś, nieco słaba.

Pierwszą rzeczą, którą należy wiedzieć o MPI-IO, jest to, że pozwala używać zwykłych typów danych MPI do określania "" widoku "pliku, a następnie odczytywania tylko danych należących do tego widoku. W naszym centrum mamy slajdy i kod źródłowy na półdniowy kurs równoległego IO; pierwsza trzecia dotyczy podstaw MPI-IO. Dostępne są slajdy here i przykładowy kod źródłowy here.

Po drugie, MPI ma wbudowany sposób tworzenia typów danych "macierzy rozproszonych", z których jedna umożliwia rozkładanie cykliczne danych; omówiono to ogólnie w mojej odpowiedzi na to pytanie: What is the difference between darray and subarray in mpi?.

Oznacza to, że jeśli masz plik binarny zawierający dużą matrycę, możesz go odczytać za pomocą mpi-io przy użyciu MPI_Type_create_darray, a będzie on dystrybuowany przez zadania w sposób cykliczny. W takim razie jest to tylko kwestia wykonania połączenia blac lub skalapack. Przykładowy program użycia psapemv skalapack do mnożenia wektorów macierzowych zamiast psgemma znajduje się w mojej odpowiedzi na na tablicy wymiany danych obliczeniowych.

Aby dać ci pojęcie, w jaki sposób elementy pasują do siebie, poniżej znajduje się prosty program, który czyta w pliku binarnym zawierającym macierz (najpierw rozmiar kwadratowej macierzy N, a następnie elementy N^2), a następnie oblicza wartości własne i wektory za pomocą procedury scalapack (new) pssyevr. Łączy w sobie elementy MPI-IO, darray i scalapack. Jest w języku Fortran, ale wywołania funkcji są takie same w językach opartych na języku 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