如何在 Fortran mex 文件中实现 A = sparse(I, J, K)来自三元组的稀疏矩阵？

问题描述

我正在尝试通过 mex 函数（用 Fortran 编写）在 Matlab 中创建一个稀疏方阵。我想要类似 A = sparse(I,J,K) 的东西。我的三胞胎是这样的，条目之间有重复

femi = [1,2,3,4,5,6,2]
femj = [2,1,4]
femk = [2,7,4]

我写了一段粗略的代码，它适用于小矩阵维度，但比固有的 Matlab 的 sparse 慢得多。因为我几乎没有编码背景，所以我不知道我做错了什么（分配变量的方式错误？循环太多了？）。任何帮助表示赞赏。谢谢你。这是 mex 计算子程序。它返回 pr,ir,jc 索引数组以提供给稀疏矩阵

subroutine new_sparse(femi,femj,femk,pr,jc,n,m)

      implicit none
      intrinsic:: SUM,COUNT,ANY
      integer :: i,j,k,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 进行排序/排序（并排序 {{1 }} 和 femi 根据 femk) 提供所需的结果。剩下的唯一事情就是处理重复项。

编辑 #2：
我用 Engblom and Lukarksi 逐行翻译了 C 代码的序列化版本。该文档非常清楚地解释了他们的推理，我认为它对像我这样的初学者很有用。但是，由于我的经验不足，我无法翻译代码的并行版本。也许这会引发另一个问题。

femj

解决方法

这应该有效：

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,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,4,5,6,2]
  femj = [2,1,4]
  femk = real([2,7,4],dp)
  
  write(*,*) 'Input:'
  do a=1,n
    write(*,'(a,i0,a,f2.0)') '(',femi(a),',femj(a),') : ',femk(a)
  enddo
  
  write(*,*)
  
  call new_sparse(femi,m)
  
  write(*,*) 'Output:'
  do a=1,m
    do b=jc(a),jc(a+1)-1
      write(*,ir(b),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) : 6.
(6,1) : 7.
(1,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 是值得的，这样就可以插入新元素而无需重新排列所有元素。