Context Navigation

source: branches/HeuristicLab.Problems.GaussianProcessTuning/ILNumerics.2.14.4735.573/Native/FFT/ILMKLFFT.cs @ 10367

Visit:

Last change on this file since 10367 was 9102, checked in by gkronber, 12 years ago
#1967: ILNumerics source for experimentation
File size: 55.9 KB

Line
1	///
2	/// This file is part of ILNumerics Community Edition.
3	///
4	/// ILNumerics Community Edition - high performance computing for applications.
5	/// Copyright (C) 2006 - 2012 Haymo Kutschbach, http://ilnumerics.net
6	///
7	/// ILNumerics Community Edition is free software: you can redistribute it and/or modify
8	/// it under the terms of the GNU General Public License version 3 as published by
9	/// the Free Software Foundation.
10	///
11	/// ILNumerics Community Edition is distributed in the hope that it will be useful,
12	/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13	/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	/// GNU General Public License for more details.
15	///
16	/// You should have received a copy of the GNU General Public License
17	/// along with ILNumerics Community Edition. See the file License.txt in the root
18	/// of your distribution package. If not, see <http://www.gnu.org/licenses/>.
19	///
20	/// In addition this software uses the following components and/or licenses:
21	///
22	/// =================================================================================
23	/// The Open Toolkit Library License
24	///
25	/// Copyright (c) 2006 - 2009 the Open Toolkit library.
26	///
27	/// Permission is hereby granted, free of charge, to any person obtaining a copy
28	/// of this software and associated documentation files (the "Software"), to deal
29	/// in the Software without restriction, including without limitation the rights to
30	/// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
31	/// the Software, and to permit persons to whom the Software is furnished to do
32	/// so, subject to the following conditions:
33	///
34	/// The above copyright notice and this permission notice shall be included in all
35	/// copies or substantial portions of the Software.
36	///
37	/// =================================================================================
38	///
39
40	#pragma warning disable 1591
41
42	using System;
43	using System.Collections.Generic;
44	using System.Text;
45	using System.Security;
46	using System.Runtime.InteropServices;
47	using ILNumerics.Exceptions;
48	using ILNumerics.Misc;
49
50	namespace ILNumerics.Native {
51
52	/// <summary>
53	/// MKL configuration parameters (constant definitions)
54	/// </summary>
55	internal sealed class MKLParameter {
56	/* Domain for forward transform. No default value */
57	public static readonly int FORWARD_DOMAIN = 0;
58
59	/* Dimensionality; or rank. No default value */
60	public static readonly int DIMENSION = 1;
61
62	/* Length(s) of transform. No default value */
63	public static readonly int LENGTHS = 2;
64
65	/* Floating point precision. No default value */
66	public static readonly int PRECISION = 3;
67
68	/* Scale factor for forward transform [1.0] */
69	public static readonly int FORWARD_SCALE = 4;
70
71	/* Scale factor for backward transform [1.0] */
72	public static readonly int BACKWARD_SCALE = 5;
73
74	/* Exponent sign for forward transform [NEGATIVE] */
75	/* FORWARD_SIGN = 6; ## NOT IMPLEMENTED */
76
77	/* Number of data sets to be transformed [1] */
78	public static readonly int NUMBER_OF_TRANSFORMS = 7;
79
80	/* Storage of finite complex-valued sequences in complex domain
81	[COMPLEX_COMPLEX] */
82	public static readonly int COMPLEX_STORAGE = 8;
83
84	/* Storage of finite real-valued sequences in real domain
85	[REAL_REAL] */
86	public static readonly int REAL_STORAGE = 9;
87
88	/* Storage of finite complex-valued sequences in conjugate-even
89	domain [COMPLEX_REAL] */
90	public static readonly int CONJUGATE_EVEN_STORAGE = 10;
91
92	/* Placement of result [INPLACE] */
93	public static readonly int PLACEMENT = 11;
94
95	/* Generalized strides for input data layout [tight; row-major for C] */
96	public static readonly int INPUT_STRIDES = 12;
97
98	/* Generalized strides for output data layout [tight; row-major
99	for C] */
100	public static readonly int OUTPUT_STRIDES = 13;
101
102	/* Distance between first input elements for multiple transforms
103	[0] */
104	public static readonly int INPUT_DISTANCE = 14;
105
106	/* Distance between first output elements for multiple transforms
107	[0] */
108	public static readonly int OUTPUT_DISTANCE = 15;
109
110	/* Ordering of the result [ORDERED] */
111	public static readonly int ORDERING = 18;
112
113	/* Possible transposition of result [NONE] */
114	public static readonly int TRANSPOSE = 19;
115
116	/* User-settable descriptor name [""] */
117	public static readonly int DESCRIPTOR_NAME = 20; /* DEPRECATED */
118
119	/* Packing format for COMPLEX_REAL storage of finite
120	conjugate-even sequences [CCS_FORMAT] */
121	public static readonly int PACKED_FORMAT = 21;
122
123	/* Commit status of the descriptor - R/O parameter */
124	public static readonly int COMMIT_STATUS = 22;
125
126	/* Version string for this DFTI implementation - R/O parameter */
127	public static readonly int VERSION = 23;
128
129	/* Number of user threads that share the descriptor [1] */
130	public static readonly int NUMBER_OF_USER_THREADS = 26;
131	}
132	/// <summary>
133	/// MKL configuration values (constant definitions)
134	/// </summary>
135	internal sealed class MKLValues {
136	public static readonly int MKL_ALL = 0;
137	public static readonly int MKL_BLAS = 1;
138	public static readonly int MKL_FFT = 2;
139	public static readonly int MKL_VML = 3;
140	/* Values of the descriptor configuration parameters */
141	/* COMMIT_STATUS */
142	public static readonly int COMMITTED = 30;
143	public static readonly int UNCOMMITTED = 31;
144
145	/* FORWARD_DOMAIN */
146	public static readonly int COMPLEX = 32;
147	public static readonly int REAL = 33;
148	/* CONJUGATE_EVEN = 34; ## NOT IMPLEMENTED */
149
150	/* PRECISION */
151	public static readonly int SINGLE = 35;
152	public static readonly int DOUBLE = 36;
153
154	/* COMPLEX_STORAGE and CONJUGATE_EVEN_STORAGE */
155	public static readonly int COMPLEX_COMPLEX = 39;
156	public static readonly int COMPLEX_REAL = 40;
157
158	/* REAL_STORAGE */
159	public static readonly int REAL_COMPLEX = 41;
160	public static readonly int REAL_REAL = 42;
161
162	/* PLACEMENT */
163	public static readonly int INPLACE = 43; /* Result overwrites input */
164	public static readonly int NOT_INPLACE = 44; /* Have another place for result */
165
166	/* ORDERING */
167	public static readonly int ORDERED = 48;
168	public static readonly int BACKWARD_SCRAMBLED = 49;
169
170	/* Allow/avoid certain usages */
171	public static readonly int ALLOW = 51; /* Allow transposition or workspace */
172	public static readonly int NONE = 53;
173
174	/* PACKED_FORMAT (for storing congugate-even finite sequence
175	in real array) */
176	public static readonly int CCS_FORMAT = 54; /* Complex conjugate-symmetric */
177	public static readonly int PACK_FORMAT = 55; /* Pack format for real DFT */
178	public static readonly int PERM_FORMAT = 56; /* Perm format for real DFT */
179	public static readonly int CCE_FORMAT = 57; /* Complex conjugate-even */
180
181	/* Error classes */
182	public static readonly int NO_ERROR = 0;
183	public static readonly int MEMORY_ERROR = 1;
184	public static readonly int INVALID_CONFIGURATION = 2;
185	public static readonly int INCONSISTENT_CONFIGURATION = 3;
186	public static readonly int MULTITHREADED_ERROR = 4;
187	public static readonly int BAD_DESCRIPTOR = 5;
188	public static readonly int UNIMPLEMENTED = 6;
189	public static readonly int MKL_INTERNAL_ERROR = 7;
190	public static readonly int NUMBER_OF_THREADS_ERROR = 8;
191	public static readonly int ONED_LENGTH_EXCEEDS_INT32 = 9;
192
193	public static readonly int MAX_MESSAGE_LENGTH = 40; /* Maximum length of error string */
194	public static readonly int MAX_NAME_LENGTH = 10; /* Maximum length of descriptor name */
195	public static readonly int VERSION_LENGTH = 198; /* Maximum length of MKL version string */
196	}
197	/// <summary>
198	/// import functions (pinvoke)
199	/// </summary>
200	internal sealed class MKLImports {
201
202	#region pinvoke definitions
203
204	#if x64
205	#region 64 bit includes
206	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
207	internal static extern int mkl_domain_set_num_threads(ref int num, ref int mask);
208	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
209	internal static extern int mkl_domain_get_max_threads(ref int mask);
210	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
211	internal static extern void mkl_set_num_threads(int num);
212	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
213	internal static extern int mkl_get_max_threads();
214	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl)]
215	internal static extern void omp_set_num_threads(int num);
216	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
217	internal static extern int omp_get_num_threads();
218
219	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
220	internal static extern void mkl_set_dynamic(ref int boolean_value);
221
222	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
223	public static extern int DftiCreateDescriptor(ref IntPtr pDescriptor, int precision, int domain, int dimCount, int dim1);
224	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
225	public static extern int DftiCreateDescriptor(ref IntPtr pDescriptor, int precision, int domain, int dimCount, int[] dims);
226	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
227	public static extern int DftiFreeDescriptor(ref IntPtr pDescriptor);
228	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
229	public static extern int DftiCommitDescriptor(IntPtr pDescriptor);
230	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
231	public static extern int DftiErrorClass(int status, int ERROR_CLASS);
232	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
233	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, int value);
234	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
235	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, double value);
236	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
237	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, int[] values);
238	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
239	public static extern int DftiGetValue(IntPtr pDescriptor, int parameter, ref int value);
240	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
241	public static extern int DftiGetValue(IntPtr pDescriptor, int parameter, ref double value);
242	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
243	public static extern int DftiComputeForward(IntPtr pDescriptor, IntPtr pArray);
244	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
245	public static extern int DftiComputeForward(IntPtr pDescriptor, IntPtr pInArray, IntPtr pOutArray);
246
247	/** DFTI native DftiComputeForward declaration */
248	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
249	internal static extern int DftiComputeForward(IntPtr desc, [In] double[] x_in, [Out] double[] x_out);
250
251	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
252	public static extern int DftiComputeBackward(IntPtr pDescriptor, IntPtr pArray);
253	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
254	public static extern int DftiComputeBackward(IntPtr pDescriptor, IntPtr pInArray, IntPtr pOutArray);
255
256	/** DFTI native DftiComputeBackward declaration */
257	[DllImport("mkl_custom64", CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
258	internal static extern int DftiComputeBackward(IntPtr desc, [In] double[] x_in, [Out] double[] x_out);
259
260	[DllImport("mkl_custom64",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
261	public static extern string DftiErrorMessage(int errID);
262	#endregion
263	#else
264	#region 32 bit includes
265	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
266	internal static extern int mkl_domain_set_num_threads(ref int num, ref int mask);
267	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
268	internal static extern int mkl_domain_get_max_threads(ref int mask);
269
270	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
271	internal static extern void mkl_set_num_threads(int num);
272	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
273	internal static extern int mkl_get_max_threads();
274	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
275	internal static extern void omp_set_num_threads(int num);
276	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
277	internal static extern void mkl_set_dynamic(ref int boolean_value);
278	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
279	public static extern int DftiCreateDescriptor(ref IntPtr pDescriptor, int precision, int domain, int dimCount, int dim1);
280	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
281	public static extern int DftiCreateDescriptor(ref IntPtr pDescriptor, int precision, int domain, int dimCount, int[] dims);
282	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
283	public static extern int DftiFreeDescriptor(ref IntPtr pDescriptor);
284	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
285	public static extern int DftiCommitDescriptor(IntPtr pDescriptor);
286	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
287	public static extern int DftiErrorClass(int status, int ERROR_CLASS);
288
289	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
290	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, int value);
291	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
292	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, double value);
293	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
294	public static extern int DftiSetValue(IntPtr pDescriptor, int parameter, int[] values);
295
296	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
297	public static extern int DftiGetValue(IntPtr pDescriptor, int parameter, ref int value);
298	[DllImport("mkl_custom32",CallingConvention =CallingConvention.Cdecl),SuppressUnmanagedCodeSecurity]
299	public static extern int DftiGetValue(IntPtr pDescriptor, int parameter, ref double value);
300
301	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
302	public static extern int DftiComputeForward(IntPtr pDescriptor, IntPtr pArray);
303	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
304	public static extern int DftiComputeForward(IntPtr pDescriptor, IntPtr pInArray, IntPtr pOutArray);
305	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
306	public static extern int DftiComputeBackward(IntPtr pDescriptor, IntPtr pArray);
307	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
308	public static extern int DftiComputeBackward(IntPtr pDescriptor, IntPtr pInArray, IntPtr pOutArray);
309	[DllImport("mkl_custom32", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
310	public static extern string DftiErrorMessage(int errID);
311	#endregion
312	#endif
313	#endregion
314	}
315
316	internal sealed class MKLDescriptor : IDisposable {
317	private int[] Stride;
318	private int Precision;
319	private int RealComplex;
320	private bool InPlace;
321	private ILSize Dimension;
322	private int AlongDimensionIndex;
323	private IntPtr Descriptor = IntPtr.Zero;
324	private int Dimensionality;
325	private int LocalLoopIncrement;
326	private int LocalLoopLength;
327
328	public IntPtr GetDescriptor(int precision, int realComplex, int dimensionality,
329	int dim, ILSize dimension, bool inplace,
330	out int localLoopLength, out int localLoopIncrement) {
331	// DEBUG: remove "false" for descriptor caching!
332	if (Descriptor != IntPtr.Zero
333	&& precision == Precision
334	&& realComplex == RealComplex
335	&& dimensionality == Dimensionality
336	&& dimension.IsSameShape(Dimension)
337	&& dim == AlongDimensionIndex
338	&& inplace == InPlace) {
339	// reuse old descriptor
340	localLoopIncrement = LocalLoopIncrement;
341	localLoopLength = LocalLoopLength;
342	return Descriptor;
343	} else {
344
345	// must create a new one
346	if (Descriptor != IntPtr.Zero) {
347	MKLImports.DftiFreeDescriptor(ref Descriptor);
348	}
349	IntPtr tmpDescriptor = IntPtr.Zero;
350	int error;
351	if (dimensionality == 1) {
352	#region one dimensional transforms
353	error = MKLImports.DftiCreateDescriptor(ref tmpDescriptor, precision ,realComplex, dimensionality, dimension[dim]);
354	if (error != 0) {
355	throw new ILInvalidOperationException ("error creating descriptor: " + MKLImports.DftiErrorMessage(error));
356	}
357
358	// storage,number of subsequent transformations
359	int nrTransGlobal = dimension.NumberOfElements / dimension[dim];
360	int nrTransPerPass;
361	if (dim == 0) {
362	// MKL does all transf. at once
363	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.NUMBER_OF_TRANSFORMS, nrTransGlobal);
364	int distances = dimension[0];
365	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.INPUT_DISTANCE, distances);
366	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.OUTPUT_DISTANCE, distances);
367	localLoopIncrement = 1;
368	localLoopLength = 1;
369	nrTransPerPass = 1;
370	} else {
371	nrTransPerPass = dimension.SequentialIndexDistance(dim);
372	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.NUMBER_OF_TRANSFORMS, nrTransPerPass);
373	int distances = 1;
374	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.INPUT_DISTANCE, distances);
375	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.OUTPUT_DISTANCE, distances);
376	localLoopIncrement = dimension.SequentialIndexDistance(dim) * dimension[dim];
377	localLoopLength = nrTransGlobal / nrTransPerPass;
378	}
379
380	// spacing between elements
381	if (Stride == null \|\| Stride.Length < 2) {
382	Stride = new int[10];
383	}
384	Stride[1] = nrTransPerPass;
385	MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.INPUT_STRIDES, Stride);
386	error = MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.OUTPUT_STRIDES, Stride);
387	if (error != 0) {
388	throw new ILInvalidOperationException("error: " + MKLImports.DftiErrorMessage(error));
389	}
390	#endregion
391	} else {
392	#region multidimensional transforms
393	int[] tmpDims = dimension.ToIntArray(10);
394	error = MKLImports.DftiCreateDescriptor(ref tmpDescriptor, precision ,realComplex, dimensionality, tmpDims);
395	// MKL will handle to do all transforms at once
396	localLoopIncrement = 1;
397	localLoopLength = 1;
398	// how many transformations ?
399	int nrElementsPerPass = dimension.SequentialIndexDistance(dimensionality);
400	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.NUMBER_OF_TRANSFORMS, dimension.NumberOfElements / nrElementsPerPass);
401	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.INPUT_DISTANCE, nrElementsPerPass);
402	MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.OUTPUT_DISTANCE, nrElementsPerPass);
403
404	// spacing between elements
405	if (Stride == null \|\| Stride.Length < dimensionality + 1) {
406	Stride = new int[dimensionality + 1];
407	}
408	int[] seqDistances = dimension.GetSequentialIndexDistances(dimensionality);
409	for (int i = dimensionality; i-- > 0; ) {
410	Stride[i + 1] = seqDistances[i];
411	}
412	//Stride[0] = 0; // should never gotten changed
413	MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.INPUT_STRIDES, Stride);
414	error = MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.OUTPUT_STRIDES, Stride);
415	if (error != 0) {
416	throw new ILInvalidOperationException("error: " + MKLImports.DftiErrorMessage(error));
417	}
418	#endregion
419	}
420	//// backward scale
421	//if (!forwardTransform) {
422	// MKLImports.DftiSetValue(tmpDescriptor,MKLParameter.BACKWARD_SCALE, __arglist(1.0/Dimension[dim]));
423	//}
424	if (!inplace) {
425	error = MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.PLACEMENT, MKLValues.NOT_INPLACE);
426	//error = MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.REAL_STORAGE, __arglist(MKLValues.REAL_REAL));
427	error = MKLImports.DftiSetValue(tmpDescriptor, MKLParameter.CONJUGATE_EVEN_STORAGE, MKLValues.COMPLEX_COMPLEX);
428	}
429
430	error = MKLImports.DftiCommitDescriptor(tmpDescriptor);
431	if (error != 0) {
432	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
433	}
434
435	Precision = precision;
436	RealComplex = realComplex;
437	Dimensionality = dimensionality;
438	Dimension = dimension;
439	AlongDimensionIndex = dim;
440	Descriptor = tmpDescriptor;
441	LocalLoopIncrement = localLoopIncrement;
442	LocalLoopLength = localLoopLength;
443	InPlace = inplace;
444	return Descriptor;
445	}
446	}
447	public string GetDescriptorInfo() {
448	if (Descriptor != IntPtr.Zero) {
449	return GetDescriptorInfo(Descriptor);
450	} else {
451	return "(no descriptor set)";
452	}
453	}
454	public string GetDescriptorInfo(IntPtr descriptor) {
455	StringBuilder ret = new StringBuilder();
456	int val = 0;
457	MKLImports.DftiGetValue(descriptor, MKLParameter.PRECISION, ref val); ret.Append("Precision: " + val + " " + Environment.NewLine);
458	MKLImports.DftiGetValue(descriptor, MKLParameter.FORWARD_DOMAIN, ref val); ret.Append("ForwDomain: " + val + " " + Environment.NewLine);
459	MKLImports.DftiGetValue(descriptor, MKLParameter.DIMENSION, ref val); ret.Append("DIMENSION: " + val + " " + Environment.NewLine);
460	MKLImports.DftiGetValue(descriptor, MKLParameter.LENGTHS, ref val); ret.Append("Length: " + val + " " + Environment.NewLine);
461	MKLImports.DftiGetValue(descriptor, MKLParameter.PLACEMENT, ref val); ret.Append("PLACEMENT: " + val + " " + Environment.NewLine);
462	MKLImports.DftiGetValue(descriptor, MKLParameter.NUMBER_OF_TRANSFORMS, ref val); ret.Append("NUMBER_OF_TRANSFORMS: " + val + " " + Environment.NewLine);
463	MKLImports.DftiGetValue(descriptor, MKLParameter.COMPLEX_STORAGE, ref val); ret.Append("COMPLEX_STORAGE: " + val + " " + Environment.NewLine);
464	MKLImports.DftiGetValue(descriptor, MKLParameter.CONJUGATE_EVEN_STORAGE, ref val); ret.Append("CONJUGATE_EVEN_STORAGE: " + val + " " + Environment.NewLine);
465	MKLImports.DftiGetValue(descriptor, MKLParameter.INPUT_DISTANCE, ref val); ret.Append("INPUT_DISTANCE: " + val + " " + Environment.NewLine);
466	MKLImports.DftiGetValue(descriptor, MKLParameter.OUTPUT_DISTANCE, ref val); ret.Append("OUTPUT_DISTANCE: " + val + " " + Environment.NewLine);
467	MKLImports.DftiGetValue(descriptor, MKLParameter.INPUT_STRIDES, ref val); ret.Append("INPUT_STRIDES: " + val + " " + Environment.NewLine);
468	MKLImports.DftiGetValue(descriptor, MKLParameter.OUTPUT_STRIDES, ref val); ret.Append("OUTPUT_STRIDES: " + val + " " + Environment.NewLine);
469	return ret.ToString();
470	}
471
472	#region IDisposable Members
473
474	public void Dispose() {
475	lock (this)
476	if (Descriptor != IntPtr.Zero) {
477	Descriptor = IntPtr.Zero;
478	MKLImports.DftiFreeDescriptor(ref Descriptor);
479	}
480	}
481	~MKLDescriptor () {
482	Dispose();
483	}
484	#endregion
485	}
486	/// <summary>
487	/// Wrapper for FFT interface using MKL 10_03
488	/// </summary>
489	public class ILMKLFFT : IILFFT {
490
491	public static void Init() {
492	try {
493	// unless the user explicitely configured the number of threads, let omp configure it!
494	if (Settings.MaxNumberThreadsConfigured) {
495	SetNumThreads(Settings.MaxNumberThreads);
496	}
497	} catch (System.IO.FileNotFoundException) {}
498	}
499	public static void SetNumThreads(int numThreads) {
500	int mklfft = MKLValues.MKL_FFT;
501	MKLImports.mkl_domain_set_num_threads(ref numThreads,ref mklfft);
502	}
503	public static int GetMaxThreads() {
504	int mklfft = MKLValues.MKL_FFT;
505	return MKLImports.mkl_domain_get_max_threads(ref mklfft);
506	}
507	private static void SetDynamicThreads(bool dynamic) {
508	int val = dynamic ? 1 : 0;
509	MKLImports.mkl_set_dynamic(ref val);
510	}
511
512	#region attributes
513	[ThreadStatic]
514	private static int[] s_strides;
515	private static int[] Strides {
516	get {
517	if (s_strides == null) {
518	s_strides = new int[10];
519	}
520	return s_strides;
521	}
522	}
523	[ThreadStatic]
524	private static MKLDescriptor s_descriptor;
525	private static MKLDescriptor DescriptorCache {
526	get {
527	if (s_descriptor == null) {
528	s_descriptor = new MKLDescriptor();
529	}
530	return s_descriptor;
531	}
532	}
533	#endregion
534
535	#region constructor
536	public ILMKLFFT () {
537	Init();
538	}
539	#endregion
540
541	#region IILFFT Member - 1-D
542
543
544
545	public ILRetArray< complex > FFTForward1D (ILInArray< double > A, int dim) {
546	using (ILScope.Enter(A)) {
547	if (object.Equals(A,null) \|\| dim < 0)
548	throw new ILArgumentException("invalid parameter");
549	if (A.IsEmpty) return ILMath.empty< complex >(A.Size);
550	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
551	return ILMath.tocomplex(A);
552	// prepare output array
553	ILArray< complex > ret = ILMath.tocomplex(A) ;
554	int localLoopLength, localLoopIncrement;
555
556	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE,MKLValues.COMPLEX,1,dim,A.Size,true,
557	out localLoopLength, out localLoopIncrement);
558	int error = 0;
559	// do the transform(s)
560	unsafe {
561	fixed ( complex * retArr = ret.GetArrayForWrite()) {
562	for (int i = 0; i < localLoopLength && error == 0; i ++)
563	error = MKLImports.DftiComputeForward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
564	}
565	}
566	if (isMKLError(error)) {
567	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
568	}
569
570
571	return ret;
572	}
573	}
574
575
576	#region HYCALPER AUTO GENERATED CODE
577
578
579	public ILRetArray< fcomplex > FFTBackward1D (ILInArray< fcomplex > A, int dim) {
580	using (ILScope.Enter(A)) {
581	if (object.Equals(A,null) \|\| dim < 0)
582	throw new ILArgumentException("invalid parameter");
583	if (A.IsEmpty) return ILMath.empty< fcomplex >(A.Size);
584	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
585	return A.C;
586	// prepare output array
587	ILArray< fcomplex > ret = A.C ;
588	int localLoopLength, localLoopIncrement;
589
590	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE,MKLValues.COMPLEX,1,dim,A.Size,true,
591	out localLoopLength, out localLoopIncrement);
592	int error = 0;
593	// do the transform(s)
594	unsafe {
595	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
596	for (int i = 0; i < localLoopLength && error == 0; i ++)
597	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
598	}
599	}
600	if (isMKLError(error)) {
601	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
602	}
603	ret = ret / new fcomplex(A.D[dim],0);
604	return ret;
605	}
606	}
607
608
609	public ILRetArray< complex > FFTBackward1D (ILInArray< complex > A, int dim) {
610	using (ILScope.Enter(A)) {
611	if (object.Equals(A,null) \|\| dim < 0)
612	throw new ILArgumentException("invalid parameter");
613	if (A.IsEmpty) return ILMath.empty< complex >(A.Size);
614	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
615	return A.C;
616	// prepare output array
617	ILArray< complex > ret = A.C ;
618	int localLoopLength, localLoopIncrement;
619
620	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE,MKLValues.COMPLEX,1,dim,A.Size,true,
621	out localLoopLength, out localLoopIncrement);
622	int error = 0;
623	// do the transform(s)
624	unsafe {
625	fixed ( complex * retArr = ret.GetArrayForWrite()) {
626	for (int i = 0; i < localLoopLength && error == 0; i ++)
627	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
628	}
629	}
630	if (isMKLError(error)) {
631	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
632	}
633	ret = ret / new complex(A.D[dim],0);
634	return ret;
635	}
636	}
637
638
639	public ILRetArray< fcomplex > FFTForward1D (ILInArray< fcomplex > A, int dim) {
640	using (ILScope.Enter(A)) {
641	if (object.Equals(A,null) \|\| dim < 0)
642	throw new ILArgumentException("invalid parameter");
643	if (A.IsEmpty) return ILMath.empty< fcomplex >(A.Size);
644	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
645	return A.C;
646	// prepare output array
647	ILArray< fcomplex > ret = A.C ;
648	int localLoopLength, localLoopIncrement;
649
650	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE,MKLValues.COMPLEX,1,dim,A.Size,true,
651	out localLoopLength, out localLoopIncrement);
652	int error = 0;
653	// do the transform(s)
654	unsafe {
655	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
656	for (int i = 0; i < localLoopLength && error == 0; i ++)
657	error = MKLImports.DftiComputeForward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
658	}
659	}
660	if (isMKLError(error)) {
661	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
662	}
663
664	return ret;
665	}
666	}
667
668
669	public ILRetArray< fcomplex > FFTForward1D (ILInArray< float > A, int dim) {
670	using (ILScope.Enter(A)) {
671	if (object.Equals(A,null) \|\| dim < 0)
672	throw new ILArgumentException("invalid parameter");
673	if (A.IsEmpty) return ILMath.empty< fcomplex >(A.Size);
674	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
675	return ILMath.tofcomplex(A);
676	// prepare output array
677	ILArray< fcomplex > ret = ILMath.tofcomplex(A) ;
678	int localLoopLength, localLoopIncrement;
679
680	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE,MKLValues.COMPLEX,1,dim,A.Size,true,
681	out localLoopLength, out localLoopIncrement);
682	int error = 0;
683	// do the transform(s)
684	unsafe {
685	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
686	for (int i = 0; i < localLoopLength && error == 0; i ++)
687	error = MKLImports.DftiComputeForward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
688	}
689	}
690	if (isMKLError(error)) {
691	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
692	}
693
694	return ret;
695	}
696	}
697
698
699	public ILRetArray< complex > FFTForward1D (ILInArray< complex > A, int dim) {
700	using (ILScope.Enter(A)) {
701	if (object.Equals(A,null) \|\| dim < 0)
702	throw new ILArgumentException("invalid parameter");
703	if (A.IsEmpty) return ILMath.empty< complex >(A.Size);
704	if (dim >= A.Size.NumberOfDimensions \|\| A.Size[dim] == 1)
705	return A.C;
706	// prepare output array
707	ILArray< complex > ret = A.C ;
708	int localLoopLength, localLoopIncrement;
709
710	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE,MKLValues.COMPLEX,1,dim,A.Size,true,
711	out localLoopLength, out localLoopIncrement);
712	int error = 0;
713	// do the transform(s)
714	unsafe {
715	fixed ( complex * retArr = ret.GetArrayForWrite()) {
716	for (int i = 0; i < localLoopLength && error == 0; i ++)
717	error = MKLImports.DftiComputeForward(descriptor, (IntPtr)(retArr + i * localLoopIncrement));
718	}
719	}
720	if (isMKLError(error)) {
721	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
722	}
723
724	return ret;
725	}
726	}
727
728
729	#endregion HYCALPER AUTO GENERATED CODE
730
731
732
733
734	public ILRetArray< double > FFTBackwSym1D(ILInArray< complex > A, int dim) {
735	using (ILScope.Enter(A)) {
736	if (object.Equals(A,null) \|\| dim < 0)
737	throw new ILArgumentException("invalid parameter");
738	if (A.IsEmpty) return ILArray< double >.empty(A.Size);
739	if (A.IsScalar \|\| A.Size[dim] == 1)
740	return ILMath.real(A.C);
741	// prepare output array, if we query the memory directly from
742	// memory pool, we prevent from clearing the elements since every
743	// single one will be overwritten by fft anyway
744	int error = 0,inDim = A.Size[dim];
745	double [] retArr = ILMemoryPool.Pool.New< double >(A.Size.NumberOfElements);
746	ILArray< double > ret = ILMath.array< double >(retArr,A.Size);
747	int localLoopLength, localLoopIncrement;
748	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE , MKLValues.REAL, 1, dim, A.S, false,
749	out localLoopLength, out localLoopIncrement);
750	// do the transform(s)
751	unsafe {
752	fixed ( double * retArrP = ret.GetArrayForWrite())
753	fixed ( complex * pA = A.GetArrayForRead()) {
754	for (int i = 0; i < (localLoopLength * localLoopIncrement) && error == 0; i += localLoopIncrement)
755	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(pA + i), (IntPtr)(retArrP + i));
756	}
757	}
758	if (isMKLError(error)) {
759	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
760	}
761	return ret / inDim;
762	}
763	}
764
765
766	#region HYCALPER AUTO GENERATED CODE
767
768
769
770	public ILRetArray< float > FFTBackwSym1D(ILInArray< fcomplex > A, int dim) {
771	using (ILScope.Enter(A)) {
772	if (object.Equals(A,null) \|\| dim < 0)
773	throw new ILArgumentException("invalid parameter");
774	if (A.IsEmpty) return ILArray< float >.empty(A.Size);
775	if (A.IsScalar \|\| A.Size[dim] == 1)
776	return ILMath.real(A.C);
777	// prepare output array, if we query the memory directly from
778	// memory pool, we prevent from clearing the elements since every
779	// single one will be overwritten by fft anyway
780	int error = 0,inDim = A.Size[dim];
781	float [] retArr = ILMemoryPool.Pool.New< float >(A.Size.NumberOfElements);
782	ILArray< float > ret = ILMath.array< float >(retArr,A.Size);
783	int localLoopLength, localLoopIncrement;
784	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE , MKLValues.REAL, 1, dim, A.S, false,
785	out localLoopLength, out localLoopIncrement);
786	// do the transform(s)
787	unsafe {
788	fixed ( float * retArrP = ret.GetArrayForWrite())
789	fixed ( fcomplex * pA = A.GetArrayForRead()) {
790	for (int i = 0; i < (localLoopLength * localLoopIncrement) && error == 0; i += localLoopIncrement)
791	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(pA + i), (IntPtr)(retArrP + i));
792	}
793	}
794	if (isMKLError(error)) {
795	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
796	}
797	return ret / inDim;
798	}
799	}
800
801
802	#endregion HYCALPER AUTO GENERATED CODE
803
804	#endregion
805
806	#region IILFFT Member - n-D
807
808
809
810	public ILRetArray< complex > FFTForward (ILInArray< double > A, int nDims) {
811	using (ILScope.Enter(A)) {
812	if (object.Equals(A, null) \|\| nDims <= 0 )
813	throw new ILArgumentException("invalid argument!");
814	if (A.IsEmpty) return ILArray< complex > .empty(A.Size);
815	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
816	return ILMath.tocomplex(A) ;
817	if (nDims > A.Size.NumberOfDimensions)
818	nDims = A.Size.NumberOfDimensions;
819
820	if (nDims > 3) return FFTForward(ILMath.tocomplex(A),nDims);
821	// prepare output array
822	ILArray< complex > ret = ILMath.tocomplex(A) ;
823	int localLoopLength, localLoopIncrement;
824	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
825	out localLoopLength, out localLoopIncrement);
826	int error = 0;
827	// do the transform(s)
828	unsafe {
829	fixed ( complex * retArr = ret.GetArrayForWrite()) {
830	error = MKLImports.DftiComputeForward (descriptor, (IntPtr)retArr);
831	}
832	}
833	if (isMKLError(error)) {
834	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
835	}
836
837
838	return ret;
839	}
840	}
841
842	#region HYCALPER AUTO GENERATED CODE
843
844
845	public ILRetArray< complex > FFTBackward (ILInArray< complex > A, int nDims) {
846	using (ILScope.Enter(A)) {
847	if (object.Equals(A, null) \|\| nDims <= 0 )
848	throw new ILArgumentException("invalid argument!");
849	if (A.IsEmpty) return ILArray< complex > .empty(A.Size);
850	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
851	return A.C ;
852	if (nDims > A.Size.NumberOfDimensions)
853	nDims = A.Size.NumberOfDimensions;
854
855	// prepare output array
856	ILArray< complex > ret = A.C ;
857	int localLoopLength, localLoopIncrement;
858	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
859	out localLoopLength, out localLoopIncrement);
860	int error = 0;
861	// do the transform(s)
862	unsafe {
863	fixed ( complex * retArr = ret.GetArrayForWrite()) {
864	error = MKLImports.DftiComputeBackward (descriptor, (IntPtr)retArr);
865	}
866	}
867	if (isMKLError(error)) {
868	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
869	}
870	ret = ret / new complex(A.D.SequentialIndexDistance(nDims),0);
871	return ret;
872	}
873	}
874
875	public ILRetArray< fcomplex > FFTBackward (ILInArray< fcomplex > A, int nDims) {
876	using (ILScope.Enter(A)) {
877	if (object.Equals(A, null) \|\| nDims <= 0 )
878	throw new ILArgumentException("invalid argument!");
879	if (A.IsEmpty) return ILArray< fcomplex > .empty(A.Size);
880	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
881	return A.C ;
882	if (nDims > A.Size.NumberOfDimensions)
883	nDims = A.Size.NumberOfDimensions;
884
885	// prepare output array
886	ILArray< fcomplex > ret = A.C ;
887	int localLoopLength, localLoopIncrement;
888	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
889	out localLoopLength, out localLoopIncrement);
890	int error = 0;
891	// do the transform(s)
892	unsafe {
893	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
894	error = MKLImports.DftiComputeBackward (descriptor, (IntPtr)retArr);
895	}
896	}
897	if (isMKLError(error)) {
898	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
899	}
900	ret = ret / new fcomplex(A.D.SequentialIndexDistance(nDims),0);
901	return ret;
902	}
903	}
904
905	public ILRetArray< fcomplex > FFTForward (ILInArray< fcomplex > A, int nDims) {
906	using (ILScope.Enter(A)) {
907	if (object.Equals(A, null) \|\| nDims <= 0 )
908	throw new ILArgumentException("invalid argument!");
909	if (A.IsEmpty) return ILArray< fcomplex > .empty(A.Size);
910	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
911	return A.C ;
912	if (nDims > A.Size.NumberOfDimensions)
913	nDims = A.Size.NumberOfDimensions;
914
915	// prepare output array
916	ILArray< fcomplex > ret = A.C ;
917	int localLoopLength, localLoopIncrement;
918	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
919	out localLoopLength, out localLoopIncrement);
920	int error = 0;
921	// do the transform(s)
922	unsafe {
923	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
924	error = MKLImports.DftiComputeForward (descriptor, (IntPtr)retArr);
925	}
926	}
927	if (isMKLError(error)) {
928	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
929	}
930
931	return ret;
932	}
933	}
934
935	public ILRetArray< fcomplex > FFTForward (ILInArray< float > A, int nDims) {
936	using (ILScope.Enter(A)) {
937	if (object.Equals(A, null) \|\| nDims <= 0 )
938	throw new ILArgumentException("invalid argument!");
939	if (A.IsEmpty) return ILArray< fcomplex > .empty(A.Size);
940	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
941	return ILMath.tofcomplex(A) ;
942	if (nDims > A.Size.NumberOfDimensions)
943	nDims = A.Size.NumberOfDimensions;
944	if (nDims > 3) return FFTForward(ILMath.tofcomplex(A),nDims);
945	// prepare output array
946	ILArray< fcomplex > ret = ILMath.tofcomplex(A) ;
947	int localLoopLength, localLoopIncrement;
948	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
949	out localLoopLength, out localLoopIncrement);
950	int error = 0;
951	// do the transform(s)
952	unsafe {
953	fixed ( fcomplex * retArr = ret.GetArrayForWrite()) {
954	error = MKLImports.DftiComputeForward (descriptor, (IntPtr)retArr);
955	}
956	}
957	if (isMKLError(error)) {
958	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
959	}
960
961	return ret;
962	}
963	}
964
965	public ILRetArray< complex > FFTForward (ILInArray< complex > A, int nDims) {
966	using (ILScope.Enter(A)) {
967	if (object.Equals(A, null) \|\| nDims <= 0 )
968	throw new ILArgumentException("invalid argument!");
969	if (A.IsEmpty) return ILArray< complex > .empty(A.Size);
970	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
971	return A.C ;
972	if (nDims > A.Size.NumberOfDimensions)
973	nDims = A.Size.NumberOfDimensions;
974
975	// prepare output array
976	ILArray< complex > ret = A.C ;
977	int localLoopLength, localLoopIncrement;
978	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE , MKLValues.COMPLEX, nDims, 0, A.Size, true,
979	out localLoopLength, out localLoopIncrement);
980	int error = 0;
981	// do the transform(s)
982	unsafe {
983	fixed ( complex * retArr = ret.GetArrayForWrite()) {
984	error = MKLImports.DftiComputeForward (descriptor, (IntPtr)retArr);
985	}
986	}
987	if (isMKLError(error)) {
988	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
989	}
990
991	return ret;
992	}
993	}
994
995	#endregion HYCALPER AUTO GENERATED CODE
996
997
998
999	public ILRetArray< double > FFTBackwSym(ILInArray< complex > A, int nDims) {
1000	using (ILScope.Enter(A)) {
1001	if (object.Equals(A, null) \|\| nDims <= 0 )
1002	throw new ILArgumentException("invalid argument!");
1003	if (A.IsEmpty) return ILMath.empty< double > (A.Size);
1004	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
1005	return ILMath.real(A.C);
1006	if (nDims > A.Size.NumberOfDimensions)
1007	nDims = A.Size.NumberOfDimensions;
1008	// MKL currently does not handle more than 3 dimensions this way
1009	if (nDims > 3) {
1010	return ILMath.real(FFTBackward(A,nDims));
1011	}
1012	// prepare output array, if we query the memory directly from
1013	// memory pool, we prevent from clearing the elements since every
1014	// single one will be overwritten by fft anyway
1015	double [] retArr = ILMemoryPool.Pool.New< double >(A.Size.NumberOfElements);
1016	ILArray< double> ret = ILMath.array< double>(retArr, A.Size);
1017	int localLoopIncrement, localLoopLength;
1018	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.DOUBLE , MKLValues.REAL,nDims, 0, A.S, false,
1019	out localLoopIncrement, out localLoopLength);
1020	int error = 0;
1021	// do the transform(s)
1022	unsafe {
1023	fixed ( double * retArrP = retArr)
1024	fixed ( complex * inArr = A.GetArrayForRead()) {
1025	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(inArr),(IntPtr)(retArrP));
1026	}
1027	}
1028	if (isMKLError(error)) {
1029	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
1030	}
1031	ret = ret / A.S.SequentialIndexDistance(nDims);
1032	return ret;
1033	}
1034	}
1035
1036	#region HYCALPER AUTO GENERATED CODE
1037
1038
1039	public ILRetArray< float > FFTBackwSym(ILInArray< fcomplex > A, int nDims) {
1040	using (ILScope.Enter(A)) {
1041	if (object.Equals(A, null) \|\| nDims <= 0 )
1042	throw new ILArgumentException("invalid argument!");
1043	if (A.IsEmpty) return ILMath.empty< float > (A.Size);
1044	if (A.IsScalar \|\| (A.Size[0] == 1 && nDims == 1))
1045	return ILMath.real(A.C);
1046	if (nDims > A.Size.NumberOfDimensions)
1047	nDims = A.Size.NumberOfDimensions;
1048	// MKL currently does not handle more than 3 dimensions this way
1049	if (nDims > 3) {
1050	return ILMath.real(FFTBackward(A,nDims));
1051	}
1052	// prepare output array, if we query the memory directly from
1053	// memory pool, we prevent from clearing the elements since every
1054	// single one will be overwritten by fft anyway
1055	float [] retArr = ILMemoryPool.Pool.New< float >(A.Size.NumberOfElements);
1056	ILArray< float> ret = ILMath.array< float>(retArr, A.Size);
1057	int localLoopIncrement, localLoopLength;
1058	IntPtr descriptor = DescriptorCache.GetDescriptor( MKLValues.SINGLE , MKLValues.REAL,nDims, 0, A.S, false,
1059	out localLoopIncrement, out localLoopLength);
1060	int error = 0;
1061	// do the transform(s)
1062	unsafe {
1063	fixed ( float * retArrP = retArr)
1064	fixed ( fcomplex * inArr = A.GetArrayForRead()) {
1065	error = MKLImports.DftiComputeBackward(descriptor, (IntPtr)(inArr),(IntPtr)(retArrP));
1066	}
1067	}
1068	if (isMKLError(error)) {
1069	throw new ILInvalidOperationException ("error: " + MKLImports.DftiErrorMessage(error));
1070	}
1071	ret = ret / A.S.SequentialIndexDistance(nDims);
1072	return ret;
1073	}
1074	}
1075
1076	#endregion HYCALPER AUTO GENERATED CODE
1077
1078	#endregion
1079
1080	#region IILFFT Member - Misc
1081
1082	public bool CachePlans {
1083	get {
1084	return true;
1085	}
1086	}
1087
1088	public void FreePlans() {
1089	DescriptorCache.Dispose();
1090	}
1091
1092	public bool SpeedyHermitian {
1093	get { return true; }
1094	}
1095
1096	#endregion
1097
1098	#region private helper
1099	private static bool isMKLError(int error) {
1100	// TODO check! only the first 32 bits seem to be relevant for error!
1101	//return MKLImports.DftiErrorClass(error, MKLValues.NO_ERROR) != 1;
1102	return error != MKLValues.NO_ERROR;
1103	}
1104	#endregion
1105
1106	}
1107	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences