Initializing sparse matrices in Matlab's Mex files (Fortran), and setting their properties (ir, jc, pr) later in the computational routine

Question

I'm trying to initialize a sparse matrix in a Mex file written in Fortran. I saw the example code fulltosparse. There, the number of nonzero elements and the matrix dimensions are set/guessed before entering the computational subroutine that modifies the matrix arrays of indices jc, ir, and pr. In my case, I don't know beforehand how many zeros there will be in the matrix. This question follows directly this one, but here the question revolves strictly around how Mex interface works. I tried coding something like this example that deals with handling numeric arrays of unknown length. After compiling the file in Matlab (Intel compiler), when I try to run it as S = mex_function_name(ii,jj,kk), Matlab crashes. This is what I coded. Thank you

  #include "fintrf.h"

  subroutine mexFunction(nlhs, plhs, nrhs, prhs)
  implicit none
  mwPointer plhs(*), prhs(*)
  integer (kind = 4) nlhs, nrhs
  mwPointer mxGetN, mxCreateSparse, mxGetPr, mxGetIr, mxGetJc
  mwPointer  :: ir, jc, pr
  mwSize     :: n, m=4, zero=0, nnew
  integer(4) :: mxREAL = 0
  
  n = mxGetN(prhs(1)) 
  
  plhs(1) = mxCreateSparse(zero,zero,zero,mxREAL) 
  pr = mxGetPr(plhs(1))
  ir = mxGetIr(plhs(1))
  jc = mxGetJc(plhs(1))
  
  call new_sparse(%val(mxGetPr(prhs(1))), %val(mxGetPr(prhs(2))), %val(mxGetPr(prhs(3))), pr, ir, jc, n, m, nnew)
  
  call mxSetM(plhs(1),m)
  call mxSetN(plhs(1),m)
  
  call mxSetNzmax(plhs(1),nnew)

  ! call mxSetPr(plhs(1), mxRealloc(pr, nnew*8));
  ! call mxSetIr(plhs(1), mxRealloc(ir, nnew*1));
  ! call mxSetJc(plhs(1), mxRealloc(jc, (m+1)*1));

  call mxSetPr(plhs(1),pr) 
  call mxSetIr(plhs(1),ir)
  call mxSetJc(plhs(1),jc)
  
  return
  end

!============================================================
subroutine new_sparse(MatI, MatJ, MatK, pr, ir, jc, n, m, nnew)

implicit none
real(kind = 8),  intent(in) :: MatK(n),     MatI(n),     MatJ(n)
mwSize    :: nnew, n,m
mwPointer :: pr, ir, jc
mwPointer :: mxMalloc
real   (kind = 8), allocatable :: pr_all(:)
integer(kind = 4), allocatable :: ir_all(:), jc_all(:)
integer(kind = 4) ::  i, k, col, row, c, r 
integer(kind = 4) ::  hcol(m+1), jcS(m+1), jrS(m+1)
integer(kind = 4) ::  ixijs, irank(n), rankk(n) 

hcol = 0
jcS  = 0
jrS  = 0    

do i = 1,n
    jrS(MatI(i)+1) = jrS(MatI(i)+1)+1
end do
do r = 2,m+1
    jrS(r) = jrS(r) + jrS(r-1)
end do

do  i = 1,n 
     rankk(jrS(MatI(i))+1) = i
     jrS(MatI(i)) = jrS(MatI(i)) + 1
end do

k = 1
do row = 1,m
    do i = k , jrS(row) 
        ixijs = rankk(i)
        col = MatJ(ixijs)
        if (hcol(col) < row) then
            hcol(col) = row
            jcS(col+1) = jcS(col+1)+1
        end if
        irank(ixijs) = jcS(col+1)
        k = k+1
    end do
end do

do c = 2,m+1  
    jcS(c) = jcS(c) + jcS(c-1)
end do

do i = 1,n
    irank(i) = irank(i) + jcS(MatJ(i))
end do    

nnew = jcS(m+1)

allocate(pr_all(nnew))
allocate(ir_all(nnew))
allocate(jc_all(m+1))

jc_all = jcS
ir_all(irank) = MatI - 1
pr_all = 0
do i = 1,n
    pr_all(irank(i)) = pr_all(irank(i)) +  MatK(i)
end do

