From ea273e8e938a261d9a364642949880e1ead17c3e Mon Sep 17 00:00:00 2001 From: Stack-1 Date: Mon, 17 Aug 2026 16:12:19 +0200 Subject: [PATCH] [ADD] Binary matrix format with a parallel reader Matrix Market is ASCII with variable-length lines, so the byte offset of nonzero k cannot be computed without scanning the file: reading it in parallel is impossible without a redistribution pass. psb_mm2bin converts it once into a row-sorted binary layout, and psb_binmat_mod reads that layout in parallel. Layout, 8-byte native-endian throughout: a header carrying a magic number, the index and real kinds the file was written with, and the dimensions; then row_ptr(nrows+1), col_idx(nnz) and values(nnz). The kinds are recorded so a mismatched reader fails loudly instead of reinterpreting the bytes. Each rank reads the header, then only its own slice of row_ptr, and from that seeks straight to the nonzeros of the rows it owns. Because the block-row partition is implicit in the file order, psb_matdist is not needed at all -- that scatter alone accounted for 28 s of MPI time at 448 ranks. Geo_1438, 1.44M rows, 63.2M nonzeros: Matrix Market read, serial 27.0 s conversion, once 27.9 s binary read at 8 ranks, including assembly and one SpMV 0.98 s Verified: byte-for-byte layout check on a 4x4 symmetric matrix with the lower triangle expanded; A*1 identical at 1, 2, 3 and 4 ranks there, the 3-rank case covering the uneven split; sum(A*1) identical to sixteen digits at 1, 4 and 8 ranks on Geo_1438. Positioned stream reads rather than MPI-IO: the ranks touch disjoint byte ranges, so there is nothing to coordinate. If collective buffering turns out to matter at high rank counts, MPI_File_read_at_all is a change local to psb_binmat_mod. --- test/comm/psb_binmat_mod.f90 | 125 +++++++++++++++++++++++++++++++++++ test/comm/psb_mm2bin.f90 | 106 +++++++++++++++++++++++++++++ 2 files changed, 231 insertions(+) create mode 100644 test/comm/psb_binmat_mod.f90 create mode 100644 test/comm/psb_mm2bin.f90 diff --git a/test/comm/psb_binmat_mod.f90 b/test/comm/psb_binmat_mod.f90 new file mode 100644 index 000000000..4bb633bf3 --- /dev/null +++ b/test/comm/psb_binmat_mod.f90 @@ -0,0 +1,125 @@ +! +! Parallel read of the binary layout written by psb_mm2bin. +! +! Every rank opens the file, reads the header and only its own slice of row_ptr, +! then seeks straight to the nonzeros of the rows it owns. Nothing is read twice +! and nothing is redistributed: the block-row partition is implicit in the file +! order, so psb_matdist -- 28 s of MPI time at 448 ranks -- is not needed at all. +! +! Positioned stream reads are used rather than MPI-IO. The ranks touch disjoint +! byte ranges, so there is nothing to coordinate; if collective buffering ever +! turns out to matter at high rank counts, swapping these for +! MPI_File_read_at_all is a change local to this file. +! +module psb_binmat_mod + use psb_base_mod + implicit none + + integer(psb_lpk_), parameter :: psb_bin_magic = 6510534141251461121_psb_lpk_ + integer(psb_lpk_), parameter :: psb_bin_version = 1_psb_lpk_ + +contains + + subroutine psb_bin_mat_read(ctxt, fname, afmt, a, desc_a, info) + type(psb_ctxt_type), intent(in) :: ctxt + character(len=*), intent(in) :: fname, afmt + type(psb_dspmat_type), intent(out) :: a + type(psb_desc_type), intent(out) :: desc_a + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: iam, np, iu, nloc, ist + integer(psb_lpk_) :: magic, vers, ipkb, rpkb, nr, nc, nz, res + integer(psb_lpk_) :: r1, r2, k1, k2, lnz, i, k, base, cbase, vbase, nb, rem + integer(psb_lpk_), allocatable :: rp(:), ia(:), ja(:) + real(psb_dpk_), allocatable :: val(:) + + info = psb_success_ + call psb_info(ctxt, iam, np) + + iu = 71 + iam + open(unit=iu, file=trim(fname), access='stream', form='unformatted', & + & status='old', action='read', iostat=ist) + if (ist /= 0) then + if (iam == psb_root_) write(psb_err_unit,*) 'cannot open ', trim(fname) + info = psb_err_internal_error_ + return + end if + + read(iu, pos=1, iostat=ist) magic, vers, ipkb, rpkb, nr, nc, nz, res + if ((ist /= 0).or.(magic /= psb_bin_magic)) then + if (iam == psb_root_) write(psb_err_unit,*) 'not a psb binary matrix: ', trim(fname) + info = psb_err_internal_error_ ; close(iu) ; return + end if + ! The header records the kinds the file was written with: a mismatch has to + ! fail here rather than silently reinterpret the bytes. + if ((ipkb /= 8).or.(rpkb /= 8)) then + if (iam == psb_root_) write(psb_err_unit,*) 'kind mismatch in ', trim(fname), ipkb, rpkb + info = psb_err_internal_error_ ; close(iu) ; return + end if + if (nr /= nc) then + if (iam == psb_root_) write(psb_err_unit,*) 'matrix is not square' + info = psb_err_internal_error_ ; close(iu) ; return + end if + + ! Block-row split, in rank order: the same assignment psb_cdall(nl=) makes, + ! which is what lets the file order stand in for the partition. + nb = nr/np + rem = mod(nr,np) + if (iam < rem) then + nloc = int(nb+1, psb_ipk_) + r1 = iam*(nb+1) + 1 + else + nloc = int(nb, psb_ipk_) + r1 = rem*(nb+1) + (iam-rem)*nb + 1 + end if + r2 = r1 + nloc - 1 + + base = 65 ! row_ptr starts here (1-based stream pos) + cbase = base + (nr+1)*8 ! col_idx + vbase = cbase + nz*8 ! values + + allocate(rp(nloc+1), stat=ist) + if (ist /= 0) then + info = psb_err_alloc_dealloc_ ; close(iu) ; return + end if + read(iu, pos=base + (r1-1)*8, iostat=ist) rp + if (ist /= 0) then + info = psb_err_internal_error_ ; close(iu) ; return + end if + + k1 = rp(1) + k2 = rp(nloc+1) - 1 + lnz = k2 - k1 + 1 + if (lnz < 0) lnz = 0 + + allocate(ia(max(lnz,1_psb_lpk_)), ja(max(lnz,1_psb_lpk_)), & + & val(max(lnz,1_psb_lpk_)), stat=ist) + if (ist /= 0) then + info = psb_err_alloc_dealloc_ ; close(iu) ; return + end if + + if (lnz > 0) then + read(iu, pos=cbase + (k1-1)*8, iostat=ist) ja(1:lnz) + if (ist == 0) read(iu, pos=vbase + (k1-1)*8, iostat=ist) val(1:lnz) + if (ist /= 0) then + info = psb_err_internal_error_ ; close(iu) ; return + end if + ! Global row index for every nonzero, from the row_ptr slice. + do i = 1, nloc + do k = rp(i), rp(i+1)-1 + ia(k-k1+1) = r1 + i - 1 + end do + end do + end if + close(iu) + + call psb_cdall(ctxt, desc_a, info, nl=nloc) + if (info == psb_success_) call psb_spall(a, desc_a, info, nnz=int(lnz,psb_ipk_)) + if ((info == psb_success_).and.(lnz > 0)) & + & call psb_spins(int(lnz,psb_ipk_), ia, ja, val, a, desc_a, info) + if (info == psb_success_) call psb_cdasb(desc_a, info) + if (info == psb_success_) call psb_spasb(a, desc_a, info, afmt=afmt) + + end subroutine psb_bin_mat_read + +end module psb_binmat_mod diff --git a/test/comm/psb_mm2bin.f90 b/test/comm/psb_mm2bin.f90 new file mode 100644 index 000000000..1606af54a --- /dev/null +++ b/test/comm/psb_mm2bin.f90 @@ -0,0 +1,106 @@ +! +! Convert a Matrix Market file into a row-sorted binary CSR-like layout that can +! be read in parallel. +! +! Why this exists: Matrix Market is ASCII with variable-length lines, so the byte +! offset of nonzero k cannot be computed without scanning the file. Reading it in +! parallel is therefore impossible without a redistribution step. Converting once +! to a fixed-record binary layout lets every rank seek straight to its own rows, +! which also removes the need for psb_matdist -- at 448 ranks that scatter cost +! 28 s of MPI time on its own. +! +! Layout, all 8-byte native-endian: +! +! header magic, version, ipk_bytes, rpk_bytes, nrows, ncols, nnz, reserved +! row_ptr nrows+1 entries, 1-based, row i occupies row_ptr(i)..row_ptr(i+1)-1 +! col_idx nnz entries, 1-based global column indices +! values nnz doubles +! +! Native endianness on purpose: this is an intermediate for local runs, and the +! header records the kinds so a mismatched reader fails loudly instead of +! producing garbage. +! +program psb_mm2bin + use psb_base_mod + use psb_util_mod + implicit none + + integer(psb_lpk_), parameter :: bin_magic = 6510534141251461121_psb_lpk_ + integer(psb_lpk_), parameter :: bin_version = 1_psb_lpk_ + + type(psb_ldspmat_type) :: aux + type(psb_ld_coo_sparse_mat) :: coo + character(len=1024) :: infile, outfile + integer(psb_ipk_) :: info, iout + integer(psb_lpk_) :: nr, nc, nz, k, irow + integer(psb_lpk_), allocatable :: row_ptr(:) + real(psb_dpk_) :: t0, t1 + + if (command_argument_count() < 2) then + write(psb_err_unit,*) 'usage: psb_mm2bin ' + stop 1 + end if + call get_command_argument(1, infile) + call get_command_argument(2, outfile) + + call cpu_time(t0) + call mm_mat_read(aux, info, filename=trim(infile)) + if (info /= psb_success_) then + write(psb_err_unit,*) 'cannot read ', trim(infile), ' info=', info + stop 2 + end if + call cpu_time(t1) + write(psb_out_unit,'("read : ",f8.2," s")') t1-t0 + + nr = aux%get_nrows() + nc = aux%get_ncols() + + ! Move the storage into a COO object and sort it row-major: fix() is what + ! guarantees the entries of a row are contiguous, which is the whole premise + ! of the layout below. + call cpu_time(t0) + call aux%mv_to(coo) + call coo%fix(info) + if (info /= psb_success_) then + write(psb_err_unit,*) 'fix failed, info=', info + stop 3 + end if + nz = coo%get_nzeros() + call cpu_time(t1) + write(psb_out_unit,'("sort : ",f8.2," s")') t1-t0 + write(psb_out_unit,'("rows ",i0," cols ",i0," nonzeros ",i0)') nr, nc, nz + + ! row_ptr from the sorted row indices + allocate(row_ptr(nr+1), stat=info) + if (info /= 0) stop 4 + row_ptr = 0 + do k = 1, nz + irow = coo%ia(k) + row_ptr(irow+1) = row_ptr(irow+1) + 1 + end do + row_ptr(1) = 1 + do k = 1, nr + row_ptr(k+1) = row_ptr(k+1) + row_ptr(k) + end do + if (row_ptr(nr+1) /= nz+1) then + write(psb_err_unit,*) 'row_ptr inconsistent: ', row_ptr(nr+1), nz+1 + stop 5 + end if + + call cpu_time(t0) + iout = 91 + open(unit=iout, file=trim(outfile), access='stream', form='unformatted', & + & status='replace', action='write', iostat=info) + if (info /= 0) stop 6 + write(iout) bin_magic, bin_version, & + & int(storage_size(nz)/8, psb_lpk_), int(storage_size(coo%val(1))/8, psb_lpk_), & + & nr, nc, nz, 0_psb_lpk_ + write(iout) row_ptr + write(iout) coo%ja(1:nz) + write(iout) coo%val(1:nz) + close(iout) + call cpu_time(t1) + write(psb_out_unit,'("write : ",f8.2," s")') t1-t0 + write(psb_out_unit,'("output : ",a)') trim(outfile) + +end program psb_mm2bin