Fortran Coder

查看: 2054|回复: 3

[数学库] MKLFFT

[复制链接]

24

帖子

7

主题

0

精华

熟手

F 币
143 元
贡献
88 点
发表于 2021-12-27 22:55:10 | 显示全部楼层 |阅读模式
[Fortran] 纯文本查看 复制代码
    program TEST_FFTwithMatlab
    USE MKL_DFTI

    IMPLICIT NONE
    INTEGER*4::X_in_Len = 7
    COMPLEX*8::X_in(7)
    INTEGER*4::Dim=1,X_Len,Status,i
    REAL*8,DIMENSION(7)::a,b
    COMPLEX*8,DIMENSION(:),ALLOCATABLE::X_inout
    Type(DFTI_DESCRIPTOR), POINTER :: Desc_handle

    !------1
    X_Len = 7
    if(allocated(X_inout)) deallocate(X_inout)
    Allocate(X_inout(X_Len))
    a = (/0.34907d0,0.54890d0,0.74776d0,0.94459d0,1.13850d0,1.32850d0,1.51370d0/)
    b = 0d0
    X_inout = CMPLX(a,b)
    !------2
    !------fft 
    X_in(1)=(0.34907D0,0D0)
    X_in(2)=(0.54890D0,0D0)
    X_in(3)=(0.74776D0,0D0)
    X_in(4)=(0.94459D0,0D0)
    X_in(5)=(1.13850D0,0D0)
    X_in(6)=(1.32850D0,0D0)
    X_in(7)=(1.51370D0,0D0)
    do i=1,7
        print*, X_in(i)
    end do
    Status = DftiCreateDescriptor (Desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, Dim, X_Len)
    !Status = DftiGetValue(Desc_handle, DFTI_placement)
    !Status = DftiGetValue(Desc_handle, DFTI_Forward_Scale)
    Status = DftiSetValue(Desc_handle, DFTI_placement, dfti_inplace)
    Status = DftiCommitDescriptor (Desc_handle)
    Status = DftiComputeForward(Desc_handle, X_in)
    Status = DftiFreeDescriptor(Desc_handle)
    if (Status .ne. 0) then
        if (.not. DftiErrorClass(status,DFTI_NO_ERROR)) then
            print *, 'Error: ', Status, DftiErrorMessage(Status)
        endif
    endif
    !------ifft   
    PRINT *,X_in
    X_in(1)=(0.34907D0,0D0)
    X_in(2)=(0.54890D0,0D0)
    X_in(3)=(0.74776D0,0D0)
    X_in(4)=(0.94459D0,0D0)
    X_in(5)=(1.13850D0,0D0)
    X_in(6)=(1.32850D0,0D0)
    X_in(7)=(1.51370D0,0D0)
    do i=1,7
        print*, X_in(i)
    end do 
    Status = DftiCreateDescriptor(Desc_handle, DFTI_DOUBLE, DFTI_COMPLEX, Dim, X_Len)
    Status = DftiSetValue(Desc_handle, DFTI_BACKWARD_SCALE, 1.0d0/X_Len)
    Status = DftiSetValue(Desc_handle, DFTI_placement, dfti_inplace)
    Status = DftiCommitDescriptor (Desc_handle)
    Status = DftiComputeBackward(Desc_handle, X_in)  !720,1
    Status = DftiFreeDescriptor(Desc_handle)
    if (Status .ne. 0) then
        if (.not. DftiErrorClass(status,DFTI_NO_ERROR)) then
            print *, 'Error: ', Status, DftiErrorMessage(Status)
        endif
    endif
    PRINT *,X_in/X_Len



    END PROGRAM


这段是使用MKL的FFT的代码,之前计算过而且FFT IFFT和matlab对比没有错误。
现在使用file:///C:\Users\Davy\AppData\Roaming\Tencent\Users\958256222\QQ\WinTemp\RichOle\$6~3JHLSA}UH__ULQEQ$3FV.png
Error:            5
Intel MKL DFTI ERROR: Descriptor is uncommitted or corrupted

会报这个错,请问是什么原因呢


1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

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

发表于 2021-12-28 08:13:30 | 显示全部楼层
complex*8 全部改成
complex(8)

complex*8 占用8个字节,实部虚部各4,是单精度!!
complex(8) 占用16个字节,实部虚部各8,是双精度

你用了 DFTI_DOUBLE,需要用双精度~!
当前不推荐使用 * 来定义变量的kind值。非常不推荐。

24

帖子

7

主题

0

精华

熟手

F 币
143 元
贡献
88 点
 楼主| 发表于 2021-12-28 13:02:44 | 显示全部楼层
fcode 发表于 2021-12-28 08:13
complex*8 全部改成
complex(8)

谢谢雪球老师!
刘涛最美!

1948

帖子

12

主题

5

精华

论坛跑堂

臭石头雪球

F 币
1298 元
贡献
547 点

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

发表于 2021-12-28 14:09:49 | 显示全部楼层
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-3-29 12:47

Powered by Tencent X3.4

© 2013-2024 Tencent

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