tópico Re: BUG for in-place forward DFTs ?? na Intel® oneAPI Math Kernel Library

BUG for in-place forward DFTs ??

OP1 — Thu, 07 Aug 2008 15:40:24 GMT

The following simple program does not produce the correct results. It's behavior actually seems to be some times correct, most of the times wrong, when you run it several times, consecutively!

It is a simple forward transform of two identical sets of datas (DFTI_NUMBER_OF_TRANSFORMS is > 1) stored consecutively in the same 1D array (TIME_DATA). The two data sets (8 time points each) correspond to a simple cosine signal.

The forward transform for thefirst data set is always properly calculated. The second transform is wrong most of the time.

Can somebody shed some light on this? I have to add that if I do not specify inplace transforms, then it works. Oh, and it seems I just can't use DFTIFREEDESCRIPTOR without having my code crashing. What is going on? I am using IVF 10.1 and the MKL 10.1 Beta (the same behavior was observed with MKL 8.1).

Thanks in advance! This is very puzzling (and frustrating...)

PROGRAM
TEST

USE
MKL_DFTI

IMPLICIT NONE

INTEGER(4)
,PARAMETER :: N_TIME = 8

INTEGER(4)
,PARAMETER :: N_SETS = 2

INTEGER(4)
I,J,S1,S2,S3,S4,S5,S6,S7,ERROR

REAL(8)
PI,T,TIME_DATA(1:2*N_TIME)

TYPE
(DFTI_DESCRIPTOR),POINTER :: DH

DO
J=1,N_SETS

DO I=1,N_TIME

T = (8.0D+0*

DATAN(1.0D+0)*DBLE(I-1)/DBLE(N_TIME))

TIME_DATA((J-1)*N_TIME+I) =

DCOS(T)

ENDDO

ENDDO

WRITE
(*,'(20F5.1)') TIME_DATA

S1 = DFTICREATEDESCRIPTOR(DH,DFTI_DOUBLE,DFTI_REAL,1,N_TIME)

S2 = DFTISETVALUE(DH,DFTI_NUMBER_OF_TRANSFORMS,N_SETS)

S3 = DFTISETVALUE(DH,DFTI_PLACEMENT,DFTI_INPLACE)

S4 = DFTISETVALUE(DH,DFTI_INPUT_DISTANCE,N_TIME)

S5 = DFTISETVALUE(DH,DFTI_OUTPUT_DISTANCE,N_TIME)

S6 = DFTICOMMITDESCRIPTOR(DH)

S7 = DFTICOMPUTEFORWARD(DH,TIME_DATA)

WRITE
(*,'(20F5.1)') TIME_DATA

READ
(*,*)

END
PROGRAM

Re: BUG for in-place forward DFTs ??

Todd_R_Intel — Thu, 07 Aug 2008 21:27:27 GMT

To use CCS storage with the INPLACE algorithm you will need to add padding. That will involve changing N_TIME to (N_TIME/2+1)*2 in several places. Give the following changes a try (- = remove the line, + = add the line):

--- y.f90   2008-08-07 23:45:42.161126300 +0700
+++ z.f90   2008-08-07 23:47:20.551121600 +0700
@@ -4,22 +4,23 @@
 INTEGER(4),PARAMETER :: N_TIME = 8
 INTEGER(4),PARAMETER :: N_SETS = 2
 INTEGER(4) I,J,S1,S2,S3,S4,S5,S6,S7,ERROR
-REAL(8) PI,T,TIME_DATA(1:2*N_TIME)
+REAL(8) PI,T,TIME_DATA(1:2*(N_TIME/2+1)*2)
 TYPE(DFTI_DESCRIPTOR),POINTER :: DH
 DO J=1,N_SETS
 DO I=1,N_TIME
 T = (8.0D+0*DATAN(1.0D+0)*DBLE(I-1)/DBLE(N_TIME))
-TIME_DATA((J-1)*N_TIME+I) = DCOS(T)
+TIME_DATA((J-1)*(N_TIME/2+1)*2+I) = DCOS(T)
 ENDDO
 ENDDO
 WRITE(*,'(20F5.1)') TIME_DATA
 S1 = DFTICREATEDESCRIPTOR(DH,DFTI_DOUBLE,DFTI_REAL,1,N_TIME)
 S2 = DFTISETVALUE(DH,DFTI_NUMBER_OF_TRANSFORMS,N_SETS)
 S3 = DFTISETVALUE(DH,DFTI_PLACEMENT,DFTI_INPLACE)
-S4 = DFTISETVALUE(DH,DFTI_INPUT_DISTANCE,N_TIME)
-S5 = DFTISETVALUE(DH,DFTI_OUTPUT_DISTANCE,N_TIME)
+S4 = DFTISETVALUE(DH,DFTI_INPUT_DISTANCE,(N_TIME/2+1)*2)
+S5 = DFTISETVALUE(DH,DFTI_OUTPUT_DISTANCE,N_TIME/2)
 S6 = DFTICOMMITDESCRIPTOR(DH)
 S7 = DFTICOMPUTEFORWARD(DH,TIME_DATA)
 WRITE(*,'(20F5.1)') TIME_DATA
-READ(*,*)
+S7 = DFTIFREEDESCRIPTOR(DH)
+PRINT *,"End"
 END PROGRAM

Output of the fixed program is:

  1.0  0.7  0.0 -0.7 -1.0 -0.7  0.0  0.7  0.0  0.0  1.0  0.7  0.0 -0.7 -1.0 -0.7  0.0  0.7  0.0  0.0
  0.0  0.0  4.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  4.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 End

Re: BUG for in-place forward DFTs ??

OP1 — Thu, 07 Aug 2008 22:45:22 GMT

Thanks a lot!! This does the trick... It probably also explains the erratic behavior of my code (since the MKL subroutines were manipulating arrays that they thought were longer than their actual size).

This being said... would it be possible to implement a test in DftiCommitDescriptor that would check the compatibility of the arguments set with DftiSetValue and the length of the arrays (when Dfti_Inplace is used)? It seems it's areally easy mistake to do!

Well, thanks again!

Re: BUG for in-place forward DFTs ??

Todd_R_Intel — Fri, 08 Aug 2008 22:17:56 GMT

Yeah, that might be worthwhile. I'll pass your suggestion along. -Todd