Questions on dnscsr and coocsr

Js23 · ‎03-29-2025

Hi,

I'm new to fortran, so apologies if I am getting something elementary
wrong!

I am trying to work with mkl_spblas and I'm having challenges with the
matrix conversion routines coocsr and dnscsr. I'm attaching a minimal
example which explains the args I'm using, why, and the output.

I've written down the scipy.sparse code I'm attempting to
replicate. The overarching question is why I'm getting different
results from mkl conversion routines than from scipy.sparse.

MKL version is 2025.0.1.

Thanks so much for your help, and please let me know if you need any
other information.

  ! The goal is to replicate this Python scipy.sparse behavior using
  ! mkl_spblas:
  !
  !   coo = coo_array(([13, 17, 19], ([0, 1, 3], [0, 1, 4])), shape=(4, 5))
  !   print(coo.data)       # [13 17 19]
  !   print(coo.coords[0])  # [0 1 3]
  !   print(coo.coords[1])  # [0 1 4]
  !   csr = csr_array(coo) 
  !   print(csr.data)       # [13 17 19]
  !   print(csr.indptr)     # [0 1 2 2 3] (equiv 1 2 3 3 4)
  !   print(csr.indices)    # [0 1 4]     (equiv 1 2 5)
  !   print(csr.toarray())
  !   # [[13  0  0  0  0]
  !   #  [ 0 17  0  0  0]
  !   #  [ 0  0  0  0  0]
  !   #  [ 0  0  0  0 19]]


  ! compiled with:
  !   ifx -o spblas_example spblas_example.f90  ${MKLROOT}/lib/libmkl_blas95_lp64.a ${MKLROOT}/lib/libmkl_lapack95_lp64.a -Wl,--start-group ${MKLROOT}/lib/libmkl_intel_lp64.a ${MKLROOT}/lib/libmkl_sequential.a ${MKLROOT}/lib/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl  -I${MKLROOT}/include/mkl/intel64/lp64 -I"${MKLROOT}/include" && ./spblas_example
  
