如何在Fortran Mex文件中实现A=Sparse(I，J，K)(由三元组生成的稀疏矩阵)？ [英] How to implement A = sparse(I, J, K) (sparse matrix from triplet) in a Fortran mex file?
本文介绍了如何在Fortran Mex文件中实现A=Sparse(I，J，K)(由三元组生成的稀疏矩阵)？的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我正在尝试通过一个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和@veryreflie在下面的评论中建议，最好是直接排序femi、femj和femk。我想，先femi排序(并根据femi排序femj和femk)，然后再排序femj(和femj和femk)才能得到想要的结果。剩下的唯一一件事就是处理重复项。
编辑#2：
我用Engblom and 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
推荐答案
这应该可以工作：
module test
implicit none
! This should probably be whatever floating point format Matlab uses.
integer, parameter :: dp = selected_real_kind(15,300)
contains
subroutine new_sparse(femi, femj, femk, pr, ir, jc, n, m)
integer, intent(in) :: n ! The size of femi, femj, femk.
integer, intent(in) :: m ! The no. of rows (and cols) in the matrix.
integer, intent(in) :: femi(n) ! The input i indices.
integer, intent(in) :: femj(n) ! The input j indices.
real(dp), intent(in) :: femk(n) ! The input values.
real(dp), intent(out) :: pr(n) ! The output values.
integer, intent(out) :: ir(n) ! The output i indices.
integer, intent(out) :: jc(m+1) ! Column j has jc(j+1)-jc(j) non-zero entries
! loop indices.
integer :: a,b
! Initialise jc.
! All elements of `jc` are `1` as the output initially contains no elements.
jc = 1
! Loop over the input elements.
do_a : do a=1,n
associate(i=>femi(a), j=>femj(a), k=>femk(a))
! Loop over the stored entries in column j of the output,
! looking for element (i,j).
do b=jc(j),jc(j+1)-1
! Element (i,j) is already in the output, update the output and cycle.
if (ir(b)==i) then
pr(b) = pr(b) + femk(a)
cycle do_a
endif
enddo
! Element (i,j) is not already in the output.
! First make room for the new element in ir and pr,
! then add the element to ir and pr,
! then update jc.
ir(jc(j+1)+1:jc(m+1)) = ir(jc(j+1):jc(m+1)-1)
pr(jc(j+1)+1:jc(m+1)) = pr(jc(j+1):jc(m+1)-1)
ir(jc(j+1)) = i
pr(jc(j+1)) = k
jc(j+1:) = jc(j+1:) + 1
end associate
enddo do_a
end subroutine
end module
program prog
use test
implicit none
integer, parameter :: n = 14
integer, parameter :: m = 6
integer :: femi(n), femj(n)
real(dp) :: femk(n)
real(dp) :: pr(n)
integer :: ir(n),jc(m+1)
integer :: a,b
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 = real([2, 1, 5, 4, 2, 4, 5, 7, 2, 1, 6, 2, 1, 4], dp)
write(*,*) 'Input:'
do a=1,n
write(*,'(a,i0,a,i0,a,f2.0)') '(',femi(a),',',femj(a),') : ',femk(a)
enddo
write(*,*)
call new_sparse(femi,femj,femk,pr,ir,jc,n,m)
write(*,*) 'Output:'
do a=1,m
do b=jc(a),jc(a+1)-1
write(*,'(a,i0,a,i0,a,f2.0)') '(',ir(b),',',a,') : ',pr(b)
enddo
enddo
end program
这写道：
Input:
(1,2) : 2.
(2,2) : 1.
(3,1) : 5.
(2,1) : 4.
(2,1) : 2.
(4,3) : 4.
(5,3) : 5.
(5,6) : 7.
(4,3) : 2.
(6,1) : 1.
(6,1) : 6.
(5,2) : 2.
(5,2) : 1.
(2,4) : 4.
Output:
(3,1) : 5.
(2,1) : 6.
(6,1) : 7.
(1,2) : 2.
(2,2) : 1.
(5,2) : 3.
(4,3) : 6.
(5,3) : 5.
(2,4) : 4.
(5,6) : 7.
算法的瓶颈来自指令indices = [femi==i .and. femj==j]、any(indices .eqv. .true.)和sum(femk, indices)。这些操作都需要O(n)操作，并且由于这些操作都在一个双循环内，所以子例程的总成本是O(m^2*n)。
我的算法分两个阶段工作。第一阶段do b=jc(j),jc(j+1)-1循环将输入中的每个元素与输出的匹配列中的每个元素进行比较，以获得O(mn)操作的最大成本。如果在输出中找到了输入元素，则会更新值，不再需要执行其他操作。
如果在输出中找不到输入元素，则需要将其添加到输出。这由第二阶段处理，即do b...循环之后的代码。由于这需要移动输出元素以便为新元素腾出空间，因此该阶段最多有O(n'^2)次操作，其中n'是输入中唯一元素的数量，对于稀疏矩阵应该满足n'<=n和n'<<m^2。
m和n应该运行得更快，但它肯定有很大的改进空间。我怀疑使用中间数据结构来存储ir和pr是值得的，这样就可以插入新元素，而不必重新排列所有元素。
这篇关于如何在Fortran Mex文件中实现A=Sparse(I，J，K)(由三元组生成的稀疏矩阵)？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
