?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed include statements such that double precision version of genborn.c
[gromacs/rigid-bodies.git] / src / mdlib / gmx_fft_mkl.c
blob8f4cd4d0dbca7c5656045f3df29f27cb88596236
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * Gromacs 4.0 Copyright (c) 1991-2003
5 * David van der Spoel, Erik Lindahl, University of Groningen.
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
14 *
15 * And Hey:
16 * Gnomes, ROck Monsters And Chili Sauce
17 */
18 #ifdef HAVE_CONFIG_H
19 #include <config.h>
20 #endif
21
22 #ifdef GMX_FFT_MKL
23
24 #include <errno.h>
25 #include <stdlib.h>
26
27 #include <mkl_dfti.h>
28
29
30 #include "gmx_fft.h"
31 #include "gmx_fatal.h"
32
33
34 /* For MKL version (<10.0), we should define MKL_LONG. */
35 #ifndef MKL_LONG
36 #define MKL_LONG long int
37 #endif
38
39
40 #ifdef GMX_DOUBLE
41 #define GMX_DFTI_PREC DFTI_DOUBLE
42 #else
43 #define GMX_DFTI_PREC DFTI_SINGLE
44 #endif
45
46 /* Contents of the Intel MKL FFT fft datatype.
47 *
48 * Note that this is one of several possible implementations of gmx_fft_t.
49 *
50 * The MKL _API_ supports 1D,2D, and 3D transforms, including real-to-complex.
51 * Unfortunately the actual library implementation does not support 3D real
52 * transforms as of version 7.2, and versions before 7.0 don't support 2D real
53 * either. In addition, the multi-dimensional storage format for real data
54 * is not compatible with our padding.
55 *
56 * To work around this we roll our own 2D and 3D real-to-complex transforms,
57 * using separate X/Y/Z handles defined to perform (ny*nz), (nx*nz), and
58 * (nx*ny) transforms at once when necessary. To perform strided multiple
59 * transforms out-of-place (i.e., without padding in the last dimension)
60 * on the fly we also need to separate the forward and backward
61 * handles for real-to-complex/complex-to-real data permutation.
62 *
63 * This makes it necessary to define 3 handles for in-place FFTs, and 4 for
64 * the out-of-place transforms. Still, whenever possible we try to use
65 * a single 3D-transform handle instead.
66 *
67 * So, the handles are enumerated as follows:
68 *
69 * 1D FFT (real too): Index 0 is the handle for the entire FFT
70 * 2D complex FFT: Index 0 is the handle for the entire FFT
71 * 3D complex FFT: Index 0 is the handle for the entire FFT
72 * 2D, inplace real FFT: 0=FFTx, 1=FFTy handle
73 * 2D, ooplace real FFT: 0=FFTx, 1=real-to-complex FFTy, 2=complex-to-real FFTy
74 * 3D, inplace real FFT: 0=FFTx, 1=FFTy, 2=FFTz handle
75 * 3D, ooplace real FFT: 0=FFTx, 1=FFTy, 2=r2c FFTz, 3=c2r FFTz
76 *
77 * Intel people reading this: Learn from FFTW what a good interface looks like :-)
78 *
79 */
80 struct gmx_fft
81 {
82 int ndim; /**< Number of dimensions in FFT */
83 int nx; /**< Length of X transform */
84 int ny; /**< Length of Y transform */
85 int nz; /**< Length of Z transform */
86 int real_fft; /**< 1 if real FFT, otherwise 0 */
87 DFTI_DESCRIPTOR * inplace[3]; /**< in-place FFT */
88 DFTI_DESCRIPTOR * ooplace[4]; /**< out-of-place FFT */
89 t_complex * work; /**< Enable out-of-place c2r FFT */
90 };
91
92
93
94 int
95 gmx_fft_init_1d(gmx_fft_t * pfft,
96 int nx,
97 gmx_fft_flag flags)
98 {
99 gmx_fft_t fft;
100 int d;
101 int status;
102
103 if(pfft==NULL)
104 {
105 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
106 return EINVAL;
107 }
108 *pfft = NULL;
109
110 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
111 {
112 return ENOMEM;
113 }
114
115 /* Mark all handles invalid */
116 for(d=0;d<3;d++)
117 {
118 fft->inplace[d] = fft->ooplace[d] = NULL;
119 }
120 fft->ooplace[3] = NULL;
121
122
123 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
124
125 if( status == 0 )
126 status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
127
128 if( status == 0 )
129 status = DftiCommitDescriptor(fft->inplace[0]);
130
131
132 if( status == 0 )
133 status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
134
135 if( status == 0)
136 DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
137
138 if( status == 0)
139 DftiCommitDescriptor(fft->ooplace[0]);
140
141
142 if( status != 0 )
143 {
144 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
145 gmx_fft_destroy(fft);
146 return status;
147 }
148
149 fft->ndim = 1;
150 fft->nx = nx;
151 fft->real_fft = 0;
152 fft->work = NULL;
153
154 *pfft = fft;
155 return 0;
156 }
157
158
159
160 int
161 gmx_fft_init_1d_real(gmx_fft_t * pfft,
162 int nx,
163 gmx_fft_flag flags)
164 {
165 gmx_fft_t fft;
166 int d;
167 int status;
168
169 if(pfft==NULL)
170 {
171 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
172 return EINVAL;
173 }
174 *pfft = NULL;
175
176 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
177 {
178 return ENOMEM;
179 }
180
181 /* Mark all handles invalid */
182 for(d=0;d<3;d++)
183 {
184 fft->inplace[d] = fft->ooplace[d] = NULL;
185 }
186 fft->ooplace[3] = NULL;
187
188 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nx);
189
190 if( status == 0 )
191 status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
192
193 if( status == 0 )
194 status = DftiCommitDescriptor(fft->inplace[0]);
195
196
197 if( status == 0 )
198 status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nx);
199
200 if( status == 0 )
201 status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
202
203 if( status == 0 )
204 status = DftiCommitDescriptor(fft->ooplace[0]);
205
206
207 if(status == DFTI_UNIMPLEMENTED)
208 {
209 gmx_fatal(FARGS,
210 "The linked Intel MKL version (<6.0?) cannot do real FFTs.");
211 gmx_fft_destroy(fft);
212 return status;
213 }
214
215
216 if( status != 0 )
217 {
218 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
219 gmx_fft_destroy(fft);
220 return status;
221 }
222
223 fft->ndim = 1;
224 fft->nx = nx;
225 fft->real_fft = 1;
226 fft->work = NULL;
227
228 *pfft = fft;
229 return 0;
230 }
231
232
233
234 int
235 gmx_fft_init_2d(gmx_fft_t * pfft,
236 int nx,
237 int ny,
238 gmx_fft_flag flags)
239 {
240 gmx_fft_t fft;
241 int d;
242 int status;
243 MKL_LONG length[2];
244
245 if(pfft==NULL)
246 {
247 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
248 return EINVAL;
249 }
250 *pfft = NULL;
251
252 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
253 {
254 return ENOMEM;
255 }
256
257 /* Mark all handles invalid */
258 for(d=0;d<3;d++)
259 {
260 fft->inplace[d] = fft->ooplace[d] = NULL;
261 }
262 fft->ooplace[3] = NULL;
263
264 length[0] = nx;
265 length[1] = ny;
266
267 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,2,length);
268
269 if( status == 0 )
270 status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
271
272 if( status == 0 )
273 status = DftiCommitDescriptor(fft->inplace[0]);
274
275
276 if( status == 0 )
277 status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,2,length);
278
279 if( status == 0 )
280 status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
281
282 if( status == 0 )
283 status = DftiCommitDescriptor(fft->ooplace[0]);
284
285
286 if( status != 0 )
287 {
288 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
289 gmx_fft_destroy(fft);
290 return status;
291 }
292
293 fft->ndim = 2;
294 fft->nx = nx;
295 fft->ny = ny;
296 fft->real_fft = 0;
297 fft->work = NULL;
298
299 *pfft = fft;
300 return 0;
301 }
302
303
304
305 int
306 gmx_fft_init_2d_real(gmx_fft_t * pfft,
307 int nx,
308 int ny,
309 gmx_fft_flag flags)
310 {
311 gmx_fft_t fft;
312 int d;
313 int status;
314 MKL_LONG stride[2];
315 MKL_LONG nyc;
316
317 if(pfft==NULL)
318 {
319 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
320 return EINVAL;
321 }
322 *pfft = NULL;
323
324 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
325 {
326 return ENOMEM;
327 }
328
329 nyc = (ny/2 + 1);
330
331 /* Mark all handles invalid */
332 for(d=0;d<3;d++)
333 {
334 fft->inplace[d] = fft->ooplace[d] = NULL;
335 }
336 fft->ooplace[3] = NULL;
337
338 /* Roll our own 2D real transform using multiple transforms in MKL,
339 * since the current MKL versions does not support our storage format,
340 * and all but the most recent don't even have 2D real FFTs.
341 */
342
343 /* In-place X FFT */
344 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
345
346 if ( status == 0 )
347 {
348 stride[0] = 0;
349 stride[1] = nyc;
350
351 status =
352 (DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE) ||
353 DftiSetValue(fft->inplace[0],DFTI_NUMBER_OF_TRANSFORMS,nyc) ||
354 DftiSetValue(fft->inplace[0],DFTI_INPUT_DISTANCE,1) ||
355 DftiSetValue(fft->inplace[0],DFTI_INPUT_STRIDES,stride) ||
356 DftiSetValue(fft->inplace[0],DFTI_OUTPUT_DISTANCE,1) ||
357 DftiSetValue(fft->inplace[0],DFTI_OUTPUT_STRIDES,stride));
358 }
359
360 if( status == 0 )
361 status = DftiCommitDescriptor(fft->inplace[0]);
362
363 /* Out-of-place X FFT */
364 if( status == 0 )
365 status = DftiCreateDescriptor(&(fft->ooplace[0]),GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
366
367 if( status == 0 )
368 {
369 stride[0] = 0;
370 stride[1] = nyc;
371
372 status =
373 (DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
374 DftiSetValue(fft->ooplace[0],DFTI_NUMBER_OF_TRANSFORMS,nyc) ||
375 DftiSetValue(fft->ooplace[0],DFTI_INPUT_DISTANCE,1) ||
376 DftiSetValue(fft->ooplace[0],DFTI_INPUT_STRIDES,stride) ||
377 DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_DISTANCE,1) ||
378 DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_STRIDES,stride));
379 }
380
381 if( status == 0 )
382 status = DftiCommitDescriptor(fft->ooplace[0]);
383
384
385 /* In-place Y FFT */
386 if( status == 0 )
387 status = DftiCreateDescriptor(&fft->inplace[1],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
388
389 if( status == 0 )
390 {
391 stride[0] = 0;
392 stride[1] = 1;
393
394 status =
395 (DftiSetValue(fft->inplace[1],DFTI_PLACEMENT,DFTI_INPLACE) ||
396 DftiSetValue(fft->inplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx) ||
397 DftiSetValue(fft->inplace[1],DFTI_INPUT_DISTANCE,2*nyc) ||
398 DftiSetValue(fft->inplace[1],DFTI_INPUT_STRIDES,stride) ||
399 DftiSetValue(fft->inplace[1],DFTI_OUTPUT_DISTANCE,2*nyc) ||
400 DftiSetValue(fft->inplace[1],DFTI_OUTPUT_STRIDES,stride) ||
401 DftiCommitDescriptor(fft->inplace[1]));
402 }
403
404
405 /* Out-of-place real-to-complex (affects output distance) Y FFT */
406 if( status == 0 )
407 status = DftiCreateDescriptor(&fft->ooplace[1],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
408
409 if( status == 0 )
410 {
411 stride[0] = 0;
412 stride[1] = 1;
413
414 status =
415 (DftiSetValue(fft->ooplace[1],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
416 DftiSetValue(fft->ooplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx) ||
417 DftiSetValue(fft->ooplace[1],DFTI_INPUT_DISTANCE,(MKL_LONG)ny) ||
418 DftiSetValue(fft->ooplace[1],DFTI_INPUT_STRIDES,stride) ||
419 DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_DISTANCE,2*nyc) ||
420 DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_STRIDES,stride) ||
421 DftiCommitDescriptor(fft->ooplace[1]));
422 }
423
424
425 /* Out-of-place complex-to-real (affects output distance) Y FFT */
426 if( status == 0 )
427 status = DftiCreateDescriptor(&fft->ooplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)ny);
428
429 if( status == 0 )
430 {
431 stride[0] = 0;
432 stride[1] = 1;
433
434 status =
435 (DftiSetValue(fft->ooplace[2],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
436 DftiSetValue(fft->ooplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx) ||
437 DftiSetValue(fft->ooplace[2],DFTI_INPUT_DISTANCE,2*nyc) ||
438 DftiSetValue(fft->ooplace[2],DFTI_INPUT_STRIDES,stride) ||
439 DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)ny) ||
440 DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_STRIDES,stride) ||
441 DftiCommitDescriptor(fft->ooplace[2]));
442 }
443
444
445 if ( status == 0 )
446 {
447 if ((fft->work = (t_complex *)malloc(sizeof(t_complex)*(nx*(ny/2+1)))) == NULL)
448 {
449 status = ENOMEM;
450 }
451 }
452
453 if( status != 0 )
454 {
455 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
456 gmx_fft_destroy(fft);
457 return status;
458 }
459
460 fft->ndim = 2;
461 fft->nx = nx;
462 fft->ny = ny;
463 fft->real_fft = 1;
464
465 *pfft = fft;
466 return 0;
467 }
468
469
470
471 int
472 gmx_fft_init_3d(gmx_fft_t * pfft,
473 int nx,
474 int ny,
475 int nz,
476 gmx_fft_flag flags)
477 {
478 gmx_fft_t fft;
479 int d;
480 MKL_LONG length[3];
481 int status;
482
483 if(pfft==NULL)
484 {
485 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
486 return EINVAL;
487 }
488 *pfft = NULL;
489
490 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
491 {
492 return ENOMEM;
493 }
494
495 /* Mark all handles invalid */
496 for(d=0;d<3;d++)
497 {
498 fft->inplace[d] = fft->ooplace[d] = NULL;
499 }
500 fft->ooplace[3] = NULL;
501
502 length[0] = nx;
503 length[1] = ny;
504 length[2] = nz;
505
506 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,(MKL_LONG)3,length);
507
508 if( status == 0 )
509 status = DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE);
510
511 if( status == 0 )
512 status = DftiCommitDescriptor(fft->inplace[0]);
513
514
515 if( status == 0 )
516 status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,(MKL_LONG)3,length);
517
518 if( status == 0 )
519 status = DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE);
520
521 if( status == 0 )
522 status = DftiCommitDescriptor(fft->ooplace[0]);
523
524
525 if( status != 0 )
526 {
527 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
528 gmx_fft_destroy(fft);
529 return status;
530 }
531
532
533 fft->ndim = 3;
534 fft->nx = nx;
535 fft->ny = ny;
536 fft->nz = nz;
537 fft->real_fft = 0;
538 fft->work = NULL;
539
540 *pfft = fft;
541 return 0;
542 }
543
544
545
546
547 int
548 gmx_fft_init_3d_real(gmx_fft_t * pfft,
549 int nx,
550 int ny,
551 int nz,
552 gmx_fft_flag flags)
553 {
554 gmx_fft_t fft;
555 int d;
556 int status;
557 MKL_LONG stride[2];
558 int nzc;
559
560 if(pfft==NULL)
561 {
562 gmx_fatal(FARGS,"Invalid opaque FFT datatype pointer.");
563 return EINVAL;
564 }
565 *pfft = NULL;
566
567 nzc = (nz/2 + 1);
568
569 if( (fft = (gmx_fft_t)malloc(sizeof(struct gmx_fft))) == NULL)
570 {
571 return ENOMEM;
572 }
573
574 /* Mark all handles invalid */
575 for(d=0;d<3;d++)
576 {
577 fft->inplace[d] = fft->ooplace[d] = NULL;
578 }
579 fft->ooplace[3] = NULL;
580
581 /* Roll our own 3D real transform using multiple transforms in MKL,
582 * since the current MKL versions does not support our storage format
583 * or 3D real transforms.
584 */
585
586 /* In-place X FFT.
587 * ny*nzc complex-to-complex transforms, length nx
588 * transform distance: 1
589 * element strides: ny*nzc
590 */
591 status = DftiCreateDescriptor(&fft->inplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
592
593 if ( status == 0)
594 {
595 stride[0] = 0;
596 stride[1] = ny*nzc;
597
598 status =
599 (DftiSetValue(fft->inplace[0],DFTI_PLACEMENT,DFTI_INPLACE) ||
600 DftiSetValue(fft->inplace[0],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)ny*nzc) ||
601 DftiSetValue(fft->inplace[0],DFTI_INPUT_DISTANCE,1) ||
602 DftiSetValue(fft->inplace[0],DFTI_INPUT_STRIDES,stride) ||
603 DftiSetValue(fft->inplace[0],DFTI_OUTPUT_DISTANCE,1) ||
604 DftiSetValue(fft->inplace[0],DFTI_OUTPUT_STRIDES,stride) ||
605 DftiCommitDescriptor(fft->inplace[0]));
606 }
607
608 /* Out-of-place X FFT:
609 * ny*nzc complex-to-complex transforms, length nx
610 * transform distance: 1
611 * element strides: ny*nzc
612 */
613 if( status == 0 )
614 status = DftiCreateDescriptor(&fft->ooplace[0],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)nx);
615
616 if( status == 0 )
617 {
618 stride[0] = 0;
619 stride[1] = ny*nzc;
620
621 status =
622 (DftiSetValue(fft->ooplace[0],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
623 DftiSetValue(fft->ooplace[0],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)ny*nzc) ||
624 DftiSetValue(fft->ooplace[0],DFTI_INPUT_DISTANCE,1) ||
625 DftiSetValue(fft->ooplace[0],DFTI_INPUT_STRIDES,stride) ||
626 DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_DISTANCE,1) ||
627 DftiSetValue(fft->ooplace[0],DFTI_OUTPUT_STRIDES,stride) ||
628 DftiCommitDescriptor(fft->ooplace[0]));
629 }
630
631
632 /* In-place Y FFT.
633 * We cannot do all NX*NZC transforms at once, so define a handle to do
634 * NZC transforms, and then execute it NX times.
635 * nzc complex-to-complex transforms, length ny
636 * transform distance: 1
637 * element strides: nzc
638 */
639 if( status == 0 )
640 status = DftiCreateDescriptor(&fft->inplace[1],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)ny);
641
642 if( status == 0 )
643 {
644 stride[0] = 0;
645 stride[1] = nzc;
646
647 status =
648 (DftiSetValue(fft->inplace[1],DFTI_PLACEMENT,DFTI_INPLACE) ||
649 DftiSetValue(fft->inplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nzc) ||
650 DftiSetValue(fft->inplace[1],DFTI_INPUT_DISTANCE,1) ||
651 DftiSetValue(fft->inplace[1],DFTI_INPUT_STRIDES,stride) ||
652 DftiSetValue(fft->inplace[1],DFTI_OUTPUT_DISTANCE,1) ||
653 DftiSetValue(fft->inplace[1],DFTI_OUTPUT_STRIDES,stride) ||
654 DftiCommitDescriptor(fft->inplace[1]));
655 }
656
657
658 /* Out-of-place Y FFT:
659 * We cannot do all NX*NZC transforms at once, so define a handle to do
660 * NZC transforms, and then execute it NX times.
661 * nzc complex-to-complex transforms, length ny
662 * transform distance: 1
663 * element strides: nzc
664 */
665 if( status == 0 )
666 status = DftiCreateDescriptor(&fft->ooplace[1],GMX_DFTI_PREC,DFTI_COMPLEX,1,(MKL_LONG)ny);
667
668 if( status == 0 )
669 {
670 stride[0] = 0;
671 stride[1] = nzc;
672
673 status =
674 (DftiSetValue(fft->ooplace[1],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
675 DftiSetValue(fft->ooplace[1],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nzc) ||
676 DftiSetValue(fft->ooplace[1],DFTI_INPUT_DISTANCE,1) ||
677 DftiSetValue(fft->ooplace[1],DFTI_INPUT_STRIDES,stride) ||
678 DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_DISTANCE,1) ||
679 DftiSetValue(fft->ooplace[1],DFTI_OUTPUT_STRIDES,stride) ||
680 DftiCommitDescriptor(fft->ooplace[1]));
681 }
682
683 /* In-place Z FFT:
684 * nx*ny real-to-complex transforms, length nz
685 * transform distance: nzc*2 -> nzc*2
686 * element strides: 1
687 */
688 if( status == 0 )
689 status = DftiCreateDescriptor(&fft->inplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
690
691 if( status == 0 )
692 {
693 stride[0] = 0;
694 stride[1] = 1;
695
696 status =
697 (DftiSetValue(fft->inplace[2],DFTI_PLACEMENT,DFTI_INPLACE) ||
698 DftiSetValue(fft->inplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
699 DftiSetValue(fft->inplace[2],DFTI_INPUT_DISTANCE,(MKL_LONG)nzc*2) ||
700 DftiSetValue(fft->inplace[2],DFTI_INPUT_STRIDES,stride) ||
701 DftiSetValue(fft->inplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nzc*2) ||
702 DftiSetValue(fft->inplace[2],DFTI_OUTPUT_STRIDES,stride) ||
703 DftiCommitDescriptor(fft->inplace[2]));
704 }
705
706
707 /* Out-of-place real-to-complex (affects distance) Z FFT:
708 * nx*ny real-to-complex transforms, length nz
709 * transform distance: nz -> nzc*2
710 * element STRIDES: 1
711 */
712 if( status == 0 )
713 status = DftiCreateDescriptor(&fft->ooplace[2],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
714
715 if( status == 0 )
716 {
717 stride[0] = 0;
718 stride[1] = 1;
719
720 status =
721 (DftiSetValue(fft->ooplace[2],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
722 DftiSetValue(fft->ooplace[2],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
723 DftiSetValue(fft->ooplace[2],DFTI_INPUT_DISTANCE,(MKL_LONG)nz) ||
724 DftiSetValue(fft->ooplace[2],DFTI_INPUT_STRIDES,stride) ||
725 DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nzc*2) ||
726 DftiSetValue(fft->ooplace[2],DFTI_OUTPUT_STRIDES,stride) ||
727 DftiCommitDescriptor(fft->ooplace[2]));
728 }
729
730
731 /* Out-of-place complex-to-real (affects distance) Z FFT:
732 * nx*ny real-to-complex transforms, length nz
733 * transform distance: nzc*2 -> nz
734 * element STRIDES: 1
735 */
736 if( status == 0 )
737 status = DftiCreateDescriptor(&fft->ooplace[3],GMX_DFTI_PREC,DFTI_REAL,1,(MKL_LONG)nz);
738
739 if( status == 0 )
740 {
741 stride[0] = 0;
742 stride[1] = 1;
743
744 status =
745 (DftiSetValue(fft->ooplace[3],DFTI_PLACEMENT,DFTI_NOT_INPLACE) ||
746 DftiSetValue(fft->ooplace[3],DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)nx*ny) ||
747 DftiSetValue(fft->ooplace[3],DFTI_INPUT_DISTANCE,(MKL_LONG)nzc*2) ||
748 DftiSetValue(fft->ooplace[3],DFTI_INPUT_STRIDES,stride) ||
749 DftiSetValue(fft->ooplace[3],DFTI_OUTPUT_DISTANCE,(MKL_LONG)nz) ||
750 DftiSetValue(fft->ooplace[3],DFTI_OUTPUT_STRIDES,stride) ||
751 DftiCommitDescriptor(fft->ooplace[3]));
752 }
753
754
755 if ( status == 0 )
756 {
757 if ((fft->work = (t_complex *)malloc(sizeof(t_complex)*(nx*ny*(nz/2+1)))) == NULL)
758 {
759 status = ENOMEM;
760 }
761 }
762
763
764 if( status != 0 )
765 {
766 gmx_fatal(FARGS,"Error initializing Intel MKL FFT; status=%d",status);
767 gmx_fft_destroy(fft);
768 return status;
769 }
770
771
772 fft->ndim = 3;
773 fft->nx = nx;
774 fft->ny = ny;
775 fft->nz = nz;
776 fft->real_fft = 1;
777
778 *pfft = fft;
779 return 0;
780 }
781
782
783
784
785 int
786 gmx_fft_1d(gmx_fft_t fft,
787 enum gmx_fft_direction dir,
788 void * in_data,
789 void * out_data)
790 {
791 int inplace = (in_data == out_data);
792 int status = 0;
793
794 if( (fft->real_fft == 1) || (fft->ndim != 1) ||
795 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
796 {
797 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
798 return EINVAL;
799 }
800
801 if(dir==GMX_FFT_FORWARD)
802 {
803 if(inplace)
804 {
805 status = DftiComputeForward(fft->inplace[0],in_data);
806 }
807 else
808 {
809 status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
810 }
811 }
812 else
813 {
814 if(inplace)
815 {
816 status = DftiComputeBackward(fft->inplace[0],in_data);
817 }
818 else
819 {
820 status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
821 }
822 }
823
824 if( status != 0 )
825 {
826 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
827 status = -1;
828 }
829
830 return status;
831 }
832
833
834
835 int
836 gmx_fft_1d_real(gmx_fft_t fft,
837 enum gmx_fft_direction dir,
838 void * in_data,
839 void * out_data)
840 {
841 int inplace = (in_data == out_data);
842 int status = 0;
843
844 if( (fft->real_fft != 1) || (fft->ndim != 1) ||
845 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
846 {
847 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
848 return EINVAL;
849 }
850
851 if(dir==GMX_FFT_REAL_TO_COMPLEX)
852 {
853 if(inplace)
854 {
855 status = DftiComputeForward(fft->inplace[0],in_data);
856 }
857 else
858 {
859 status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
860 }
861 }
862 else
863 {
864 if(inplace)
865 {
866 status = DftiComputeBackward(fft->inplace[0],in_data);
867 }
868 else
869 {
870 status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
871 }
872 }
873
874 if( status != 0 )
875 {
876 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
877 status = -1;
878 }
879
880 return status;
881 }
882
883
884 int
885 gmx_fft_2d(gmx_fft_t fft,
886 enum gmx_fft_direction dir,
887 void * in_data,
888 void * out_data)
889 {
890 int inplace = (in_data == out_data);
891 int status = 0;
892
893 if( (fft->real_fft == 1) || (fft->ndim != 2) ||
894 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
895 {
896 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
897 return EINVAL;
898 }
899
900 if(dir==GMX_FFT_FORWARD)
901 {
902 if(inplace)
903 {
904 status = DftiComputeForward(fft->inplace[0],in_data);
905 }
906 else
907 {
908 status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
909 }
910 }
911 else
912 {
913 if(inplace)
914 {
915 status = DftiComputeBackward(fft->inplace[0],in_data);
916 }
917 else
918 {
919 status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
920 }
921 }
922
923 if( status != 0 )
924 {
925 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
926 status = -1;
927 }
928
929 return status;
930 }
931
932
933 int
934 gmx_fft_2d_real(gmx_fft_t fft,
935 enum gmx_fft_direction dir,
936 void * in_data,
937 void * out_data)
938 {
939 int inplace = (in_data == out_data);
940 int status = 0;
941
942 if( (fft->real_fft != 1) || (fft->ndim != 2) ||
943 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
944 {
945 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
946 return EINVAL;
947 }
948
949 if(dir==GMX_FFT_REAL_TO_COMPLEX)
950 {
951 if(inplace)
952 {
953 /* real-to-complex in Y dimension, in-place */
954 status = DftiComputeForward(fft->inplace[1],in_data);
955
956 /* complex-to-complex in X dimension, in-place */
957 if ( status == 0 )
958 status = DftiComputeForward(fft->inplace[0],in_data);
959 }
960 else
961 {
962 /* real-to-complex in Y dimension, in_data to out_data */
963 status = DftiComputeForward(fft->ooplace[1],in_data,out_data);
964
965 /* complex-to-complex in X dimension, in-place to out_data */
966 if ( status == 0 )
967 status = DftiComputeForward(fft->inplace[0],out_data);
968 }
969 }
970 else
971 {
972 if(inplace)
973 {
974 /* complex-to-complex in X dimension, in-place */
975 status = DftiComputeBackward(fft->inplace[0],in_data);
976
977 /* complex-to-real in Y dimension, in-place */
978 if ( status == 0 )
979 status = DftiComputeBackward(fft->inplace[1],in_data);
980
981 }
982 else
983 {
984 /* complex-to-complex in X dimension, from in_data to work */
985 status = DftiComputeBackward(fft->ooplace[0],in_data,fft->work);
986
987 /* complex-to-real in Y dimension, from work to out_data */
988 if ( status == 0 )
989 status = DftiComputeBackward(fft->ooplace[1],fft->work,out_data);
990
991 }
992 }
993
994 if( status != 0 )
995 {
996 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
997 status = -1;
998 }
999
1000 return status;
1001 }
1002
1003
1004 int
1005 gmx_fft_3d(gmx_fft_t fft,
1006 enum gmx_fft_direction dir,
1007 void * in_data,
1008 void * out_data)
1009 {
1010 int inplace = (in_data == out_data);
1011 int status = 0;
1012
1013 if( (fft->real_fft == 1) || (fft->ndim != 3) ||
1014 ((dir != GMX_FFT_FORWARD) && (dir != GMX_FFT_BACKWARD)) )
1015 {
1016 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1017 return EINVAL;
1018 }
1019
1020 if(dir==GMX_FFT_FORWARD)
1021 {
1022 if(inplace)
1023 {
1024 status = DftiComputeForward(fft->inplace[0],in_data);
1025 }
1026 else
1027 {
1028 status = DftiComputeForward(fft->ooplace[0],in_data,out_data);
1029 }
1030 }
1031 else
1032 {
1033 if(inplace)
1034 {
1035 status = DftiComputeBackward(fft->inplace[0],in_data);
1036 }
1037 else
1038 {
1039 status = DftiComputeBackward(fft->ooplace[0],in_data,out_data);
1040 }
1041 }
1042
1043 if( status != 0 )
1044 {
1045 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
1046 status = -1;
1047 }
1048
1049 return status;
1050 }
1051
1052
1053 int
1054 gmx_fft_3d_real(gmx_fft_t fft,
1055 enum gmx_fft_direction dir,
1056 void * in_data,
1057 void * out_data)
1058 {
1059 int inplace = (in_data == out_data);
1060 int status = 0;
1061 int i;
1062 int nx,ny,nzc;
1063
1064 nx = fft->nx;
1065 ny = fft->ny;
1066 nzc = fft->nz/2 + 1;
1067
1068 if( (fft->real_fft != 1) || (fft->ndim != 3) ||
1069 ((dir != GMX_FFT_REAL_TO_COMPLEX) && (dir != GMX_FFT_COMPLEX_TO_REAL)) )
1070 {
1071 gmx_fatal(FARGS,"FFT plan mismatch - bad plan or direction.");
1072 return EINVAL;
1073 }
1074
1075 if(dir==GMX_FFT_REAL_TO_COMPLEX)
1076 {
1077 if(inplace)
1078 {
1079 /* real-to-complex in Z dimension, in-place */
1080 status = DftiComputeForward(fft->inplace[2],in_data);
1081
1082 /* complex-to-complex in Y dimension, in-place */
1083 for(i=0;i<nx;i++)
1084 {
1085 if ( status == 0 )
1086 status = DftiComputeForward(fft->inplace[1],(t_complex *)in_data+i*ny*nzc);
1087 }
1088
1089 /* complex-to-complex in X dimension, in-place */
1090 if ( status == 0 )
1091 status = DftiComputeForward(fft->inplace[0],in_data);
1092 }
1093 else
1094 {
1095 /* real-to-complex in Z dimension, from in_data to out_data */
1096 status = DftiComputeForward(fft->ooplace[2],in_data,out_data);
1097
1098 /* complex-to-complex in Y dimension, in-place */
1099 for(i=0;i<nx;i++)
1100 {
1101 if ( status == 0 )
1102 status = DftiComputeForward(fft->inplace[1],(t_complex *)out_data+i*ny*nzc);
1103 }
1104
1105 /* complex-to-complex in X dimension, in-place */
1106 if ( status == 0 )
1107 status = DftiComputeForward(fft->inplace[0],out_data);
1108 }
1109 }
1110 else
1111 {
1112 if(inplace)
1113 {
1114 /* complex-to-complex in X dimension, in-place */
1115 status = DftiComputeBackward(fft->inplace[0],in_data);
1116
1117 /* complex-to-complex in Y dimension, in-place */
1118 for(i=0;i<nx;i++)
1119 {
1120 if ( status == 0 )
1121 status = DftiComputeBackward(fft->inplace[1],(t_complex *)in_data+i*ny*nzc);
1122 }
1123
1124 /* complex-to-real in Z dimension, in-place */
1125 if ( status == 0 )
1126 status = DftiComputeBackward(fft->inplace[2],in_data);
1127 }
1128 else
1129 {
1130 /* complex-to-complex in X dimension, from in_data to work */
1131 status = DftiComputeBackward(fft->ooplace[0],in_data,fft->work);
1132
1133 /* complex-to-complex in Y dimension, in-place */
1134 for(i=0;i<nx;i++)
1135 {
1136 if ( status == 0 )
1137 status = DftiComputeBackward(fft->inplace[1],fft->work+i*ny*nzc);
1138 }
1139
1140 /* complex-to-real in Z dimension, work to out_data */
1141 if ( status == 0 )
1142 status = DftiComputeBackward(fft->ooplace[2],fft->work,out_data);
1143 }
1144 }
1145
1146 if( status != 0 )
1147 {
1148 gmx_fatal(FARGS,"Error executing Intel MKL FFT.");
1149 status = -1;
1150 }
1151
1152 return status;
1153 }
1154
1155
1156
1157 void
1158 gmx_fft_destroy(gmx_fft_t fft)
1159 {
1160 int d;
1161
1162 if(fft != NULL)
1163 {
1164 for(d=0;d<3;d++)
1165 {
1166 if(fft->inplace[d] != NULL)
1167 {
1168 DftiFreeDescriptor(&fft->inplace[d]);
1169 }
1170 if(fft->ooplace[d] != NULL)
1171 {
1172 DftiFreeDescriptor(&fft->ooplace[d]);
1173 }
1174 }
1175 if(fft->ooplace[3] != NULL)
1176 {
1177 DftiFreeDescriptor(&fft->ooplace[3]);
1178 }
1179 free(fft);
1180 }
1181 }
1182
1183 #else
1184 int
1185 gmx_fft_mkl_empty;
1186 #endif /* GMX_FFT_MKL */