我正在尝试通过 mex 函数(用 Fortran 编写)在 Matlab 中创建一个稀疏方阵。我想要类似的东西A = sparse(I,J,K)
。我的三胞胎长这样,条目之间有重复
femi = [1, 2, 3, 2, 2, 4, 5, 5, 4, 6, 6, 5, 5, 2]
femj = [2, 2, 1, 1, 1, 3, 3, 6, 3, 1, 1, 2, 2, 4]
femk = [2, 1, 5, 4, 2, 4, 5, 7, 2, 1, 6, 2, 1, 4]
我写了一段粗略的代码,它适用于小矩阵维度,但它比固有的 Matlab 的sparse
. 由于我几乎没有编码背景,所以我不知道我做错了什么(分配变量的方法错误?太多的循环?)。任何帮助表示赞赏。谢谢你。这是 mex 计算子程序。它返回 pr、ir、jc 索引数组以提供给稀疏矩阵
subroutine new_sparse(femi, femj, femk, pr, ir, jc, n, m)
implicit none
intrinsic:: SUM, COUNT, ANY
integer :: i, j, k, n, indjc, m
real*8 :: femi(n), femj(n), femk(n)
real*8 :: pr(n)
integer :: ir(n),jc(m+1)
logical :: indices(n)
indices = .false.
k = 1
indjc = 0
jc(1) = 0
do j=1,m
do i =1,m
indices = [femi==i .and. femj==j]
if (ANY(indices .eqv. .true.)) then
ir(k) = i-1
pr(k) = SUM(femk, indices)
k = k+1
indjc = indjc + 1
end if
end do
if (indjc/=0) then
jc(j+1) = jc(j) + indjc
indjc = 0
else
jc(j+1) = jc(j)
end if
end do
return
end
编辑:
正如用户@jack 和@veryreverie 在下面的评论中所建议的那样,最好直接排序femi
,femj
和femk
. 我猜femi
首先排名/排序(以及排序femj
和femk
根据femi
)然后排名/排序femj
(以及排序femi
和femk
根据femj
)提供了所需的结果。剩下的唯一事情就是处理重复项。
编辑#2:
我逐行翻译了Engblom 和 Lukarksi的 C 代码的序列化版本。该文档非常清楚地解释了他们的推理,我认为它对像我这样的初学者很有用。但是,由于我缺乏经验,我无法翻译代码的并行版本。也许这会引发另一个问题。
subroutine new_sparse(ir, jcS, pr, MatI, MatJ, MatK, n, m)
! use omp_lib
implicit none
integer, parameter :: dp = selected_real_kind(15,300)
integer, intent(in) :: n, m
real(dp), intent(in) :: MatK(n), MatI(n), MatJ(n)
! integer*8, intent(out) :: nnew
integer :: i, k, col, row, c, r !, nthreads
integer :: hcol(m+1), jcS(m+1), jrS(m+1)
integer :: ixijs, irank(n), rank(n)
real*8 :: pr(*)
integer :: ir(*)
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
rank(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 = rank(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
ir(irank) = MatI-1
do i = 1,n
pr(irank(i)) = pr(irank(i)) + MatK(i)
end do
return
end