matlab - 在 Matlab 的 Mex 文件 (Fortran) 中初始化稀疏矩阵，并稍后在计算例程中设置它们的属性 (ir, jc, pr)

Question

我正在尝试在用 Fortran 编写的 Mex 文件中初始化一个稀疏矩阵。我看到了示例代码fulltosparse。jc在那里，在进入修改索引、ir和的矩阵数组的计算子例程之前设置/猜测非零元素的数量和矩阵维数pr。就我而言，我事先不知道矩阵中有多少个零。这个问题直接跟在这个问题之后，但这里的问题完全围绕着 Mex 界面的工作原理。我尝试编写类似this example的代码来处理未知长度的数字数组。在 Matlab（英特尔编译器）中编译文件后，当我尝试将其作为 运行时S = mex_function_name(ii,jj,kk)，Matlab 崩溃。这是我编码的。谢谢

  #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

编辑＃1：
我在开始时练习编写一个函数，该函数接收三个长度的索引（ii和jj）和值（kk），n并产生另一个长度的三元组nnew作为输出，它们的值重新排列并消除了重复。该功能[iinew, jjnew, kknew]=mex_function_name(ii,jj,kk)工作正常。例如，如果我原来的三胞胎是

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]

输出是

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]

这是代码：

  #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
    

标签： matlabfortransparse-matrixmex