pr = mxMalloc( 8 * nnew )
ir = mxMalloc( 8 * nnew )
jc = mxMalloc( 8 * (m+1) )
call mxCopyReal8ToPtr   (pr_all, pr, nnew)
call mxCopyInteger4ToPtr(ir_all, ir, nnew)
call mxCopyInteger4ToPtr(jc_all, jc, (m+1))

deallocate(pr_all)
deallocate(ir_all)
deallocate(jc_all)

return    
end

Edit #1:
I practiced at the beginning writing a function that receives a triplet of indices (ii and jj) and values (kk) of length n and produces another triplet of length nnew as an output, with their values rearranged and with the repetitions eliminated. The function [iinew, jjnew, kknew]=mex_function_name(ii,jj,kk) works fine. For example, if my original triplet is

ii = [3, 4, 1, 3, 2, 1, 4, 4, 4, 3, 2, 3, 1]
jj = [3, 3, 1, 4, 1, 1, 4, 3, 1, 3, 2, 2, 4]
kk = [4, 4, 5, 7, 3, 5, 5, 4, 3, 4, 9, 7, -2]

The output is

iinew = [1, 2, 4, 2, 3, 3, 4, 1, 3, 4]
jjnew = [1, 1, 1, 2, 2, 3, 3, 4, 4, 4]
kknew = [10, 3, 3, 9, 7, 8, 8, -2, 7, 5]

And this is the code:

  #include "fintrf.h"

  subroutine mexFunction(nlhs, plhs, nrhs, prhs)

  implicit none

  mwPointer plhs(*), prhs(*)
  integer (kind = 8) nlhs, nrhs

  mwPointer mxGetN, mxCreateDoubleMatrix, mxGetPr     

  mwPointer MatI_out, MatJ_out, MatK_out
  mwSize     :: n, m, nnew, zero=0
  integer(4) :: mxREAL = 0
  
  n = mxGetN(prhs(1)) 
  m = 4   

  plhs(1) = mxCreateDoubleMatrix(zero,zero,mxREAL)     
  plhs(2) = mxCreateDoubleMatrix(zero,zero,mxREAL)         
  plhs(3) = mxCreateDoubleMatrix(zero,zero,mxREAL) 

  MatI_out = mxGetPr(plhs(1))
  MatJ_out = mxGetPr(plhs(2))  
  MatK_out = mxGetPr(plhs(3))

  call new_sparse(%val(mxGetPr(prhs(1))), %val(mxGetPr(prhs(2))), %val(mxGetPr(prhs(3))), MatI_out, MatJ_out, MatK_out, n, m, nnew)
  
  call mxSetPr(plhs(1),MatI_out) 
  call mxSetM(plhs(1),1)
  call mxSetN(plhs(1),nnew)
  call mxSetPr(plhs(2),MatJ_out) 
  call mxSetM(plhs(2),1)
  call mxSetN(plhs(2),nnew)
  call mxSetPr(plhs(3),MatK_out) 
  call mxSetM(plhs(3),1)
  call mxSetN(plhs(3),nnew)

  return
  end

 !======================================================================
 !    Computational subroutine

subroutine new_sparse(MatI, MatJ, MatK, MatI_out, MatJ_out, MatK_out, n, m, nnew)

implicit none
real(kind = 8),  intent(in) :: MatK(n),     MatI(n),     MatJ(n)
mwPointer :: MatK_out, MatI_out, MatJ_out
mwPointer :: mxMalloc
real(kind = 8), allocatable :: MatK_all(:), MatI_all(:), MatJ_all(:)
mwSize  n, m, nnew
integer ::  i, k, col, row, c, r 
integer ::  hcol(m+1), jcS(m+1), jrS(m+1)
integer ::  ixijs, irank(n), rankk(n)   

... ! same computational subroutine as before

MatI_out = mxMalloc( 8 * nnew )
MatJ_out = mxMalloc( 8 * nnew )
MatK_out = mxMalloc( 8 * nnew )
call mxCopyReal8ToPtr(MatI_all, MatI_out, nnew)
call mxCopyReal8ToPtr(MatJ_all, MatJ_out, nnew)
call mxCopyReal8ToPtr(MatK_all, MatK_out, nnew)
deallocate(MatI_all)
deallocate(MatJ_all)
deallocate(MatK_all)

return
  
end