Fortran Coder

查看: 7125|回复: 1
打印 上一主题 下一主题

[数学库] mkl中f95接口在64位下结果出错

[复制链接]

35

帖子

17

主题

0

精华

熟手

F 币
136 元
贡献
240 点
跳转到指定楼层
楼主
发表于 2018-9-7 17:01:10 | 只看该作者 回帖奖励 |正序浏览 |阅读模式
问题:使用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

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2015

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1574 元
贡献
676 点

美女勋章热心勋章星光勋章新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

沙发
发表于 2018-9-7 20:43:50 | 只看该作者
你忘了重新给 b 赋值了。

b = [9.52, 24.35, 0.77, -6.22]
!---------------------------------------------------------------------------
! mkl solver
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-11-11 02:12

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表