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.
load_external_matrix was called by every rank and had no rank guard, so all of
them read the whole file. psb_matdist expects the global matrix on the root and
scatters from there, which is the pattern the fileread samples use.
For a 60M nonzero matrix that is over 1 GB of COO per rank: at 448 ranks over
4 nodes, more than 100 GB on each node, and the job was OOM-killed before it
reached the SpMV. Geo_1438, Flan_1565 and com-Orkut all died this way.
The global right-hand side and solution arrays are allocated full size on the
root only for the same reason: psb_scatter reads them nowhere else.
Paired A/B at 448 ranks, three replicas per arm inside one allocation: with
no_locks=true the aggregated MPI_Win_create cost 1463 s against 1450 s
without it, the two ranges fully overlapping. A probe with MPI_Win_get_info
confirmed this Open MPI retains the hint, so the null result belongs to the
hint and not to a failed measurement; same_disp_unit was silently ignored.
Removes the cached MPI_Info handles and their helpers, and returns the four
PSCW win_create calls to MPI_INFO_NULL.
Epoch handling. Both one-sided schemes opened and closed their access epoch
on every halo swap: push with lock_all/unlock_all over the window's whole
communicator, which is O(P) and not O(neighbours), pull with a
Win_lock/Win_unlock pair around every single Get. The epoch is now opened
once at window creation and closed by the handle teardown, which already
unlocked whenever window_open was set. Completion is explicit instead:
flush_all on the push side before the notifications go out, flush_local_all
on the pull side since a Get only needs to have landed in our own combuf, and
Win_sync before the scatter to reconcile the window's public and private
copies -- combuf is itself the exposed memory, so ordinary loads on it need
what the per-swap unlock used to provide implicitly.
448 ranks, dim 200, one scheme per process, 500 iterations, 3 replicas
rma_push 6.48e-2 -> 3.46e-4 (188x)
rma_pull 6.09e-4 -> 2.29e-4 (2.7x)
isend_irecv unchanged at 1.9e-4 (control)
rma_pull is now level with the point-to-point baseline while issuing no
two-sided traffic at all. rma_push stays at 1.8x because its completion
notification is still a point-to-point exchange.
The push loop is also split so the Puts are issued before any completion: a
flush per neighbour inside the loop serialised one round trip each. Measured
neutral on its own, kept because it removes one MPI call per neighbour.
Bound checks. y_nrows was assigned inside the do_send branch but read from
do_recv, so with split start/wait phases the second entry compared halo
indices against uninitialised stack memory: at 8 ranks and dim 200 that read
~32600 against legitimate indices up to 1040000 and aborted the run, while at
dim 40 the same garbage exceeded the halo and the check passed. The
minval/maxval index scans now run under the existing debug guard; they walked
the whole index list per neighbour per swap. The O(1) checks stay.
Dynamic windows (2a767a40) and the window pool (a3249c50) are kept in
history as measured negative results: no gain on window creation,
MPI_Get 49x slower, and a remote-access crash at 448 ranks. Back to
mpi_win_create, now with cached MPI_Info handles: same_disp_unit=true
always, no_locks=true only on the PSCW path in psi_dswapdata, the one
file where no passive-target call survives.
Validated at 8 ranks: psb_d_comm_test 9/9, and the AMG driver converges
in the same iterations to the same residual on all five schemes.
The four drivers the template generates in each of test/comm/{swapdata,
cg,spmv} were copied into the tree but appeared in no make rule, so they
were never compiled and never run. Only the generic double-precision
driver was built. That is why the complex halo exchange could stay broken
without anyone noticing.
Each Makefile now builds the four typed executables alongside the generic
one, which is kept: psb_spmv_kernel and psb_comm_cg_test are referenced by
the sbatch scripts. clean covers the new objects and executables, and
there are run/run-typed targets with NP and IDIM as before.
In cg and spmv the typed drivers are built CPU-only, without
-DPSB_HAVE_CUDA and without libpsb_cuda. The CUDA layer has no complex
support upstream: psb_?_cuda_vect_mod declares an external cdot/zdot,
which are not BLAS symbols (the complex dots are cdotc/cdotu), and
cuda/Makefile never compiles the c/z hdiag modules. Both are present in
development too, so this is not something the branch introduced. Building
all four the same way also keeps the cross-type comparison honest.
For the two types the CUDA layer does support, 'make gpu' builds
psb_{s,d}_*_gpu variants from the same sources with the CUDA defines and
libraries, into separate objects so switching between the two does not
force a full rebuild; 'make both' does both. The spmv drivers declare a
module, so the GPU build writes its module files to gpumod/ to avoid the
two variants overwriting each other.
The complex halo exchange declared real MPI datatypes on complex buffers,
with counts in elements and never doubled, so it moved the wrong number
of bytes. On z every scheme transferred half the payload; on c the
baseline and persistent paths happened to work because complex(spk) and
real(dpk) are both 8 bytes, while the neighbour collective and the RMA
get/put moved half.
Regenerated from template-psblas after fixing @RMPI_TYPE@ to @MPI_TYPE@
there; the only change in these files is the datatype handle, no other
line differs.
Only c and z are regenerated on purpose. psi_sswapdata.F90 carries the
index bounds checks and psi_dswapdata.F90 the PSCW optimization, neither
of which has been folded back into the template yet, so regenerating
those two would silently drop them. Neither needs this fix anyway: for
s and d the real and complex macros resolve to the same handle.
Per-type halo test on 4 ranks, all five schemes:
before s 9/9 d 9/9 c 3/9 z 2/9
after s 9/9 d 9/9 c 9/9 z 9/9
Confirmed end to end by the CG driver, where all five schemes now reach
the same final error to the last digit in the same 62 iterations.
GitHub's inline PDF viewer failed with 'Error loading PDF page number 1'.
Linearize the PDF (fast web view) so a streaming viewer can load page 1
before fetching the whole file. As a side effect qpdf recompresses the
streams (2.8M -> 0.68M). Fonts stay embedded (26/27; the remaining one is
a base-14 standard font), 188 pages, PDF 1.5.
- communication.md: replace the Mermaid state diagram with a plain-text
diagram (Mermaid showed 'Error rendering embedded code' on GitHub)
- regenerate psblas-3.9.pdf as PDF 1.5 (matching the previous build)
instead of TeX Live 2026's default 1.7, which GitHub's viewer failed
to load
Introduce docs/internals/ as a two-layer complement to the user manual:
- README.md: index plus the user-manual vs developer-guide rationale
- communication.md: the psb_halo/psb_ovrl -> psi_swapdata dispatch path,
the psb_comm_handle_type hierarchy and factory, the five MPI schemes
(isend/irecv, neighbor alltoallv, persistent, RMA pull/push), and the
start/wait/sync swap-status state machine
- psb_halo: document tran and mode arguments; note work is array-only
- psb_ovrl: document mode argument
- add psb_comm_status_{sync,start,wait}_ to the Named Constants list
- fix psb_add_ -> psb_sum_ in the psb_ovrl update operator list
- remove duplicated psb_gather synopsis line
- drop spurious alpha from the psb_halo data-type table header
- regenerate psblas-3.9.pdf
Rewrite psb_d_nest_cg_test to build the operator through the psb_d_nest_matrix
utility (init/ins/asb + get_owned_rows) instead of the low-level path, so no
per-field descriptor or l2g idiom appears in user-facing test code; x_exact=1
is set with x%set(done) rather than an l2g loop.
With this change psb_d_nest_cg_test fully subsumes psb_d_nest_builder_test
(same operator via the same builder, NONE plus DIAG/BJAC), so the latter is
removed. The test suite is now glob (square matvec), rect (rectangular
matvec) and cg (builder + preconditioned CG). Build hooks and README updated.
Author: Simone Staccone (Stack-1)
Extend the nested (MATNEST) matrix support to all the arithmetics: the
psb_{s,c,z}_nest_{mat,base_mat,tools,builder}_mod modules and the
psb_{s,c,z}_nest_mod umbrellas are generated from the template-psblas
X_nest_* templates; the d sources are regenerated byte-identical.
Preparatory changes to the d sources for clean templating: rowsum/arwsum and
colsum/aclsum no longer share a helper (for the complex arithmetics the
absolute sums are real-valued while the plain sums are complex-valued), the
transposed kernel forwards the actual 'T'/'C' character to the blocks
(conjugate transpose for the complex types), and the capacity helper takes a
type-neutral name.
Build hooks (autotools Makefile and CMakeLists) updated with the per-arith
objects, compile rules and dependencies. All four d tests keep passing.
Author: Simone Staccone (Stack-1)
Add get_owned_rows(i_field) and get_owned_row_count(i_field) to
psb_d_nest_matrix: the list of GLOBAL row indices of a field owned by the
calling process (i.e. the rows it is expected to insert through ins) and
their count. They replace the descriptor-level idiom
field_desc(i)%get_local_rows() / field_desc(i)%l2g(...) in user code, which
leaked descriptor jargon into the build loop.
The high-level tests (glob, rect, builder) are rewritten on the new queries;
the low-level CG test intentionally keeps the descriptor path. README updated
with the new queries and an example.
Author: Simone Staccone (Stack-1)
Complete the integration of the nested (MATNEST) operator into the standard
PSBLAS infrastructure:
- Preconditioners: implement get_diag and csgetrow on psb_d_nest_base_mat so
the stock one-level preconditioners build directly on the nested operator
(DIAG through the concatenated block diagonals, BJAC through the
format-agnostic csget path used by the ILU factorizations).
- Configurable block storage: psb_d_nest_rect_block and psb_d_nest_matrix%asb
accept an optional type ('CSR' default, 'CSC', 'COO') or mold (any class
extending psb_d_base_sparse_mat, e.g. the psb_ext ELL/HLL formats); the
operator is format-agnostic since every operation delegates to the blocks.
- Device-capable matvec: override vect_mv to gather/scatter through the
vectors' own gth/sct with encapsulated index vectors (device kernels on
device vectors) and to run each block through its vect_mv, so device block
formats execute their native kernels; bit-equivalent to csmv on host.
- Full psb_d_base_sparse_mat contract by delegation to the blocks: transposed
csmv (dedicated kernel, ghost contributions left to the transposed halo
exchange), multi-RHS csmm, cp_to_coo/mv_to_coo (unlocking cscnv, csclip,
tril/triu through the base generics), rowsum/arwsum/colsum/aclsum,
maxval/spnmi/spnm1, scal (left/right) and scals, clone (view semantics:
shared blocks, re-owned index maps), mold, sizeof. cp_from_coo/mv_from_coo,
csput and cssv/cssm are intentionally left to the base error (meaningless
for a block-operator view), documented in the type and in the README.
Tests: glob assembles the blocks in HLL (psb_ext) and rect in CSC, both still
bit-identical to the monolithic CSR oracle; the CG test solves under NONE,
DIAG and BJAC/ILU(0), requiring convergence to the exact solution for all of
them and DIAG bit-identical to NONE (exactness check of the nested get_diag).
README updated with the user API reference, the preconditioner section and
the implemented-contract section.
Author: Simone Staccone (Stack-1)
The MINRES implementation merged from development calls the vector psb_spmm
and prec%apply with work=aux, but communication_v2 removed the work argument
from the vector interfaces. Remove work=aux from psb_{c,d,s,z}minres so the
vector solvers match the communication_v2 interface; the library now builds.
The MINRES implementation merged from development calls the vector psb_spmm
and prec%apply with work=aux, but communication_v2 removed the work argument
from the vector interfaces. Remove work=aux from psb_{c,d,s,z}minres so the
vector solvers match the communication_v2 interface; the library builds and
the nested tests pass.
Propagate the latest development (via communication_v2) onto the nested branch:
brings the GMRES refactor, the stopping-criterion change and the restored work
parameter on top of the nested (MATNEST) matrix support. Clean merge, no
conflicts.
Realign communication_v2 with the latest development (10 commits, including the
GMRES refactor and the stopping-criterion change), keeping the communication_v2
work intact.
Conflict resolution:
- base/modules/Makefile (veryclean): keep communication_v2's '/bin/rm -f *.h'.
- linsolve/impl/psb_{c,d,s,z}rgmres.f90: keep development's variable
declaration (itmax_, naux), consistent with development's refactored GMRES
body which references itmax_.
Add a block-structured distributed operator that presents itself to Krylov
solvers and preconditioners as a single ordinary distributed matrix (the
PSBLAS analogue of PETSc MATNEST), targeting saddle-point systems
M = [[A, B^T], [B, 0]] with possibly rectangular sub-blocks.
Library (base/modules):
- psb_desc_nest_mod, psb_d_nest_mat_mod: grid of per-field descriptors and
per-block sparse storage.
- psb_d_nest_base_mat_mod: psb_d_nest_base_mat, the operator extending
psb_d_base_sparse_mat (local csmv, free, field-split hooks for a future
block preconditioner).
- psb_cd_nest_tools_mod / psb_d_nest_tools_mod: composed global descriptor
with union halo (psb_cd_nest_compose) and rectangular local block builder
(psb_d_nest_rect_block), plus the per-block assembly wrappers.
- psb_d_nest_builder_mod: psb_d_nest_matrix, the user frontend with the
init/ins/asb/free pattern hiding all descriptor/halo/compose/setup
boilerplate.
- psb_d_nest_mod: umbrella module (use psb_d_nest_mod).
Remove the earlier bespoke per-block prototype (comm/psblas/vect modules and
the pde_nest_psblas test) superseded by the single MATNEST design.
Tests (test/nested): glob (square operator vs monolithic CSR oracle), rect
(genuinely rectangular blocks), cg (low-level path, ill-conditioned SPD
red-black Laplacian solved with standard CG), builder (same solve via the
utility), plus a README describing the design and usage. All pass serially
and in parallel, with results invariant to the process count.
Build hooks updated (autotools Makefiles + CMakeLists); the nested tests are
relocated out of test/pdegen into test/nested.
Author: Simone Staccone (Stack-1)
The configure step only searched MPI_Fortran_INCLUDE_PATH for mpi.mod, which
is empty for OpenMPI (the wrapper carries includes internally and ships mpi.mod
under lib/). mpi.mod was thus never copied into the module dir, and the
mpi.mod.stamp rule failed on a clean build. Also probe the wrapper
(--showme:inc/libdirs) and lib/ + include/ relative to the Fortran MPI
compiler. MPICH (used in CI) is unaffected.
The nested layer was imported from an older base where the vector psb_spsm
still took a work buffer. communication_v2 removed work from the psb_x_vect_type
routines, so psb_dspsv_vect has no work argument and the call failed generic
resolution. Remove work from psb_d_nest_spsm (signature, declaration, call).
The configure step only searched MPI_Fortran_INCLUDE_PATH for mpi.mod, which
is empty for OpenMPI (the wrapper carries includes internally and ships mpi.mod
under lib/). mpi.mod was thus never copied into the module dir, and the
mpi.mod.stamp rule failed on a clean build. Also probe the wrapper
(--showme:inc/libdirs) and lib/ + include/ relative to the Fortran MPI
compiler. MPICH (used in CI) is unaffected.
base/CMakeLists.txt listed only psb_comm_rma_mod.F90 from comm_schemes,
omitting psb_comm_schemes_mod, psb_comm_baseline_mod, psb_comm_neighbor_impl_mod
and psb_comm_factory_mod, so their .mod files were never generated under CMake.
Also psb_comm_mod was listed as .F90 while the file is .f90 (breaks on
case-sensitive filesystems). The autotools Makefile was already correct.
Stack-1
sfilippone
federicamontes
fdurastante
sfilippone
jalmerol