问题描述
在知道矩阵稀疏的情况下,如何在Fortran中进行矩阵向量乘法?
例如,我想计算3x3矩阵之间的乘积
0 1 2
3 0 0
4 0 5
和向量
1
2
3
我正在考虑存储行和列以及值(当它不等于0却没有让我到任何地方时)。
解决方法
首先,您需要确定稀疏或稠密代数会更快:
- 实际计算中涉及的矩阵大小是多少?如果它是某个
dimension(m,n)
,例如m<=10
和n<=10
,那么无论如何,使用密集代数会更快:matmul(A,b)
或矩阵向量乘法的另一个库实现将完成这项工作。 - 矩阵的稀疏度是多少?在您的示例中,矩阵的稀疏度约为〜44%(
sparsity=count(matrix==0)/product(shape(matrix))
)。通常,如果矩阵的稀疏度至少约为50%,那么使用稀疏代数将比稠密代数更有利。
一旦决定使用稀疏代数,就不要定义自己的稀疏矩阵存储!只需使用plenty available out there之一。如果矩阵向量乘法最多是与稀疏矩阵有关,那么我将使用CSR
。一些简单的实现如下所示:
module sparse_CSR
use iso_fortran_env
implicit none
public
type,public :: CSRMatrix
integer :: m = 0
integer :: n = 0
integer,allocatable :: rowPtr(:)
integer,allocatable :: colPtr(:)
real(real64),allocatable :: aij(:)
end type CSRMatrix
contains
! Convert a dense matrix to a sparse one.
pure type(CSRMatrix) function dense_to_sparse(matrix) result(sparse)
real(real64),intent(in) :: matrix(:,:)
integer :: ierr,m,n,nnz,row,nnz_row,j1,j2,colIndex(size(matrix,2))
logical :: row_nonZeroes(size(matrix,2))
! Get matrix size
m = size(matrix,1)
n = size(matrix,2)
nnz = count(matrix/=0.0_real64)
! Create an array of indices
forall(j1=1:n) colIndex(j1) = j1
! Set matrix size
sparse%m = m
sparse%n = n
! Allocate data
allocate(sparse%rowPtr(m+1),sparse%colPtr(nnz),sparse%aij(nnz))
! Fill rows
sparse%rowPtr(1) = 1
fill_by_rows: do row=1,m
row_nonZeroes = matrix(row,:)/=0.d0
nnz_row = count(row_nonZeroes)
j1 = sparse%rowPtr(row)
j2 = j1+nnz_row-1
sparse%colPtr(j1:j2) = pack(colIndex,row_nonZeroes)
sparse%aij (j1:j2) = pack(matrix(row,:),row_nonZeroes)
! Set pointer for beginning of next line
sparse%rowPtr(row+1) = j2+1
end do fill_by_rows
end function dense_to_sparse
! Compute matrix-vector product
pure function sparse_matvec(matrix,vector) result(MxV)
type(CSRMatrix),intent(in) :: matrix
real(real64),intent(in) :: vector(:)
real(real64) :: MxV(matrix%m)
integer :: row,j,col
real(real64) :: aij
! Check size. Cannot stop in a pure procedure.
! So,return NaNs
if (size(vector,1)/=matrix%n) then
MxV(:) = transfer(-2251799813685248_int64,1._real64)
return
endif
! Compute product by rows
do row=1,matrix%m
j1 = matrix%rowPtr(row)
j2 = matrix%rowPtr(row+1)-1;
MxV(row) = 0.0_real64
do j=j1,j2
col = matrix%colPtr(j)
aij = matrix%aij(j)
MxV(row) = MxV(row) + aij*vector(col)
end do
end do
end function sparse_matvec
end module sparse_CSR
program test_sparse
use sparse_CSR
use iso_fortran_env
implicit none
real(real64) :: dense(3,3) = reshape(real([0,3,4,1,2,5],real64),[3,3])
real(real64) :: vector(3) = real([1,3],real64)
type(CSRMatrix) :: sparse
real(real64) :: MxV(3)
sparse = dense_to_sparse(dense)
MxV = sparse_matvec(sparse,vector)
print *,'MxV = ',MxV
end program test_sparse
这将返回:
MxV = 8.0000000000000000 3.0000000000000000 19.000000000000000
,
除非您想学习如何实现它,否则可以使用稀疏矩阵库。例如。 SPARSKIT或suitesparse。