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

ksfengjia · 发表于 2018-9-7 17:01:10

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

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

01program main
02    use f95_precision
03    use lapack95
04 
05    implicit none
06    !---------------------------------------------------------------------------
07    integer(4)            :: info
08    integer(4), parameter :: m = 4
09    integer(4)            :: i,j
10    real(8)               :: a(m,m),aa(m,m),b(m)
11    integer(4)            :: p(m)
12    real(8)               :: cb1,ce1,cb2,ce2
13    !---------------------------------------------------------------------------
14    data ((a(i,j),j=1,m),i=1,m) &
15        / 1.8 ,  2.88,  2.05, -0.89, &
16          5.25, -2.95, -0.95, -3.8 , &
17          1.58, -2.69, -2.9 , -1.04, &
18         -1.11, -0.66, -0.59,  0.8 /
19    data b /9.52, 24.35, 0.77, -6.22/
20     
21    !---------------------------------------------------------------------------
22    ! output the A, b
23    write(6,*) 'A ='
24    do i = 1,m
25        write(6,'(1x,<m-1>(f5.2,2x),f5.2)') (a(i,j),j=1,m)
26    end do
27    write(6,*) 'b ='
28    do i = 1,m
29        write(6,'(1x,f5.2)') b(i)
30    end do
31     
32    !---------------------------------------------------------------------------
33    ! guass solver
34    write(6,*) '=============================='
35    write(6,*) 'User defined solver:'
36    aa = a
37    call cpu_time(cb1)
38    call solve(m,aa,b)
39    call cpu_time(ce1)
40    write(6,*) 'x ='
41    do i = 1,m
42        write(6,'(1x,f5.2)') b(i)
43    end do
44    !
45!    b = matmul(a,b)
46!    write(6,*) 'Ax ='
47!    do i = 1,m
48!        write(6,'(1x,f5.2)') b(i)
49!    end do
50    write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce1-cb1
51     
52     
53    !---------------------------------------------------------------------------
54    ! mkl solver
55    write(6,*) '=============================='
56    write(6,*) 'MKL solver:'
57    call cpu_time(cb2)
58! f77 solution    
59!    call dgetrf(4, 4, a, 4, p, info )
60!    call dgetrs('N', m, 1, a, m, p, b, m, info)
61    ! f95 solution
62    call getrf(a,p)
63    call getrs(a,p,b)
64     
65    call cpu_time(ce2)
66     
67    write(6,*) 'x ='
68    do i = 1,m
69        write(6,'(1x,f5.2)') b(i)
70    end do
71!    !
72!    b = matmul(a,b)
73!    write(6,*) 'Ax ='
74!    do i = 1,m
75!        write(6,'(1x,f5.2)') b(i)
76!    end do
77    write(6,'(1x,a,f20.10)') 'Spent CPU time:', ce2-cb2
78    !---------------------------------------------------------------------------
79     
80    continue
81     
82end program

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

01subroutine solve(n,a,b)
02    !
03    implicit none
04    !
05    integer(kind=4) :: i,j,k,n
06    real(kind=8)    :: a(n,n),b(n)
07    real(kind=8)    :: x
08    !
09    do i = 1,n
10        continue
11        do j = i+1,n
12            x = a(j,i)/a(i,i)
13            a(j,i) = x
14            do k = i+1,n
15                a(j,k) = a(j,k) - x*a(i,k)
16            end do
17        end do
18    end do
19    !
20    do i = 2,n
21        do j = 1,i-1
22            b(i) = b(i) - b(j)*a(i,j)
23        end do
24    end do
25    !
26    b(n) = b(n)/a(n,n)
27    do i = n-1,1,-1
28        do j = n,i+1,-1
29            b(i) = b(i) - b(j)*a(i,j)
30        end do
31        b(i) = b(i)/a(i,i)
32    end do
33    !
34    return
35    !
36end subroutine