Performing storage of banded matrix in Fortran

Question

I wrote a derived data type to store banded matrices in Compressed Diagonal Storage format; in particular I store each diagonal of the banded matrix in a column of the 2D array cds(1:N,-L:U), where N is the number of rows of the full matrix and L and U are the number of lower and upper diagonals (this question includes the definition of the type).

I also wrote a function to perform the product between a matrix in this CDS format and a full vector. To obtain each element of the product vector, the elements of the corresponding row of cds are used, which are not contiguous in memory, since the language is Fortran. Because of this I was wandering if a better solution would be to store the diagonals in the rows of a 2D array cds2(-L:U,1:N), which seems pretty reasonable to me.

On the contrary here I read

we can allocate for the matrix A an array val(1:n,-p:q). The declaration with reversed dimensions (-p:q,n) corresponds to the LINPACK band format [132], which, unlike compressed diagonal storage (CDS), does not allow for an efficiently vectorizable matrix-vector multiplication if p + q is small.

Which is just what seems appropriate to C in my opinion. What am I missing?

EDIT

The core of the routine performing matrix vector products is the following

DO i = A%lb(1), A%ub(1)
    CDS_mat_x_full_vec(i) = DOT_PRODUCT(A%matrix(i,max(-lband,lv-i):min(uband,uv-i)), &
                                      & v(max(i-lband,lv):min(i+uband,uv)))
END DO

(Where lv and uv are used to take into account the case of the vector indexed from an index other than 1.) The matrix A is then accessed by rows.

Answer 1

I implemented the derived type which stores the diagonals in an array val(-p:q,1:n) and it is faster, as I supposed. So I think that the link I referenced refers to a row major storage language as C and not a column major one as Fortran. (Or it implements the matrix product in a way I can't even imagine.)