Fortran Coder

查看: 58|回复: 1

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

[复制链接]

30

帖子

13

主题

0

精华

熟手

F 币
120 元
贡献
198 点
发表于 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

回复

使用道具 举报

1286

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂伙计

F 币
415 元
贡献
138 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2018-9-7 20:43:50 | 显示全部楼层
你忘了重新给 b 赋值了。

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

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2018-9-25 15:22

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

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