在 Matlab 的 Mex 文件Fortran中初始化稀疏矩阵，然后在计算例程中设置它们的属性ir、jc、pr

问题描述

我正在尝试在用 Fortran 编写的 Mex 文件中初始化一个稀疏矩阵。我看到了示例代码 fulltosparse。在那里，在进入修改索引 jc、ir 和 pr 的矩阵数组的计算子程序之前，设置/猜测非零元素的数量和矩阵维度。就我而言，我事先不知道矩阵中有多少个零。这个问题直接跟在 this one 后面，但这里的问题严格围绕 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,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,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,nnew)

implicit none
real(kind = 8),intent(in) :: MatK(n),MatI(n),MatJ(n)
mwSize    :: nnew,m
mwPointer :: pr,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,nnew)
call mxCopyInteger4ToPtr(ir_all,nnew)
call mxCopyInteger4ToPtr(jc_all,(m+1))

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

return    
end

编辑#1：
我在开始时练习编写一个函数，该函数接收一个长度为 ii 的索引（jj 和 kk）和值（n）的三元组，并生成另一个长度为 { 的三元组{1}} 作为输出，它们的值重新排列并消除了重复。函数 nnew 工作正常。 例如，如果我原来的三元组是

[iinew,jjnew,kknew]=mex_function_name(ii,kk)

输出是

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

这是代码：

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

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）