mkl中f95接口在64位下结果出错

ksfengjia

问题：使用mkl的f95接口在32位计算结果正确，64位计算结果就错误了
环境：vs2015+ivf2017,环境配置正确，均可以编译通过
plus：采用mkl的f77接口在32位和64位均计算正确。
代码：main主程序，solve为高斯解法作为对比

[Fortran] 纯文本查看 复制代码

program main
    use f95_precision
    use lapack95

    implicit none
    !---------------------------------------------------------------------------
    integer(4)            :: info
    integer(4), parameter :: m = 4
    integer(4)            :: i,j
    real(8)               :: a(m,m),aa(m,m),b(m)
    integer(4)            :: p(m)
    real(8)               :: cb1,ce1,cb2,ce2
    !---------------------------------------------------------------------------
    data ((a(i,j),j=1,m),i=1,m) &
        / 1.8 ,  2.88,  2.05, -0.89, &
          5.25, -2.95, -0.95, -3.8 , &
          1.58, -2.69, -2.9 , -1.04, &
         -1.11, -0.66, -0.59,  0.8 /
    data b /9.52, 24.35, 0.77, -6.22/
    
    !---------------------------------------------------------------------------
    ! output the A, b
    write(6,*) 'A ='
    do i = 1,m
        write(6,'(1x,<m-1>(f5.2,2x),f5.2)') (a(i,j),j=1,m)
    end do
    write(6,*) 'b ='
    do i = 1,m
        write(6,'(1x,f5.2)') b(i)
    end do
    
    !---------------------------------------------------------------------------
    ! guass solver
    write(6,*) '=============================='
    write(6,*) 'User defined solver:'
    aa = a
    call cpu_time(cb1)
    call solve(m,aa,b)
    call cpu_time(ce1)
    write(6,*) 'x ='
    do i = 1,m
        write(6,'(1x,f5.2)') b(i)
    end do
    !
!    b = matmul(a,b)
!    write(6,*) 'Ax ='
!    do i = 1,m
!        write(6,'(1x,f5.2)') b(i)
!    end do
    write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce1-cb1
    
    
    !---------------------------------------------------------------------------
    ! mkl solver
    write(6,*) '=============================='
    write(6,*) 'MKL solver:'
    call cpu_time(cb2)
! f77 solution    
!    call dgetrf(4, 4, a, 4, p, info )
!    call dgetrs('N', m, 1, a, m, p, b, m, info)
    ! f95 solution
    call getrf(a,p)
    call getrs(a,p,b)
    
    call cpu_time(ce2)
    
    write(6,*) 'x ='
    do i = 1,m
        write(6,'(1x,f5.2)') b(i)
    end do
!    !
!    b = matmul(a,b)
!    write(6,*) 'Ax ='
!    do i = 1,m
!        write(6,'(1x,f5.2)') b(i)
!    end do
    write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce2-cb2
    !---------------------------------------------------------------------------
    
    continue
    
end program

[Fortran] 纯文本查看 复制代码

subroutine solve(n,a,b)
    !
    implicit none
    !
    integer(kind=4) :: i,j,k,n
    real(kind=8)    :: a(n,n),b(n)
    real(kind=8)    :: x
    !
    do i = 1,n
        continue
        do j = i+1,n
            x = a(j,i)/a(i,i)
            a(j,i) = x
            do k = i+1,n
                a(j,k) = a(j,k) - x*a(i,k)
            end do
        end do
    end do
    !
    do i = 2,n
        do j = 1,i-1
            b(i) = b(i) - b(j)*a(i,j)
        end do
    end do
    !
    b(n) = b(n)/a(n,n)
    do i = n-1,1,-1
        do j = n,i+1,-1
            b(i) = b(i) - b(j)*a(i,j)
        end do
        b(i) = b(i)/a(i,i)
    end do
    !
    return
    !
end subroutine

fcode

你忘了重新给 b 赋值了。

b = [9.52, 24.35, 0.77, -6.22]
!---------------------------------------------------------------------------
! mkl solver