Took time to study the materials in the shared folder on Google Drive to understand the H matrix, which I now know is a linear, forward interpolation operator matrix generator which should return a sparse block diagonal matrix of one or two partitions (for Infection and Death, respectively).
The rows of the partitions in the matrix are the row-major order, linear form of a sparse matrix with dimensions equal to the square matrix of cells of a SpatRaster used as input.
The cell number in the raster, corresponding to a Euclidean coordinate pair point CRS(lat, lon) of the health zone on Earth's surface, is used to populate the Moore neighbourhood of the sparse matrix, using a Chess queen neighbourhood of zero, one, and two orders (self, first order, second order)
Finished the first refactoring of the linear forward interpolation operator matrix generator, but it has a bug in it; the general solution is complete, however.
I read this for an introduction to the concepts after reading the old code.