稀疏矩阵存储方式¶
约 405 个字 151 行代码 预计阅读时间 3 分钟
这里介绍常用的三种格式：CSR、CSC、COO。
- CSR (Compressed Sparse Row, 压缩行存储)
- 核心思想：只存储非零元素，并记录每行起始位置
- 数组结构：
val(nnz)：所有非零元值，按 行优先 顺序存储。col_ind(nnz)：每个非零元的 列下标。row_ptr(n+1)：第i行的第一个非零元在val中的位置（1-based），row_ptr(n+1)为总非零元数+1。
- CSC (Compressed Sparse Column, 压缩列存储)
- 核心思想：类似CSR，但数据是 按列优先 顺序存储，并记录每列起始位置。
- 数组结构：
val(nnz)：所有非零元值，按 列优先 顺序存储。row_ind(nnz)：每个非零元的 行下标。col_ptr(n+1)：第j列的第一个非零元在val中的位置（1-based），col_ptr(n+1)为总非零元数+1。
- COO (Coordinate, 坐标格式)
- 核心思想：直接记录每个非零元的位置（行、列）和值。
- 数组结构：
row(nnz)：每个非零元的行号。col(nnz)：每个非零元的列号。val(nnz)：每个非零元的值。
以下面这个矩阵为例：
\[
a = \begin{bmatrix} 1 & 0 & 3 \\ 0 & 5 & 0 \\ 7 & 0 & 9 \end{bmatrix}
\]
可分别编制三个子程序用于稠密矩阵的转换
subroutine dense_to_csr(a, n, m, val, col_ind, row_ptr, nnz)
! 稠密矩阵转换为CSR格式
implicit none
integer, intent(in) :: n, m
real(8), intent(in) :: a(n, m)
real(8), allocatable, intent(out) :: val(:)
integer, allocatable, intent(out) :: col_ind(:), row_ptr(:)
integer, intent(out) :: nnz ! 非零元的个数
integer :: i, j, k
nnz = count(abs(a) > 0.0d0)
allocate(val(nnz), col_ind(nnz), row_ptr(n+1))
k = 1
row_ptr(1) = 1
do i = 1, n
do j = 1, m
if (a(i, j) /= 0.0d0) then
val(k) = a(i, j)
col_ind(k) = j
k = k + 1
end if
end do
row_ptr(i+1) = k
end do
end subroutine dense_to_csr
subroutine dense_to_csc(a, n, m, val, row_ind, col_ptr, nnz)
! 稠密矩阵转换为CSC格式
implicit none
integer, intent(in) :: n, m
real(8), intent(in) :: a(n, m)
real(8), allocatable, intent(out) :: val(:)
integer, allocatable, intent(out) :: row_ind(:), col_ptr(:)
integer, intent(out) :: nnz
integer :: i, j, k
nnz = count(abs(a) > 0.0d0)
allocate(val(nnz), row_ind(nnz), col_ptr(m+1))
k = 1
col_ptr(1) = 1
do j = 1, m
do i = 1, n
if (a(i, j) /= 0.0d0) then
val(k) = a(i, j)
row_ind(k) = i
k = k + 1
end if
end do
col_ptr(j+1) = k
end do
end subroutine dense_to_csc
subroutine dense_to_coo(a, n, m, val, row, col, nnz)
! 稠密矩阵转换为COO格式
implicit none
integer, intent(in) :: n, m
real(8), intent(in) :: a(n, m)
real(8), allocatable, intent(out) :: val(:)
integer, allocatable, intent(out) :: row(:), col(:)
integer, intent(out) :: nnz
integer :: i, j, k
! 统计非零元个数
nnz = count(abs(a) > 0.0d0)
allocate(val(nnz), row(nnz), col(nnz))
k = 1
do i = 1, n
do j = 1, m
if (a(i, j) /= 0.0d0) then
val(k) = a(i, j) ! 非零元值
row(k) = i ! 行号
col(k) = j ! 列号
k = k + 1
end if
end do
end do
end subroutine dense_to_coo
主程序调用：
program main
implicit none
integer, parameter :: n = 3, m = 3
real(8), allocatable :: a(:,:)
real(8), allocatable :: val_csr(:), val_csc(:), val_coo(:)
integer, allocatable :: col_ind(:), row_ptr(:)
integer, allocatable :: row_ind(:), col_ptr(:)
integer, allocatable :: row(:), col(:)
integer :: nnz_csr, nnz_csc, nnz_coo, i
! 初始化全矩阵
allocate(a(n,m))
a = reshape([1d0,0d0,3d0,0d0,5d0,0d0,7d0,0d0,9d0],[n,m],order=[2,1])
! a =
! 1 0 3
! 0 5 0
! 7 0 9
! 转为三种稀疏格式
call dense_to_csr(a, n, m, val_csr, col_ind, row_ptr, nnz_csr)
call dense_to_csc(a, n, m, val_csc, row_ind, col_ptr, nnz_csc)
call dense_to_coo(a, n, m, val_coo, row, col, nnz_coo)
! 打印 CSR
print *, "==== CSR ===="
write(*, '(A, *(F6.1))') "val : ", val_csr
write(*, '(A, *(I3))') "col_ind : ", col_ind
write(*, '(A, *(I3))') "row_ptr : ", row_ptr
! 打印 CSC
print *, "==== CSC ===="
write(*, '(A, *(F6.1))') "val : ", val_csc
write(*, '(A, *(I3))') "row_ind : ", row_ind
write(*, '(A, *(I3))') "col_ptr : ", col_ptr
! 打印 COO
print *, "==== COO ===="
write(*, '(A, *(F6.1))') "val : ", val_coo
write(*, '(A, *(I3))') "row : ", row
write(*, '(A, *(I3))') "col : ", col
print *, "==== COO list ===="
write(*,'(A)') " row col value"
do i = 1, nnz_coo
write(*,'(I4, I4, F10.3)') row(i), col(i), val_coo(i)
end do
end program main
输出：
==== CSR ====
val : 1.0 3.0 5.0 7.0 9.0
col_ind : 1 3 2 1 3
row_ptr : 1 3 4 6
==== CSC ====
val : 1.0 7.0 5.0 3.0 9.0
row_ind : 1 3 2 1 3
col_ptr : 1 3 4 6
==== COO ====
val : 1.0 3.0 5.0 7.0 9.0
row : 1 1 2 3 3
col : 1 3 2 1 3
==== COO list ====
row col value
1 1 1.000
1 3 3.000
2 2 5.000
3 1 7.000
3 3 9.000
注意：以上三个子程序仅用于稀疏矩阵存储方式的演示，在实际运用中，如有限元的整体刚度方程求解中，可以预先统计每行/每列非零元数量分配空间，边装配边填COO/CSR/三元组，组装完一次性排序/去重/压缩即可。