Skip to content

Commit f46dc9c

Browse files
author
Joan Massich
committed
First commit
0 parents  commit f46dc9c

File tree

3 files changed

+374
-0
lines changed

3 files changed

+374
-0
lines changed

‎.travis.yml

+23
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
sudo: false
2+
language:
3+
- cpp
4+
5+
addons:
6+
apt:
7+
sources:
8+
- ubuntu-toolchain-r-test
9+
10+
before_install:
11+
- wget -q -O /dev/stdout 'https://raw.githubusercontent.com/nemequ/icc-travis/master/install-icc.sh'
12+
- ./install-icc.sh
13+
14+
script:
15+
# - mkdir build
16+
- source ~/.bashrc
17+
- CC=icc CXX=icpc ./configure && make
18+
19+
after_success:
20+
- '[[ ! -z "${INTEL_INSTALL_PATH}" ]] && uninstall_intel_software'
21+
22+
notifications:
23+
email: false

‎README.md

+3
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
# Test project
2+
3+
This project tests compiling against MKL

‎basic_sp_real_dft_3d.c

+348
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,348 @@
1+
/*******************************************************************************
2+
* Copyright 2011-2017 Intel Corporation All Rights Reserved.
3+
*
4+
* The source code, information and material ("Material") contained herein is
5+
* owned by Intel Corporation or its suppliers or licensors, and title to such
6+
* Material remains with Intel Corporation or its suppliers or licensors. The
7+
* Material contains proprietary information of Intel or its suppliers and
8+
* licensors. The Material is protected by worldwide copyright laws and treaty
9+
* provisions. No part of the Material may be used, copied, reproduced,
10+
* modified, published, uploaded, posted, transmitted, distributed or disclosed
11+
* in any way without Intel's prior express written permission. No license under
12+
* any patent, copyright or other intellectual property rights in the Material
13+
* is granted to or conferred upon you, either expressly, by implication,
14+
* inducement, estoppel or otherwise. Any license under such intellectual
15+
* property rights must be express and approved by Intel in writing.
16+
*
17+
* Unless otherwise agreed by Intel in writing, you may not remove or alter this
18+
* notice or any other notice embedded in Materials by Intel or Intel's
19+
* suppliers or licensors in any way.
20+
*******************************************************************************/
21+
22+
/*
23+
! Content:
24+
! A simple example of single-precision real-to-complex in-place 3D
25+
! FFT using Intel(R) MKL DFTI
26+
!
27+
!****************************************************************************/
28+
29+
#include <stdio.h>
30+
#include <stdlib.h>
31+
#include <math.h>
32+
#include <float.h>
33+
#include "mkl_service.h"
34+
#include "mkl_dfti.h"
35+
36+
/* Define the format to printf MKL_LONG values */
37+
#if !defined(MKL_ILP64)
38+
#define LI "%li"
39+
#else
40+
#define LI "%lli"
41+
#endif
42+
43+
static void init_r(float *x, int N1, int N2, int N3, int H1, int H2, int H3);
44+
static void init_c(MKL_Complex8 *x, int N1, int N2, int N3, int H1, int H2,
45+
int H3);
46+
static int verify_r(float *x, int N1, int N2, int N3, int H1, int H2, int H3);
47+
static int verify_c(MKL_Complex8 *x, int N1, int N2, int N3, int H1, int H2,
48+
int H3);
49+
50+
/* Define the format to printf MKL_LONG values */
51+
#if !defined(MKL_ILP64)
52+
#define LI "%li"
53+
#else
54+
#define LI "%lli"
55+
#endif
56+
57+
int main(void)
58+
{
59+
/* Size of 3D transform */
60+
int N1 = 4, N2 = 5, N3 = 13;
61+
62+
/* Arbitrary harmonic used to verify FFT */
63+
int H1 = -1, H2 = 2, H3 = 4;
64+
65+
/* Execution status */
66+
MKL_LONG status = 0;
67+
68+
/* Pointer to input/output data */
69+
float *x = 0;
70+
71+
/* Strides describe data layout in real and conjugate-even domain */
72+
MKL_LONG rs[4], cs[4];
73+
74+
DFTI_DESCRIPTOR_HANDLE hand = 0;
75+
76+
char version[DFTI_VERSION_LENGTH];
77+
78+
DftiGetValue(0, DFTI_VERSION, version);
79+
printf("%s\n", version);
80+
81+
printf("Example basic_sp_real_dft_3d\n");
82+
printf("Forward-Backward single-precision in-place 3D real FFT\n");
83+
printf("Configuration parameters:\n");
84+
printf(" DFTI_PRECISION = DFTI_SINGLE\n");
85+
printf(" DFTI_FORWARD_DOMAIN = DFTI_REAL\n");
86+
printf(" DFTI_DIMENSION = 3\n");
87+
printf(" DFTI_LENGTHS = {%d, %d, %d}\n", N1, N2, N3);
88+
printf(" DFTI_PLACEMENT = DFTI_INPLACE\n");
89+
printf(" DFTI_CONJUGATE_EVEN_STORAGE = DFTI_COMPLEX_COMPLEX\n");
90+
91+
92+
printf("Create DFTI descriptor\n");
93+
{
94+
MKL_LONG N[3]; N[0] = N1; N[1] = N2; N[2] = N3;
95+
status = DftiCreateDescriptor(&hand, DFTI_SINGLE, DFTI_REAL, 3, N);
96+
if (0 != status) goto failed;
97+
}
98+
99+
/* This is not needed, default setting */
100+
/* status = DftiSetValue(hand, DFTI_PLACEMENT, DFTI_INPLACE); */
101+
/* if (0 != status) goto failed; */
102+
103+
104+
printf("Set configuration: CCE storage\n");
105+
status = DftiSetValue(hand, DFTI_CONJUGATE_EVEN_STORAGE,
106+
DFTI_COMPLEX_COMPLEX);
107+
if (0 != status) goto failed;
108+
109+
/* Compute strides */
110+
rs[3] = 1; cs[3] = 1;
111+
rs[2] = (N3/2+1)*2; cs[2] = (N3/2+1);
112+
rs[1] = N2*(N3/2+1)*2; cs[1] = N2*(N3/2+1);
113+
rs[0] = 0; cs[0] = 0;
114+
115+
printf("Set input strides = { "LI", "LI", "LI", "LI" }\n",
116+
rs[0],rs[1],rs[2],rs[3]);
117+
status = DftiSetValue(hand, DFTI_INPUT_STRIDES, rs);
118+
if (0 != status) goto failed;
119+
120+
printf("Set output strides = { "LI", "LI", "LI", "LI" }\n",
121+
cs[0],cs[1],cs[2],cs[3]);
122+
status = DftiSetValue(hand, DFTI_OUTPUT_STRIDES, cs);
123+
if (0 != status) goto failed;
124+
125+
printf("Commit the descriptor\n");
126+
status = DftiCommitDescriptor(hand);
127+
if (0 != status) goto failed;
128+
129+
printf("Allocate data array\n");
130+
x = (float *)mkl_malloc(N1*N2*(N3/2+1)*2*sizeof(float), 64);
131+
if (0 == x) goto failed;
132+
133+
printf("Initialize data for r2c transform\n");
134+
init_r(x, N1, N2, N3, H1, H2, H3);
135+
136+
printf("Compute real-to-complex in-place transform\n");
137+
status = DftiComputeForward(hand, x);
138+
if (0 != status) goto failed;
139+
140+
printf("Verify the result\n");
141+
status = verify_c((MKL_Complex8*)x, N1, N2, N3, H1, H2, H3);
142+
if (0 != status) goto failed;
143+
144+
printf("Change strides to compute backward transform\n");
145+
status = DftiSetValue(hand, DFTI_INPUT_STRIDES, cs);
146+
if (0 != status) goto failed;
147+
status = DftiSetValue(hand, DFTI_OUTPUT_STRIDES, rs);
148+
if (0 != status) goto failed;
149+
150+
printf("Commit the descriptor\n");
151+
status = DftiCommitDescriptor(hand);
152+
if (0 != status) goto failed;
153+
154+
printf("Initialize data for c2r transform\n");
155+
init_c((MKL_Complex8*)x, N1, N2, N3, H1, H2, H3);
156+
157+
printf("Compute backward transform\n");
158+
status = DftiComputeBackward(hand, x);
159+
if (0 != status) goto failed;
160+
161+
printf("Verify the result\n");
162+
status = verify_r(x, N1, N2, N3, H1, H2, H3);
163+
if (0 != status) goto failed;
164+
165+
cleanup:
166+
167+
printf("Free DFTI descriptor\n");
168+
DftiFreeDescriptor(&hand);
169+
170+
printf("Free data array\n");
171+
mkl_free(x);
172+
173+
printf("TEST %s\n",0==status ? "PASSED" : "FAILED");
174+
return status;
175+
176+
failed:
177+
printf(" ERROR, status = "LI"\n", status);
178+
status = 1;
179+
goto cleanup;
180+
}
181+
182+
/* Compute (K*L)%M accurately */
183+
static float moda(int K, int L, int M)
184+
{
185+
return (float)(((long long)K * L) % M);
186+
}
187+
188+
/* Initialize array x(N) to produce unit peaks at x(H) and x(N-H) */
189+
static void init_r(float *x, int N1, int N2, int N3, int H1, int H2, int H3)
190+
{
191+
float TWOPI = 6.2831853071795864769f, phase, factor;
192+
int n1, n2, n3, S1, S2, S3, index;
193+
194+
/* Generalized strides for row-major addressing of x */
195+
S3 = 1;
196+
S2 = (N3/2+1)*2;
197+
S1 = N2*(N3/2+1)*2;
198+
199+
factor = (2*(N1-H1)%N1==0 && 2*(N2-H2)%N2==0 && 2*(N3-H3)%N3==0) ? 1.0f : 2.0f;
200+
for (n1 = 0; n1 < N1; n1++)
201+
{
202+
for (n2 = 0; n2 < N2; n2++)
203+
{
204+
for (n3 = 0; n3 < N3; n3++)
205+
{
206+
phase = moda(n1,H1,N1) / N1;
207+
phase += moda(n2,H2,N2) / N2;
208+
phase += moda(n3,H3,N3) / N3;
209+
index = n1*S1 + n2*S2 + n3*S3;
210+
x[index] = factor * cosf( TWOPI * phase ) / (N1*N2*N3);
211+
}
212+
}
213+
}
214+
}
215+
216+
/* Verify that x has unit peak at H */
217+
static int verify_c(MKL_Complex8 *x,int N1,int N2,int N3,int H1,int H2,int H3)
218+
{
219+
float err, errthr, maxerr;
220+
int n1, n2, n3, S1, S2, S3, index;
221+
222+
/* Generalized strides for row-major addressing of x */
223+
S3 = 1;
224+
S2 = N3/2+1;
225+
S1 = N2*(N3/2+1);
226+
227+
/*
228+
* Note, this simple error bound doesn't take into account error of
229+
* input data
230+
*/
231+
errthr = 2.5f * logf( (float)N1*N2*N3 ) / logf(2.0f) * FLT_EPSILON;
232+
printf(" Check if err is below errthr %.3lg\n", errthr);
233+
234+
maxerr = 0;
235+
for (n1 = 0; n1 < N1; n1++)
236+
{
237+
for (n2 = 0; n2 < N2; n2++)
238+
{
239+
for (n3 = 0; n3 < N3/2+1; n3++)
240+
{
241+
float re_exp = 0.0f, im_exp = 0.0f, re_got, im_got;
242+
243+
if ((( n1-H1)%N1==0 && ( n2-H2)%N2==0 && ( n3-H3)%N3==0) ||
244+
((-n1-H1)%N1==0 && (-n2-H2)%N2==0 && (-n3-H3)%N3==0))
245+
{
246+
re_exp = 1;
247+
}
248+
249+
index = n1*S1 + n2*S2 + n3*S3;
250+
re_got = x[index].real;
251+
im_got = x[index].imag;
252+
err = fabsf(re_got - re_exp) + fabsf(im_got - im_exp);
253+
if (err > maxerr) maxerr = err;
254+
if (!(err < errthr))
255+
{
256+
printf(" x[%i][%i][%i]: ",n1,n2,n3);
257+
printf(" expected (%.7g,%.7g), ",re_exp,im_exp);
258+
printf(" got (%.7g,%.7g), ",re_got,im_got);
259+
printf(" err %.3lg\n", err);
260+
printf(" Verification FAILED\n");
261+
return 1;
262+
}
263+
}
264+
}
265+
}
266+
printf(" Verified, maximum error was %.3lg\n", maxerr);
267+
return 0;
268+
}
269+
270+
/* Initialize array x(N) to produce unit peak at x(H) */
271+
static void init_c(MKL_Complex8 *x,int N1,int N2,int N3,int H1,int H2,int H3)
272+
{
273+
float TWOPI = 6.2831853071795864769f, phase;
274+
int n1, n2, n3, S1, S2, S3, index;
275+
276+
/* Generalized strides for row-major addressing of x */
277+
S3 = 1;
278+
S2 = N3/2+1;
279+
S1 = N2*(N3/2+1);
280+
281+
for (n1 = 0; n1 < N1; n1++)
282+
{
283+
for (n2 = 0; n2 < N2; n2++)
284+
{
285+
for (n3 = 0; n3 < N3/2+1; n3++)
286+
{
287+
phase = moda(n1,H1,N1) / N1;
288+
phase += moda(n2,H2,N2) / N2;
289+
phase += moda(n3,H3,N3) / N3;
290+
index = n1*S1 + n2*S2 + n3*S3;
291+
x[index].real = cosf( TWOPI * phase ) / (N1*N2*N3);
292+
x[index].imag = -sinf( TWOPI * phase ) / (N1*N2*N3);
293+
}
294+
}
295+
}
296+
}
297+
298+
/* Verify that x has unit peak at H */
299+
static int verify_r(float *x, int N1, int N2, int N3, int H1, int H2, int H3)
300+
{
301+
float err, errthr, maxerr;
302+
int n1, n2, n3, S1, S2, S3, index;
303+
304+
/* Generalized strides for row-major addressing of x */
305+
S3 = 1;
306+
S2 = (N3/2+1)*2;
307+
S1 = N2*(N3/2+1)*2;
308+
309+
/*
310+
* Note, this simple error bound doesn't take into account error of
311+
* input data
312+
*/
313+
errthr = 2.5f * logf( (float)N1*N2*N3 ) / logf(2.0f) * FLT_EPSILON;
314+
printf(" Check if err is below errthr %.3lg\n", errthr);
315+
316+
maxerr = 0;
317+
for (n1 = 0; n1 < N1; n1++)
318+
{
319+
for (n2 = 0; n2 < N2; n2++)
320+
{
321+
for (n3 = 0; n3 < N3; n3++)
322+
{
323+
float re_exp = 0.0f, re_got;
324+
325+
if ((n1-H1)%N1==0 && (n2-H2)%N2==0 && (n3-H3)%N3==0)
326+
{
327+
re_exp = 1;
328+
}
329+
330+
index = n1*S1 + n2*S2 + n3*S3;
331+
re_got = x[index];
332+
err = fabsf(re_got - re_exp);
333+
if (err > maxerr) maxerr = err;
334+
if (!(err < errthr))
335+
{
336+
printf(" x[%i][%i][%i]: ",n1,n2,n3);
337+
printf(" expected %.7g, ",re_exp);
338+
printf(" got %.7g, ",re_got);
339+
printf(" err %.3lg\n", err);
340+
printf(" Verification FAILED\n");
341+
return 1;
342+
}
343+
}
344+
}
345+
}
346+
printf(" Verified, maximum error was %.3lg\n", maxerr);
347+
return 0;
348+
}

0 commit comments

Comments
 (0)