Fortran中的稀疏矩阵向量积

如何解决Fortran中的稀疏矩阵向量积

在知道矩阵稀疏的情况下，如何在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。