從COO格式轉換爲CSR稀疏矩陣格式時總計重複值

Question

從COO格式轉換爲CSR時，如何高效地總結重複值。是否存在類似於scipy實現（http://docs.scipy.org/doc/scipy-0.9.0/reference/sparse.html）的Fortran子程序？我正在使用英特爾的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 

這應該做你想做的。