这实际上是相对简单的（如果你已经知道一些关于 blacs/scalapack 和 mpi-io 的知识！）但即便如此，文档 - 甚至在线 - 正如你所发现的那样，有些差。
首先要了解 MPI-IO 是它允许您使用普通的 MPI 数据类型来指定文件的每个进程的“视图”，然后只读取属于该视图的数据。在我们的中心，我们有关于并行 IO 的半天课程的幻灯片和源代码；前三分之一左右是关于 MPI-IO 的基础知识。这里有幻灯片，这里有示例源代码。
第二件要知道的是，MPI 有一种内置的方式来创建“分布式数组”数据类型，其中的一种组合可以让您布局块循环数据分布；这在我对这个问题的回答中进行了一般性讨论： mpi 中的 darray 和 subarray 有什么区别？.
所以这意味着如果你有一个包含大矩阵的二进制文件，你可以使用 mpi-io 读取它，MPI_Type_create_darray它会以块循环的方式由任务分发。然后只需进行 blacs 或 scalapack 调用即可。我对计算科学堆栈交换问题的回答中列出了使用 scalapack psgemv 进行矩阵向量乘法而不是 psgemm 的示例程序。
为了让您了解这些部分是如何组合在一起的，下面是一个简单的程序，它读取包含矩阵的二进制文件（首先是方阵 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