program spblas_example
  use mkl_spblas
  implicit none
  character*198 :: version_string
  integer, parameter :: nrows = 6
  integer, parameter :: ncols = 4
  integer, parameter :: nnz = 3
  integer :: info, stat, i, j, lda
  integer :: coo_row_idx(nnz), coo_col_idx(nnz)
  real :: coo_values(nnz), csr_values(nnz)
  integer :: job(6)
  ! from dnscsr doc about ja: "Its length is equal to the length of the array acsr."
  integer :: csr_ja(nnz)
  ! from dnscsr doc about ia: "Array of length m + 1"; m is nrows
  integer :: csr_ia(nrows + 1)
  real :: dns(nrows, ncols)
  integer :: dim_mat(2)

  call mkl_get_version_string(version_string)
  print *, version_string
  !  Intel(R) oneAPI Math Kernel Library Version 2025.0.1-Product Build 20241031 for
  !   Intel(R) 64 architecture applications                                         
  
  coo_row_idx = [1, 2, 4]
  coo_col_idx = [1, 2, 5]
  coo_values = [13, 17, 19]

  !! Convert from coo to csr
  !   https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-0/mkl-csrcoo.html
  !     if job(1)=1, the matrix in the coordinate format is converted to the CSR format.
  !     if job(2)=1, one-based indexing for the matrix in CSR format is used.
  !     if job(3)=1, one-based indexing for the matrix in coordinate format is used.
  !     job(5)=nzmax - maximum number of the non-zero elements allowed **if job(1)=0**.
  !     job(6) - job indicator.
  !   For conversion to the CSR format:
  !     If job(6)=0, all arrays acsr, ja, ia are filled in for the output storage.
  
  job = [1, 1, 1, 0, 0, 0]
  dim_mat = [ncols, nrows]
  call mkl_scsrcoo(job, dim_mat, csr_values, csr_ja, csr_ia, nnz, coo_values, &
                   coo_row_idx, coo_col_idx, info)
  print *, 'info', info
  print *, 'csr_ja', csr_ja
  print *, 'csr_ia', csr_ia
  print *, 'csr_values', csr_values
  ! info           0
  ! csr_ja           1           2           5           0           0           0
  ! csr_ia           1           2           3           3
  ! csr_values   13.00000       17.00000       19.00000    

  ! NOTE: These are different (leaving aside the one-based indexing)
  !       from the way the Python CSR indices are encoded, see initial
  !       comment. I'm not sure why this would be the case.
  
  !! Convert from csr to dns
  !   https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-0/mkl-dnscsr.html
  !     if job(1)=1, the rectangular matrix A is restored from the CSR format.
  !     if job(2)=1, one-based indexing for the rectangular matrix A is used.
  !     if job(3)=1, one-based indexing for the matrix in CSR format is used.
  !     If job(4)=2, adns is a whole matrix A.
  !   These should be irrelevant:
  !     job(5)=nzmax: maximum number of the non-zero elements allowed **if job(1)=0.**
  !     job(6): job indicator **for conversion to CSR format.**
  
  job = [1, 1, 1, 2, 0, 0]
  ! job = [1, 1, 1, 2, nnz, 0] ! this generates the same results
  
  ! from doc: Specifies the leading dimension of adns as declared in
  ! the calling (sub)program.  For one-based indexing of A, lda must
  ! be at least max(1, m).
  
  lda = 1 ! NOTE: not sure that this is the right value
  
  call mkl_sdnscsr(job, nrows, ncols, dns, lda, csr_values, csr_ja, csr_ia, info)

  ! from doc:
  ! info
  !   INTEGER. Integer info indicator **only for restoring the matrix A from the CSR format.**
  !   If info=0, the execution is successful.
  !   If info=i, the routine is interrupted processing the i-th row
  !   because there is no space in the arrays acsr and ja according to
  !   the value nzmax.
  !
  ! NOTE: Documentation states that "info" is for *restoring* from CSR
  !       (which is what we are doing). info>0 implies the routine is
  !       interrupted because there is no room in acsr and ja
  !       according to nzmax. However, nzmax is set in job(5), and
  !       this only relevant "if job(1)=0", which is coo->csr
  !       conversion (the other way). In practice, setting the 5th
  !       element of job to nnz doesn't change the output. In my view
  !       (possibly mistaken) this would make more sense if we were
  !       converting *to* CSR because then we could plausibly run out
  !       of room.
  
  print *, 'info', info
  ! info           4

  ! the resulting dense matrix
  
  print *, 'dns'
  do i = 1, nrows
    write(*, *) (dns(i, j), j = 1, ncols)
  end do
  ! dns
  !   13.00000      0.0000000E+00  0.0000000E+00  0.0000000E+00
  !  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  !   17.00000      0.0000000E+00  0.0000000E+00  0.0000000E+00
  !  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  !  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  !  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  
end program spblas_example

noffermans · ‎04-07-2025

Hi Js23,

Thanks for reaching out. I think that the code needs a few fixes:

I believe that there is a mismatch in the matrix dimensions. The COO matrix in the Python example has dimensions (4,5), which seems to indicate 4 rows and 5 columns. Yet, in the Fortran example, the dimensions are nrows=6 and and ncols=4.
The n argument for the mkl_?csrcoo should be a single integer, not an array like dim_mat is. The documentation is a bit ambiguous but n should really be the number of rows here (nrows). The logic being that a matrix does not strictly require knowledge of ncols to be stored in CSR format. This is because the array dimensions (csr_ia(nrows+1), csr_ja(nnz) and csr_a(nnz)) do not depend on ncols.
The leading dimension for the A matrix (lda) must be at least max(1, m) = max(1, nrows) = nrows here. The leading dimension of a matrix is a way to let the library know about the memory layout of the arrays that are passed to / returned by it. In the case of 1-indexing, the Fortran convention where matrix elements are stored in memory contiguously along the columns (column-major order) is assumed. Therefore, the leading dimension for column-major layout should be at least the number of rows.

Please let me know if those fixes resolve your issues.
Note that the info output should be equal to 0 upon successful completion of a conversion. A return value of 4 indicates some error, which the documentation does not specify unfortunately. In this case, this is because the routine encounters a column index of 5 when the matrix only has 4 columns (because of the mismatch in nrows/ncols).

Best,
Nicolas