从 COO 转换为 CSR 稀疏矩阵格式时对重复值求和

Question

从 COO 格式转换为 CSR 时，如何有效地总结重复值。 fortran 子例程中是否存在类似于 scipy 实现（http://docs.scipy.org/doc/scipy-0.9.0/reference/sparse.html）的东西？我正在使用英特尔的 MKL 辅助例程从 COO 转换为 CSR，但它似乎不适用于重复值。

Answer 1

在我的代码中，我使用我编写的这个子例程：

subroutine csr_sum_duplicates(Ap, Aj, Ax)
! Sum together duplicate column entries in each row of CSR matrix A
! The column indicies within each row must be in sorted order.
! Explicit zeros are retained.
! Ap, Aj, and Ax will be modified *inplace*
integer, intent(inout) :: Ap(:), Aj(:)
real(dp), intent(inout) :: Ax(:)
integer :: nnz, r1, r2, i, j, jj
real(dp) :: x
nnz = 1
r2 = 1
do i = 1, size(Ap) - 1
    r1 = r2
    r2 = Ap(i+1)
    jj = r1
    do while (jj < r2)
        j = Aj(jj)
        x = Ax(jj)
        jj = jj + 1
        do while (jj < r2)
            if (Aj(jj) == j) then
                x = x + Ax(jj)
                jj = jj + 1
            else
                exit
            end if
        end do
        Aj(nnz) = j
        Ax(nnz) = x
        nnz = nnz + 1
    end do
    Ap(i+1) = nnz
end do
end subroutine

并且您可以使用这个子例程对索引进行排序：

subroutine csr_sort_indices(Ap, Aj, Ax)
! Sort CSR column indices inplace
integer, intent(inout) :: Ap(:), Aj(:)
real(dp), intent(inout) :: Ax(:)
integer :: i, r1, r2, l, idx(size(Aj))
do i = 1, size(Ap)-1
    r1 = Ap(i)
    r2 = Ap(i+1)-1
    l = r2-r1+1
    idx(:l) = argsort(Aj(r1:r2))
    Aj(r1:r2) = Aj(r1+idx(:l)-1)
    Ax(r1:r2) = Ax(r1+idx(:l)-1)
end do
end subroutine

其中 argsort 是

function iargsort(a) result(b)
! Returns the indices that would sort an array.
!
! Arguments
! ---------
!
integer, intent(in):: a(:)    ! array of numbers
integer :: b(size(a))         ! indices into the array 'a' that sort it
!
! Example
! -------
!
! iargsort([10, 9, 8, 7, 6])   ! Returns [5, 4, 3, 2, 1]

integer :: N                           ! number of numbers/vectors
integer :: i,imin,relimin(1)           ! indices: i, i of smallest, relative imin
integer :: temp                        ! temporary
integer :: a2(size(a))
a2 = a
N=size(a)
do i = 1, N
    b(i) = i
end do
do i = 1, N-1
    ! find ith smallest in 'a'
    relimin = minloc(a2(i:))
    imin = relimin(1) + i - 1
    ! swap to position i in 'a' and 'b', if not already there
    if (imin /= i) then
        temp = a2(i); a2(i) = a2(imin); a2(imin) = temp
        temp = b(i); b(i) = b(imin); b(imin) = temp
    end if
end do
end function

这应该做你想要的。