HPC 分散システムで行列乗算を行う学校のプロジェクトがあります。
並列 IO システムから行列を読み込み、pblacs を使用して多くの計算ノード (プロセッサ) で並列に行列乗算を実行する必要があります。MPI IO コマンドを使用してデータを読み込む必要があります。PBlacs がブロック循環分布を使用して乗算を実行していることは知っています。
教授は MPI IO についてあまり情報を提供してくれませんでした。私はそれに関する多くの情報/リソースを見つけるのに苦労しています。具体的には、ブロック循環方式で並列 io システムから行列を読み取り、それを pblacs pdgemm に簡単にプラグインする方法はありますか?
有用なリソースへのポインタは大歓迎です。私は少し時間が足りず、このプロジェクトの方向性の欠如に不満を感じています.
これは実際には比較的簡単に実行できます (blacs/scalapack と mpi-io について既に知っている場合)。
MPI-IO について最初に知っておくべきことは、通常の MPI データ型を使用して各プロセスのファイルの「ビュー」を指定し、そのビューに該当するデータのみを読み取ることができるということです。私たちのセンターには、並列 IO に関する半日コースのスライドとソース コードがあります。最初の 3 分の 1 程度は、MPI-IO の基本に関するものです。ここにスライドがあり、ここにサンプル ソース コードがあります。
2 番目に知っておくべきことは、MPI には「分散配列」データ型を作成する組み込みの方法があり、その組み合わせの 1 つでブロック循環データ分散をレイアウトできることです。これは、この質問に対する私の回答で一般的な用語で説明されています: mpi の darray と subarray の違いは何ですか? .
つまり、大きな行列を含むバイナリ ファイルがある場合は、mpi-io を使用して読み込むことができMPI_Type_create_darray、タスクによってブロック循環方式で分散されます。あとは、blacs または scalapack 呼び出しを実行するだけです。行列とベクトルの乗算に psgemm ではなく scalapack psgemv を使用するサンプル プログラムは、Computational Science スタック交換に関する質問への回答に記載されています。
ピースがどのように組み合わされるかを理解するために、次の単純なプログラムを使用して、行列 (最初に正方行列 N のサイズ、次に N^2 要素) を含むバイナリ ファイルを読み取り、固有値とscalapack の (新しい)pssyevrルーチンを使用したベクトル。MPI-IO、darray、および scalapack を組み合わせたものです。これは Fortran ですが、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