使用英特尔MKL进行3D卷积

我编写了一个C/C++代码，该代码使用Intel MKL来计算具有大约 300×200×200元素的数组的3D卷积。我想应用一个3×3×3或的内核5×5×5。3D输入数组和内核都有真实值。

该3D数组以列方式存储为类型的1D数组double。同样，内核是类型，double并按列保存。例如，

for( int k = 0; k < nk; k++ ) // Loop through the height.
    for( int j = 0; j < nj; j++ ) // Loop through the rows.
        for( int i = 0; i < ni; i++ ) // Loop through the columns.
        {
            ijk = i + ni * j + ni * nj * k;
            my3Darray[ ijk ] = 1.0;
        }

对于卷积的计算，我想对not-in-place输入数组和内核执行FFT并防止它们被修改（我稍后需要在我的代码中使用它们），然后进行向后计算in-place。

当我比较从我的代码获得的结果与MATLAB它们获得的结果时，它们有很大的不同。有人可以帮我解决问题吗？我的代码中缺少什么？

这是MATLAB我使用的代码：

a = ones( 10, 10, 10 );
kernel = ones( 3, 3, 3 );
aconvolved = convn( a, kernel, 'same' );

这是我的C/C++代码：

#include <stdio.h>
#include "mkl.h"

void Conv3D(
    double *in, double *ker, double *out,
    int nRows, int nCols, int nHeights)
{

    int NI = nRows;
    int NJ = nCols;
    int NK = nHeights;

    double *in_fft  = new double [NI*NJ*NK];
    double *ker_fft = new double [NI*NJ*NK];

    DFTI_DESCRIPTOR_HANDLE fft_desc = 0;
    MKL_LONG sizes[]   = { NK, NJ, NI };
    MKL_LONG strides[] = { 0, NJ*NI, NI, 1 };

    DftiCreateDescriptor( &fft_desc, DFTI_DOUBLE, DFTI_REAL, 3, sizes     );
    DftiSetValue        (  fft_desc, DFTI_PLACEMENT     , DFTI_NOT_INPLACE);   // Out-of-place computation.
    DftiSetValue        (  fft_desc, DFTI_INPUT_STRIDES , strides         );
    DftiSetValue        (  fft_desc, DFTI_OUTPUT_STRIDES, strides         );
    DftiSetValue        (  fft_desc, DFTI_BACKWARD_SCALE, 1/NI/NJ/NK      );
    DftiCommitDescriptor(  fft_desc );

    DftiComputeForward  (  fft_desc, in , in_fft  );
    DftiComputeForward  (  fft_desc, ker, ker_fft );

    for (long long i = 0; i < (long long)NI*NJ*NK; ++i )
        out[i] = in_fft[i]*ker_fft[i];

    // In-place computation.
    DftiSetValue        (  fft_desc, DFTI_PLACEMENT, DFTI_INPLACE );
    DftiCommitDescriptor(  fft_desc      );
    DftiComputeBackward (  fft_desc, out );

    DftiFreeDescriptor  ( &fft_desc );

    delete[] in_fft;
    delete[] ker_fft;

}

int main(int argc, char* argv[])
{
    int n = 10;
    int nkernel = 3;

    double *a          = new double [n*n*n]; // This array is real.
    double *aconvolved = new double [n*n*n]; // The convolved array is also real.
    double *kernel     = new double [nkernel*nkernel*nkernel]; // kernel is real.

    // Fill the array with some 'real' numbers.
    for( int i = 0; i < n*n*n; i++ )
        a[ i ] = 1.0;

    // Fill the kernel with some 'real' numbers.
    for( int i = 0; i < nkernel*nkernel*nkernel; i++ )
        kernel[ i ] = 1.0;

    // Calculate the convolution.
    Conv3D( a, kernel, aconvolved, n, n, n );

    printf("Convolved:\n");
    for( int i = 0; i < n*n*n; i++ )
        printf( "%15.8f\n", aconvolved[i] );

    delete[] a;
    delete[] kernel;
    delete[] aconvolved;

    return 0;
}