Context Navigation

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

Visit:

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

Rev	Line
[9102]	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