如何在Fortran中高效对COO矩阵的重复条目求和？

刚好我之前在处理有限元稀疏矩阵组装时碰到过一模一样的问题，给你几个实用且高效的解决方案：

1. 排序+线性归并（最经典、易实现） #

这是处理COO重复项的标准方法，时间复杂度为O(n log n)，对于有限元场景下的矩阵规模完全够用。核心思路是：

• 把COO的行、列、值打包成一个派生类型数组
• 按「行索引升序，行相同则按列索引升序」排序
• 遍历排序后的数组，将相同(row, col)的数值累加合并

示例代码片段：

module COO_Utils
    implicit none

    type :: COO_Entry
        integer :: row, col
        real(8) :: val  ! 用双精度适配有限元计算
    end type COO_Entry

contains

    ! 自定义排序比较函数（Fortran 2008+支持内置sort调用）
    logical function compare_entries(a, b)
        type(COO_Entry), intent(in) :: a, b
        if (a%row /= b%row) then
            compare_entries = (a%row < b%row)
        else
            compare_entries = (a%col < b%col)
        end if
    end function compare_entries

    function merge_coo_entries(entries) result(merged)
        type(COO_Entry), intent(in) :: entries(:)
        type(COO_Entry), allocatable :: merged(:)
        integer :: i, idx, n
        integer :: current_row, current_col
        real(8) :: current_val

        n = size(entries)
        if (n == 0) return

        ! 先排序
        allocate(merged(n))  ! 先分配最大可能空间，最后裁剪
        merged = entries
        call sort(merged, compare_entries)

        ! 合并重复项
        idx = 1
        current_row = merged(1)%row
        current_col = merged(1)%col
        current_val = merged(1)%val

        do i = 2, n
            if (merged(i)%row == current_row .and. merged(i)%col == current_col) then
                current_val = current_val + merged(i)%val
            else
                merged(idx) = COO_Entry(current_row, current_col, current_val)
                idx = idx + 1
                current_row = merged(i)%row
                current_col = merged(i)%col
                current_val = merged(i)%val
            end if
        end do
        ! 添加最后一个条目
        merged(idx) = COO_Entry(current_row, current_col, current_val)
        ! 裁剪到实际大小
        merged = merged(1:idx)
    end function merge_coo_entries

end module COO_Utils

如果你的编译器不支持Fortran 2008的内置sort，可以自己实现快速排序，或者用Sparse Kit里的排序子程序（它虽然没有合并功能，但排序工具是完备的）。

2. 哈希表快速查找（适合超大规模矩阵） #

如果你的矩阵规模特别大（比如百万级以上的条目），排序的O(n log n)开销可能会比较明显，这时可以用哈希表实现平均O(n)复杂度的合并。核心是用(row, col)作为键，遍历过程中直接查找是否已存在，存在则累加，不存在则插入。

示例简化版实现（用链表解决哈希碰撞）：

function hash_coo_merge(entries, table_size) result(merged)
    type(COO_Entry), intent(in) :: entries(:)
    integer, intent(in) :: table_size  ! 建议选大于最大列数的质数，比如100003
    type(COO_Entry), allocatable :: merged(:)
    type :: Hash_Node
        integer :: row, col
        real(8) :: val
        integer :: next
    end type Hash_Node
    type(Hash_Node), allocatable :: hash_table(:)
    integer, allocatable :: head(:)
    integer :: i, hash_key, current, idx, n
    logical :: found

    n = size(entries)
    allocate(hash_table(table_size), head(table_size))
    hash_table%next = 0
    head = 0

    do i = 1, n
        hash_key = mod(entries(i)%row * 997 + entries(i)%col, table_size) + 1  ! 保证索引从1开始
        current = head(hash_key)
        found = .false.
        ! 遍历链表查找重复项
        do while (current /= 0 .and. .not. found)
            if (hash_table(current)%row == entries(i)%row .and. hash_table(current)%col == entries(i)%col) then
                hash_table(current)%val = hash_table(current)%val + entries(i)%val
                found = .true.
            end if
            current = hash_table(current)%next
        end do
        ! 不存在则插入哈希表
        if (.not. found) then
            ! 找空节点
            do current = 1, table_size
                if (hash_table(current)%next == 0 .and. hash_table(current)%row == 0) exit
            end do
            hash_table(current)%row = entries(i)%row
            hash_table(current)%col = entries(i)%col
            hash_table(current)%val = entries(i)%val
            hash_table(current)%next = head(hash_key)
            head(hash_key) = current
        end if
    end do

    ! 导出哈希表中的条目到COO数组
    idx = 0
    do i = 1, table_size
        if (hash_table(i)%row /= 0) idx = idx + 1
    end do
    allocate(merged(idx))
    idx = 1
    do i = 1, table_size
        if (hash_table(i)%row /= 0) then
            merged(idx) = COO_Entry(hash_table(i)%row, hash_table(i)%col, hash_table(i)%val)
            idx = idx + 1
        end if
    end do
end function hash_coo_merge

注意：哈希表的大小要选合适，一般选条目数的1.5-2倍的质数，避免因碰撞过多降低效率。

3. 借助成熟稀疏矩阵库（减少重复造轮子） #

如果你的项目已经在使用第三方数值库，直接调用库函数会更省心：

• PETSc：通过MatCreateCOO创建矩阵时，设置MAT_NEW_NONZERO_ALLOCATION_ERR为false，或者调用MatSetOption(mat, MAT_KEEP_NONZERO_PATTERN, false)，PETSc会自动合并重复的(row, col)条目并求和。
• Trilinos/Epetra：使用Epetra_CrsMatrix的SumIntoGlobalValues直接累加重复项；如果已有COO数据，可先构造Epetra_COO再转换为Epetra_CrsMatrix，转换过程中会自动合并。

额外优化建议 #

在有限元组装阶段就尽量避免产生重复条目：比如每个单元刚度矩阵的(i,j)项只写入一次，而不是多次追加相同的(row, col)。这样可以从源头减少合并的工作量，提升整体效率。

内容的提问来源于stack exchange，提问作者bharat