fortran使用mkl函数库中的mkl_dcsrsymv计算矩阵与向量的乘积

用fortran语言编写数值程序时，如果要计算一个大型稀疏矩阵与一个向量的乘积，可以使用下面高效的方法
1. 首先使用CSR格式大型稀疏矩阵进行存储
2. 调用mkl函数库中的mkl_dcsrsymv计算
这样做的目的不仅可以节约内存，而且计算速度也比较快

示例代码如下：
program mkl_symv
    implicit none
    integer, parameter        :: n = 5  !// depend on your equation
    integer                   :: i, j, mm, tmp, fileid, first, num
    real(kind=8)              :: matrix(n,n), x(n), y(n)
    real(kind=8), allocatable :: csra(:)
    integer, allocatable      :: ia(:), ja(:)
    character(len=1)          :: uplo = 'u'
  
    type                      :: node
        real(kind=8)          :: s
        integer               :: k1, k2
        type (node), pointer  :: next
    end type node 
    type (node), pointer      :: head, tail, p, q  

    open ( newunit = fileid, file = 'SparseMatrix.txt' )  !// The sparse matrix data needs to be given by itself
    do i = 1, n
            read ( fileid,* ) matrix(i,:)
    end do
    read( fileid,* ) x
    close ( fileid )

    !// =========================store data by CSR format=======================
    first = 1
    if ( .not.associated(p) ) then
        allocate( p )
        q    => p
        nullify( p%next )
        q%k2 = first
        tmp  = q%k2
    end if
  
    num = 0  !// calculate the number of no-zero
    do i =  1, n
        mm = 0  !// calculate the position of no-zero in every row
        do j = i, n
            If ( matrix(i,j) /= 0.0 ) then
                if ( .not.associated(head) ) then
                    allocate( head )
                    tail    => head
                    nullify( tail%next )
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1
                    mm      = mm + 1
                else
                    allocate( tail%next )
                    tail    => tail%next
                    tail%s  = matrix(i,j)
                    tail%k1 = j
                    num     = num + 1  
                    mm      = mm + 1  
                    nullify( tail%next )
                end if
            End if
        end do
        if ( associated(p) ) then
            allocate( q%next )
            q    => q%next
            q%k2 = tmp + mm
            tmp  = q%k2
            nullify( q%next )
        end if
    end do
  
    allocate( csra(num), ja(num), ia(n+1) )  !// store the no-zero element
    !// output a and ja
    tail => head
    j = 1
    do 
        if ( .not.associated(tail) ) exit
        csra(j) = tail%s
        ja(j) = tail%k1
        j     = j + 1
        tail  => tail%next
    end do
  
    !// output ia
    q => p
    j = 1
    do 
        if ( .not.associated(q) ) exit
        ia(j) = q%k2
        j     = j + 1
        q     => q%next
    end do
  
    nullify( head, tail, p, q )

    write ( *,'(a)' ) '>> Original data'
    write ( *,'(5f7.2)' ) matrix
    write ( *,'(a)' ) '>> CSR data'
    write ( *,'(*(f7.2))' ) csra
    write ( *,'(a)' ) '>> X data'
    write ( *,'(*(f7.2))' ) x
    write ( *,'(a)' ) '>> Colnum of CSR data'
    write ( *,'(*(i4))' ) ja
    write ( *,'(a)' ) '>> The position of CSR data'
    write ( *,'(*(i4))' ) ia
    
    write ( *,'(1x,a)' ) "The result of matmul function is..."
    write ( *,'(*(f9.3))' ) matmul(matrix,x)
    call mkl_dcsrsymv( uplo, n, csra, ia, ja, x, y )
    print*
    write ( *,'(1x,a)' ) "The result of mkl_dcsrsymv is..."
    write ( *,'(*(f9.3))' ) y
  
end program mkl_symv

结果如下：
>> Original data
   1.00  -1.00   0.00  -3.00   0.00
  -1.00   5.00   0.00   0.00   0.00
   0.00   0.00   4.00   6.00   4.00
  -3.00   0.00   6.00   7.00   0.00
   0.00   0.00   4.00   0.00  -5.00
>> CSR data
   1.00  -1.00  -3.00   5.00   4.00   6.00   4.00   7.00  -5.00
>> X data
 -13.00   9.00  56.00  43.00 -13.00
>> Colnum of CSR data
   1   2   4   2   3   4   5   4   5
>> The position of CSR data
   1   4   5   8   9  10
 The result of matmul function is...
 -151.000   58.000  430.000  676.000  289.000

 The result of mkl_dcsrsymv is...
 -151.000   58.000  430.000  676.000  289.000
从结果可以看出，计算正确

代码中用到的示例数据存储在SparseMatrix.txt中，如下所示：
   1.00  -1.00   0.00  -3.00   0.00
  -1.00   5.00   0.00   0.00   0.00
   0.00   0.00   4.00   6.00   4.00
  -3.00   0.00   6.00   7.00   0.00
   0.00   0.00   4.00   0.00  -5.00
  -13.00   9.00  56.00  43.00 -13.00
前五行为矩阵m，最后一行为向量v