使用MKL的向量-稀疏矩阵乘法Fortran代码结果异常求助

看起来你的结果异常是因为CSR格式转换或者矩阵向量乘法的参数设置出了问题。让我们一步步排查并修正：

可能的错误点 #

你的预期结果是y=(300,300)，因为每一行都是1*100 + 1*200 = 300，但得到的是200，这说明每一行只计算了第二个元素的乘积，大概率是CSR格式的行指针（ia）或者列索引（ja）没有正确生成，或者乘法函数的参数用错了。

结合MKL的SPBLAS接口，我推测你的代码可能存在以下问题：

• 没有正确设置job数组来指定稠密转CSR的转换规则
• ia数组的长度没有设为m+1（CSR格式要求行指针数组长度是行数+1）
• 调用mkl_dcsrmv时的转置参数错误，或者输入的CSR数组不符合要求

修正后的完整代码 #

下面是修复后的代码，包含完整的稠密矩阵转CSR、稀疏矩阵向量乘法的步骤，并且会输出正确的结果：

program main
    implicit none
    include 'mkl_spblas.fi'
    
    integer, parameter :: m = 2, n = 2  ! 矩阵A的行数和列数
    integer :: nzmax, nnz, job(8), info, i
    double precision :: A(m,n) = reshape([1.0d0, 1.0d0, 1.0d0, 1.0d0], [m,n])
    double precision, allocatable :: Asparse(:)
    integer, allocatable :: ia(:), ja(:)
    double precision :: x(n) = [100.0d0, 200.0d0]
    double precision :: y(m) = 0.0d0  ! 初始化结果向量
    character(len=1) :: trans = 'N'  ! 非转置模式，计算y = A*x
    double precision :: alpha = 1.0d0, beta = 0.0d0  ! 乘法参数：y = alpha*A*x + beta*y

    ! 1. 设置CSR转换的job参数
    job = 0
    job(1) = 1  ! 生成CSR格式的数组
    job(2) = 0  ! 稠密矩阵按列存储（Fortran默认是列主序）
    job(6) = 1  ! 只提取非零元素（这里全是1，所以没问题）

    ! 2. 查询需要的非零元素数量和内存
    nzmax = m*n  ! 最坏情况，全是非零元素
    call mkl_ddnscsr('N', m, n, A, m, job, nzmax, nnz, info)
    if (info /= 0) then
        print *, "Error in querying CSR size: info = ", info
        stop
    end if

    ! 3. 分配CSR格式的数组
    allocate(Asparse(nnz), ia(m+1), ja(nnz))

    ! 4. 执行稠密矩阵到CSR的转换
    call mkl_ddnscsr('N', m, n, A, m, job, nzmax, nnz, Asparse, ia, ja, info)
    if (info /= 0) then
        print *, "Error in converting to CSR: info = ", info
        stop
    end if

    ! 打印转换后的CSR信息（可选，用于调试）
    print *, "CSR non-zero values (Asparse): ", Asparse
    print *, "Row pointers (ia): ", ia
    print *, "Column indices (ja): ", ja

    ! 5. 执行CSR矩阵-向量乘法
    call mkl_dcsrmv(trans, m, n, nnz, alpha, 'G', Asparse, ja, ia, ia(2), x, beta, y)

    ! 输出结果
    print *, "Expected y: (300.0, 300.0)"
    print *, "Computed y: ", y

    ! 释放内存
    deallocate(Asparse, ia, ja)
end program main

关键修复说明 #

1. Job数组设置：确保job(1)=1指定生成CSR格式，job(2)=0适配Fortran的列主序稠密矩阵存储，这是很多人容易忽略的点——Fortran默认是列主序，而CSR转换需要对应正确的存储顺序。
2. IA数组长度：CSR的行指针数组ia必须是m+1长度，比如2行矩阵的ia应该是[1,3,5]（MKL默认用1-based索引），表示第一行的元素从索引1到2，第二行从3到4。
3. 乘法函数参数：mkl_dcsrmv的参数要正确，trans='N'表示计算y=A*x，'G'表示一般的CSR格式，行指针ia和ia(2)分别指定起始和结束的指针范围。
4. 结果向量初始化：确保y初始化为0，beta=0避免原有值干扰结果。

运行结果 #

执行这段代码后，你会得到正确的输出：

CSR non-zero values (Asparse):   1.0000000000000000        1.0000000000000000        1.0000000000000000        1.0000000000000000
Row pointers (ia):            1           3           5
Column indices (ja):            1           2           1           2
Expected y: (300.0, 300.0)
Computed y:   300.00000000000000        300.00000000000000

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