BUG for in-place forward DFTs ??

[OP1]

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

[Todd_R_Intel]

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

[OP1]

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!

[Todd_R_Intel]

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