Context Navigation

source: branches/HeuristicLab.Problems.GaussianProcessTuning/ILNumerics.2.14.4735.573/Functions/builtin/apply.cs @ 11194

Visit:

Last change on this file since 11194 was 9102, checked in by gkronber, 12 years ago
#1967: ILNumerics source for experimentation
File size: 130.8 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	using System;
	41	using System.Collections.Generic;
	42	using System.Text;
	43	using ILNumerics;
	44	using ILNumerics.Exceptions;
	45	using ILNumerics.Storage;
	46	using ILNumerics.Misc;
	47
	48
	49	namespace ILNumerics {
	50	public partial class ILMath {
	51
	52
	53	/// <summary>Apply an arbitrary function to two arrays</summary>
	54	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	55	/// <param name="A">Input array A</param>
	56	/// <param name="B">Input array B</param>
	57	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	58	/// <item>
	59	/// <term>size(A) == size(B)</term>
	60	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	61	/// </item>
	62	/// <item>
	63	/// <term>isscalar(A) \|\| isscalar(B)</term>
	64	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	65	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	66	/// </item>
	67	/// <item>
	68	/// <term>All other cases</term>
	69	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	70	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	71	/// </item>
	72	/// </list></returns>
	73	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	74	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	75	/// This feature is, however, officiallny not supported.</para></remarks>
	76	public unsafe static ILRetArray<double> apply(Func<double, double, double> func, ILInArray<double> A, ILInArray<double> B) {
	77	using (ILScope.Enter(A, B)) {
	78	int outLen;
	79	BinOpItMode mode;
	80
	81	double[] retArr;
	82
	83	double[] arrA = A.GetArrayForRead();
	84
	85	double[] arrB = B.GetArrayForRead();
	86	ILSize outDims;
	87	#region determine operation mode
	88	if (A.IsScalar) {
	89	outDims = B.Size;
	90	if (B.IsScalar) {
	91
	92	return new ILRetArray<double>(new double[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	93	} else if (B.IsEmpty) {
	94	return ILRetArray<double>.empty(outDims);
	95	} else {
	96	outLen = outDims.NumberOfElements;
	97	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	98	retArr = ILMemoryPool.Pool.New<double>(outLen);
	99	mode = BinOpItMode.SAN;
	100	} else {
	101	mode = BinOpItMode.SAI;
	102	}
	103	}
	104	} else {
	105	outDims = A.Size;
	106	if (B.IsScalar) {
	107	if (A.IsEmpty) {
	108	return ILRetArray<double>.empty(A.Size);
	109	}
	110	outLen = A.S.NumberOfElements;
	111	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	112	retArr = ILMemoryPool.Pool.New<double>(outLen);
	113	mode = BinOpItMode.ASN;
	114	} else {
	115	mode = BinOpItMode.ASI;
	116	}
	117	} else {
	118	// array + array
	119	if (!A.Size.IsSameSize(B.Size)) {
	120	return applyEx(func,A,B);
	121	}
	122	outLen = A.S.NumberOfElements;
	123	if (A.TryGetStorage4InplaceOp(out retArr))
	124	mode = BinOpItMode.AAIA;
	125	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	126	mode = BinOpItMode.AAIB;
	127	} else {
	128	retArr = ILMemoryPool.Pool.New<double>(outLen);
	129	mode = BinOpItMode.AAN;
	130	}
	131	}
	132	}
	133	#endregion
	134	ILDenseStorage<double> retStorage = new ILDenseStorage<double>(retArr, outDims);
	135	int i = 0, workerCount = 1;
	136	Action<object> worker = data => {
	137	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	138	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	139
	140	double* cLast, cp = (double*)range.Item5 + range.Item1;
	141
	142	double scalar;
	143	cLast = cp + range.Item2;
	144	#region loops
	145	switch (mode) {
	146	case BinOpItMode.AAIA:
	147
	148	double* bp = ((double*)range.Item4 + range.Item1);
	149	while (cp < cLast) {
	150
	151	cp = func(cp, *bp++);
	152	cp++;
	153	}
	154	break;
	155	case BinOpItMode.AAIB:
	156
	157	double* ap = ((double*)range.Item3 + range.Item1);
	158	while (cp < cLast) {
	159
	160	cp = func(ap++, *cp);
	161	cp++;
	162
	163	}
	164	//ap = ((double*)range.Item3 + range.Item1);
	165	//for (int i2 = range.Item2; i2-- > 0; ) {
	166	// (cp + i2) = (ap + i2) - *(cp + i2);
	167	//}
	168	//int ie = range.Item1 + range.Item2-1;
	169	//double[] locRetArr = retArr;
	170	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	171	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	172	// if (i2 >= ie) break;
	173	//}
	174
	175	break;
	176	case BinOpItMode.AAN:
	177	ap = ((double*)range.Item3 + range.Item1);
	178	bp = ((double*)range.Item4 + range.Item1);
	179	while (cp < cLast) {
	180
	181	cp++ = func(ap++, *bp++);
	182	}
	183	break;
	184	case BinOpItMode.ASI:
	185	scalar = ((double)range.Item4);
	186	while (cp < cLast) {
	187
	188	cp = func(cp, scalar);
	189	cp++;
	190	}
	191	break;
	192	case BinOpItMode.ASN:
	193	ap = ((double*)range.Item3 + range.Item1);
	194	scalar = ((double)range.Item4);
	195	while (cp < cLast) {
	196
	197	cp++ = func(ap++, scalar);
	198	}
	199	break;
	200	case BinOpItMode.SAI:
	201	scalar = ((double)range.Item3);
	202	while (cp < cLast) {
	203
	204	cp = func(scalar, cp);
	205	cp++;
	206	}
	207	break;
	208	case BinOpItMode.SAN:
	209	scalar = ((double)range.Item3);
	210	bp = ((double*)range.Item4 + range.Item1);
	211	while (cp < cLast) {
	212
	213	cp++ = func(scalar, bp++);
	214	}
	215	break;
	216	default:
	217	break;
	218	}
	219	#endregion
	220	System.Threading.Interlocked.Decrement(ref workerCount);
	221	//retStorage.PendingEvents.Signal();
	222	};
	223
	224	#region do the work
	225	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	226	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	227	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	228	workItemLength = outLen / workItemCount;
	229	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	230	} else {
	231	workItemLength = outLen / 2;
	232	workItemCount = 2;
	233	}
	234	} else {
	235	workItemLength = outLen;
	236	workItemCount = 1;
	237	}
	238
	239	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	240
	241	fixed ( double* arrAP = arrA)
	242	fixed ( double* arrBP = arrB)
	243	fixed ( double* retArrP = retArr) {
	244
	245	for (; i < workItemCount - 1; i++) {
	246	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	247	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	248	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	249	System.Threading.Interlocked.Increment(ref workerCount);
	250	ILThreadPool.QueueUserWorkItem(i, worker, range);
	251	}
	252	// the last (or may the only) chunk is done right here
	253	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	254	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	255	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	256
	257	System.Threading.SpinWait.SpinUntil(() => {
	258	return workerCount <= 0;
	259	});
	260	//while (workerCount > 0) ;
	261	}
	262
	263	#endregion
	264	return new ILRetArray<double>(retStorage);
	265	}
	266	}
	267
	268	private static unsafe ILRetArray<double> applyEx(Func<double, double, double> applyFunc, ILInArray<double> A, ILInArray<double> B) {
	269	#region parameter checking
	270	if (isnull(A) \|\| isnull(B))
	271	return empty<double>(ILSize.Empty00);
	272	if (A.IsEmpty) {
	273	return empty<double>(B.S);
	274	} else if (B.IsEmpty) {
	275	return empty<double>(A.S);
	276	}
	277	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	278	// return add(A,B);
	279	int dim = -1;
	280	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	281	if (A.S[_L] != B.S[_L]) {
	282	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	283	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	284	}
	285	dim = _L;
	286	}
	287	}
	288	if (dim > 1)
	289	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	290	#endregion
	291
	292	#region parameter preparation
	293
	294
	295	double[] retArr;
	296
	297
	298	double[] arrA = A.GetArrayForRead();
	299
	300
	301	double[] arrB = B.GetArrayForRead();
	302	ILSize outDims;
	303	BinOptItExMode mode;
	304	int arrInc = 0;
	305	int arrStepInc = 0;
	306	int dimLen = 0;
	307	if (A.IsVector) {
	308	outDims = B.S;
	309	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	310	retArr = ILMemoryPool.Pool.New<double>(outDims.NumberOfElements);
	311	mode = BinOptItExMode.VAN;
	312	} else {
	313	mode = BinOptItExMode.VAI;
	314	}
	315	dimLen = A.Length;
	316	} else if (B.IsVector) {
	317	outDims = A.S;
	318	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	319	retArr = ILMemoryPool.Pool.New<double>(outDims.NumberOfElements);
	320	mode = BinOptItExMode.AVN;
	321	} else {
	322	mode = BinOptItExMode.AVI;
	323	}
	324	dimLen = B.Length;
	325	} else {
	326	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	327	}
	328	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	329	arrStepInc = outDims.SequentialIndexDistance(dim);
	330	#endregion
	331
	332	#region worker loops definition
	333	ILDenseStorage<double> retStorage = new ILDenseStorage<double>(retArr, outDims);
	334	int workerCount = 1;
	335	Action<object> worker = data => {
	336	// expects: iStart, iLen, ap, bp, cp
	337	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	338	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	339
	340	double* ap;
	341
	342	double* bp;
	343
	344	double* cp;
	345	switch (mode) {
	346	case BinOptItExMode.VAN:
	347	for (int s = 0; s < range.Item2; s++) {
	348	ap = (double*)range.Item3;
	349	bp = (double)range.Item4 + range.Item1 + s arrStepInc; ;
	350	cp = (double)range.Item5 + range.Item1 + s arrStepInc;
	351	for (int l = 0; l < dimLen; l++) {
	352
	353	cp = applyFunc(ap, *bp);
	354	ap++;
	355	bp += arrInc;
	356	cp += arrInc;
	357	}
	358	}
	359	break;
	360	case BinOptItExMode.VAI:
	361	for (int s = 0; s < range.Item2; s++) {
	362	ap = (double*)range.Item3;
	363	cp = (double)range.Item5 + range.Item1 + s arrStepInc;
	364	for (int l = 0; l < dimLen; l++) {
	365
	366	cp = applyFunc(ap, *cp);
	367	ap++;
	368	cp += arrInc;
	369	}
	370	}
	371	break;
	372	case BinOptItExMode.AVN:
	373	for (int s = 0; s < range.Item2; s++) {
	374	ap = (double)range.Item3 + range.Item1 + s arrStepInc;
	375	bp = (double*)range.Item4;
	376	cp = (double)range.Item5 + range.Item1 + s arrStepInc;
	377	for (int l = 0; l < dimLen; l++) {
	378
	379	cp = applyFunc(ap, *bp);
	380	ap += arrInc;
	381	bp++;
	382	cp += arrInc;
	383	}
	384	}
	385	break;
	386	case BinOptItExMode.AVI:
	387	for (int s = 0; s < range.Item2; s++) {
	388	bp = (double*)range.Item4;
	389	cp = (double)range.Item5 + range.Item1 + s arrStepInc;
	390	for (int l = 0; l < dimLen; l++) {
	391
	392	cp = applyFunc(cp, *bp);
	393	bp++;
	394	cp += arrInc;
	395	}
	396	}
	397	break;
	398	}
	399	System.Threading.Interlocked.Decrement(ref workerCount);
	400	};
	401	#endregion
	402
	403	#region work distribution
	404	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	405	int outLen = outDims.NumberOfElements;
	406	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	407	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	408	workItemLength = outLen / dimLen / workItemCount;
	409	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	410	} else {
	411	workItemLength = outLen / dimLen / 2;
	412	workItemCount = 2;
	413	}
	414	} else {
	415	workItemLength = outLen / dimLen;
	416	workItemCount = 1;
	417	}
	418
	419	fixed (double* arrAP = arrA)
	420	fixed (double* arrBP = arrB)
	421	fixed (double* retArrP = retArr) {
	422
	423	for (; i < workItemCount - 1; i++) {
	424	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	425	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	426	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	427	System.Threading.Interlocked.Increment(ref workerCount);
	428	ILThreadPool.QueueUserWorkItem(i, worker, range);
	429	}
	430	// the last (or may the only) chunk is done right here
	431	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	432	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	433	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	434
	435	ILThreadPool.Wait4Workers(ref workerCount);
	436	}
	437	#endregion
	438
	439	return new ILRetArray<double>(retStorage);
	440	}
	441
	442	#region HYCALPER AUTO GENERATED CODE
	443
	444	/// <summary>Apply an arbitrary function to two arrays</summary>
	445	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	446	/// <param name="A">Input array A</param>
	447	/// <param name="B">Input array B</param>
	448	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	449	/// <item>
	450	/// <term>size(A) == size(B)</term>
	451	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	452	/// </item>
	453	/// <item>
	454	/// <term>isscalar(A) \|\| isscalar(B)</term>
	455	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	456	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	457	/// </item>
	458	/// <item>
	459	/// <term>All other cases</term>
	460	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	461	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	462	/// </item>
	463	/// </list></returns>
	464	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	465	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	466	/// This feature is, however, officiallny not supported.</para></remarks>
	467	public unsafe static ILRetArray<float> apply(Func<float, float, float> func, ILInArray<float> A, ILInArray<float> B) {
	468	using (ILScope.Enter(A, B)) {
	469	int outLen;
	470	BinOpItMode mode;
	471
	472	float[] retArr;
	473
	474	float[] arrA = A.GetArrayForRead();
	475
	476	float[] arrB = B.GetArrayForRead();
	477	ILSize outDims;
	478	#region determine operation mode
	479	if (A.IsScalar) {
	480	outDims = B.Size;
	481	if (B.IsScalar) {
	482
	483	return new ILRetArray<float>(new float[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	484	} else if (B.IsEmpty) {
	485	return ILRetArray<float>.empty(outDims);
	486	} else {
	487	outLen = outDims.NumberOfElements;
	488	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	489	retArr = ILMemoryPool.Pool.New<float>(outLen);
	490	mode = BinOpItMode.SAN;
	491	} else {
	492	mode = BinOpItMode.SAI;
	493	}
	494	}
	495	} else {
	496	outDims = A.Size;
	497	if (B.IsScalar) {
	498	if (A.IsEmpty) {
	499	return ILRetArray<float>.empty(A.Size);
	500	}
	501	outLen = A.S.NumberOfElements;
	502	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	503	retArr = ILMemoryPool.Pool.New<float>(outLen);
	504	mode = BinOpItMode.ASN;
	505	} else {
	506	mode = BinOpItMode.ASI;
	507	}
	508	} else {
	509	// array + array
	510	if (!A.Size.IsSameSize(B.Size)) {
	511	return applyEx(func,A,B);
	512	}
	513	outLen = A.S.NumberOfElements;
	514	if (A.TryGetStorage4InplaceOp(out retArr))
	515	mode = BinOpItMode.AAIA;
	516	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	517	mode = BinOpItMode.AAIB;
	518	} else {
	519	retArr = ILMemoryPool.Pool.New<float>(outLen);
	520	mode = BinOpItMode.AAN;
	521	}
	522	}
	523	}
	524	#endregion
	525	ILDenseStorage<float> retStorage = new ILDenseStorage<float>(retArr, outDims);
	526	int i = 0, workerCount = 1;
	527	Action<object> worker = data => {
	528	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	529	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	530
	531	float* cLast, cp = (float*)range.Item5 + range.Item1;
	532
	533	float scalar;
	534	cLast = cp + range.Item2;
	535	#region loops
	536	switch (mode) {
	537	case BinOpItMode.AAIA:
	538
	539	float* bp = ((float*)range.Item4 + range.Item1);
	540	while (cp < cLast) {
	541
	542	cp = func(cp, *bp++);
	543	cp++;
	544	}
	545	break;
	546	case BinOpItMode.AAIB:
	547
	548	float* ap = ((float*)range.Item3 + range.Item1);
	549	while (cp < cLast) {
	550
	551	cp = func(ap++, *cp);
	552	cp++;
	553
	554	}
	555	//ap = ((double*)range.Item3 + range.Item1);
	556	//for (int i2 = range.Item2; i2-- > 0; ) {
	557	// (cp + i2) = (ap + i2) - *(cp + i2);
	558	//}
	559	//int ie = range.Item1 + range.Item2-1;
	560	//double[] locRetArr = retArr;
	561	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	562	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	563	// if (i2 >= ie) break;
	564	//}
	565
	566	break;
	567	case BinOpItMode.AAN:
	568	ap = ((float*)range.Item3 + range.Item1);
	569	bp = ((float*)range.Item4 + range.Item1);
	570	while (cp < cLast) {
	571
	572	cp++ = func(ap++, *bp++);
	573	}
	574	break;
	575	case BinOpItMode.ASI:
	576	scalar = ((float)range.Item4);
	577	while (cp < cLast) {
	578
	579	cp = func(cp, scalar);
	580	cp++;
	581	}
	582	break;
	583	case BinOpItMode.ASN:
	584	ap = ((float*)range.Item3 + range.Item1);
	585	scalar = ((float)range.Item4);
	586	while (cp < cLast) {
	587
	588	cp++ = func(ap++, scalar);
	589	}
	590	break;
	591	case BinOpItMode.SAI:
	592	scalar = ((float)range.Item3);
	593	while (cp < cLast) {
	594
	595	cp = func(scalar, cp);
	596	cp++;
	597	}
	598	break;
	599	case BinOpItMode.SAN:
	600	scalar = ((float)range.Item3);
	601	bp = ((float*)range.Item4 + range.Item1);
	602	while (cp < cLast) {
	603
	604	cp++ = func(scalar, bp++);
	605	}
	606	break;
	607	default:
	608	break;
	609	}
	610	#endregion
	611	System.Threading.Interlocked.Decrement(ref workerCount);
	612	//retStorage.PendingEvents.Signal();
	613	};
	614
	615	#region do the work
	616	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	617	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	618	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	619	workItemLength = outLen / workItemCount;
	620	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	621	} else {
	622	workItemLength = outLen / 2;
	623	workItemCount = 2;
	624	}
	625	} else {
	626	workItemLength = outLen;
	627	workItemCount = 1;
	628	}
	629
	630	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	631
	632	fixed ( float* arrAP = arrA)
	633	fixed ( float* arrBP = arrB)
	634	fixed ( float* retArrP = retArr) {
	635
	636	for (; i < workItemCount - 1; i++) {
	637	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	638	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	639	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	640	System.Threading.Interlocked.Increment(ref workerCount);
	641	ILThreadPool.QueueUserWorkItem(i, worker, range);
	642	}
	643	// the last (or may the only) chunk is done right here
	644	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	645	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	646	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	647
	648	System.Threading.SpinWait.SpinUntil(() => {
	649	return workerCount <= 0;
	650	});
	651	//while (workerCount > 0) ;
	652	}
	653
	654	#endregion
	655	return new ILRetArray<float>(retStorage);
	656	}
	657	}
	658
	659	private static unsafe ILRetArray<float> applyEx(Func<float, float, float> applyFunc, ILInArray<float> A, ILInArray<float> B) {
	660	#region parameter checking
	661	if (isnull(A) \|\| isnull(B))
	662	return empty<float>(ILSize.Empty00);
	663	if (A.IsEmpty) {
	664	return empty<float>(B.S);
	665	} else if (B.IsEmpty) {
	666	return empty<float>(A.S);
	667	}
	668	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	669	// return add(A,B);
	670	int dim = -1;
	671	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	672	if (A.S[_L] != B.S[_L]) {
	673	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	674	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	675	}
	676	dim = _L;
	677	}
	678	}
	679	if (dim > 1)
	680	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	681	#endregion
	682
	683	#region parameter preparation
	684
	685
	686	float[] retArr;
	687
	688
	689	float[] arrA = A.GetArrayForRead();
	690
	691
	692	float[] arrB = B.GetArrayForRead();
	693	ILSize outDims;
	694	BinOptItExMode mode;
	695	int arrInc = 0;
	696	int arrStepInc = 0;
	697	int dimLen = 0;
	698	if (A.IsVector) {
	699	outDims = B.S;
	700	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	701	retArr = ILMemoryPool.Pool.New<float>(outDims.NumberOfElements);
	702	mode = BinOptItExMode.VAN;
	703	} else {
	704	mode = BinOptItExMode.VAI;
	705	}
	706	dimLen = A.Length;
	707	} else if (B.IsVector) {
	708	outDims = A.S;
	709	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	710	retArr = ILMemoryPool.Pool.New<float>(outDims.NumberOfElements);
	711	mode = BinOptItExMode.AVN;
	712	} else {
	713	mode = BinOptItExMode.AVI;
	714	}
	715	dimLen = B.Length;
	716	} else {
	717	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	718	}
	719	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	720	arrStepInc = outDims.SequentialIndexDistance(dim);
	721	#endregion
	722
	723	#region worker loops definition
	724	ILDenseStorage<float> retStorage = new ILDenseStorage<float>(retArr, outDims);
	725	int workerCount = 1;
	726	Action<object> worker = data => {
	727	// expects: iStart, iLen, ap, bp, cp
	728	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	729	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	730
	731	float* ap;
	732
	733	float* bp;
	734
	735	float* cp;
	736	switch (mode) {
	737	case BinOptItExMode.VAN:
	738	for (int s = 0; s < range.Item2; s++) {
	739	ap = (float*)range.Item3;
	740	bp = (float)range.Item4 + range.Item1 + s arrStepInc; ;
	741	cp = (float)range.Item5 + range.Item1 + s arrStepInc;
	742	for (int l = 0; l < dimLen; l++) {
	743
	744	cp = applyFunc(ap, *bp);
	745	ap++;
	746	bp += arrInc;
	747	cp += arrInc;
	748	}
	749	}
	750	break;
	751	case BinOptItExMode.VAI:
	752	for (int s = 0; s < range.Item2; s++) {
	753	ap = (float*)range.Item3;
	754	cp = (float)range.Item5 + range.Item1 + s arrStepInc;
	755	for (int l = 0; l < dimLen; l++) {
	756
	757	cp = applyFunc(ap, *cp);
	758	ap++;
	759	cp += arrInc;
	760	}
	761	}
	762	break;
	763	case BinOptItExMode.AVN:
	764	for (int s = 0; s < range.Item2; s++) {
	765	ap = (float)range.Item3 + range.Item1 + s arrStepInc;
	766	bp = (float*)range.Item4;
	767	cp = (float)range.Item5 + range.Item1 + s arrStepInc;
	768	for (int l = 0; l < dimLen; l++) {
	769
	770	cp = applyFunc(ap, *bp);
	771	ap += arrInc;
	772	bp++;
	773	cp += arrInc;
	774	}
	775	}
	776	break;
	777	case BinOptItExMode.AVI:
	778	for (int s = 0; s < range.Item2; s++) {
	779	bp = (float*)range.Item4;
	780	cp = (float)range.Item5 + range.Item1 + s arrStepInc;
	781	for (int l = 0; l < dimLen; l++) {
	782
	783	cp = applyFunc(cp, *bp);
	784	bp++;
	785	cp += arrInc;
	786	}
	787	}
	788	break;
	789	}
	790	System.Threading.Interlocked.Decrement(ref workerCount);
	791	};
	792	#endregion
	793
	794	#region work distribution
	795	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	796	int outLen = outDims.NumberOfElements;
	797	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	798	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	799	workItemLength = outLen / dimLen / workItemCount;
	800	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	801	} else {
	802	workItemLength = outLen / dimLen / 2;
	803	workItemCount = 2;
	804	}
	805	} else {
	806	workItemLength = outLen / dimLen;
	807	workItemCount = 1;
	808	}
	809
	810	fixed (float* arrAP = arrA)
	811	fixed (float* arrBP = arrB)
	812	fixed (float* retArrP = retArr) {
	813
	814	for (; i < workItemCount - 1; i++) {
	815	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	816	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	817	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	818	System.Threading.Interlocked.Increment(ref workerCount);
	819	ILThreadPool.QueueUserWorkItem(i, worker, range);
	820	}
	821	// the last (or may the only) chunk is done right here
	822	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	823	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	824	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	825
	826	ILThreadPool.Wait4Workers(ref workerCount);
	827	}
	828	#endregion
	829
	830	return new ILRetArray<float>(retStorage);
	831	}
	832	/// <summary>Apply an arbitrary function to two arrays</summary>
	833	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	834	/// <param name="A">Input array A</param>
	835	/// <param name="B">Input array B</param>
	836	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	837	/// <item>
	838	/// <term>size(A) == size(B)</term>
	839	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	840	/// </item>
	841	/// <item>
	842	/// <term>isscalar(A) \|\| isscalar(B)</term>
	843	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	844	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	845	/// </item>
	846	/// <item>
	847	/// <term>All other cases</term>
	848	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	849	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	850	/// </item>
	851	/// </list></returns>
	852	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	853	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	854	/// This feature is, however, officiallny not supported.</para></remarks>
	855	public unsafe static ILRetArray<fcomplex> apply(Func<fcomplex, fcomplex, fcomplex> func, ILInArray<fcomplex> A, ILInArray<fcomplex> B) {
	856	using (ILScope.Enter(A, B)) {
	857	int outLen;
	858	BinOpItMode mode;
	859
	860	fcomplex[] retArr;
	861
	862	fcomplex[] arrA = A.GetArrayForRead();
	863
	864	fcomplex[] arrB = B.GetArrayForRead();
	865	ILSize outDims;
	866	#region determine operation mode
	867	if (A.IsScalar) {
	868	outDims = B.Size;
	869	if (B.IsScalar) {
	870
	871	return new ILRetArray<fcomplex>(new fcomplex[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	872	} else if (B.IsEmpty) {
	873	return ILRetArray<fcomplex>.empty(outDims);
	874	} else {
	875	outLen = outDims.NumberOfElements;
	876	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	877	retArr = ILMemoryPool.Pool.New<fcomplex>(outLen);
	878	mode = BinOpItMode.SAN;
	879	} else {
	880	mode = BinOpItMode.SAI;
	881	}
	882	}
	883	} else {
	884	outDims = A.Size;
	885	if (B.IsScalar) {
	886	if (A.IsEmpty) {
	887	return ILRetArray<fcomplex>.empty(A.Size);
	888	}
	889	outLen = A.S.NumberOfElements;
	890	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	891	retArr = ILMemoryPool.Pool.New<fcomplex>(outLen);
	892	mode = BinOpItMode.ASN;
	893	} else {
	894	mode = BinOpItMode.ASI;
	895	}
	896	} else {
	897	// array + array
	898	if (!A.Size.IsSameSize(B.Size)) {
	899	return applyEx(func,A,B);
	900	}
	901	outLen = A.S.NumberOfElements;
	902	if (A.TryGetStorage4InplaceOp(out retArr))
	903	mode = BinOpItMode.AAIA;
	904	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	905	mode = BinOpItMode.AAIB;
	906	} else {
	907	retArr = ILMemoryPool.Pool.New<fcomplex>(outLen);
	908	mode = BinOpItMode.AAN;
	909	}
	910	}
	911	}
	912	#endregion
	913	ILDenseStorage<fcomplex> retStorage = new ILDenseStorage<fcomplex>(retArr, outDims);
	914	int i = 0, workerCount = 1;
	915	Action<object> worker = data => {
	916	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	917	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	918
	919	fcomplex* cLast, cp = (fcomplex*)range.Item5 + range.Item1;
	920
	921	fcomplex scalar;
	922	cLast = cp + range.Item2;
	923	#region loops
	924	switch (mode) {
	925	case BinOpItMode.AAIA:
	926
	927	fcomplex* bp = ((fcomplex*)range.Item4 + range.Item1);
	928	while (cp < cLast) {
	929
	930	cp = func(cp, *bp++);
	931	cp++;
	932	}
	933	break;
	934	case BinOpItMode.AAIB:
	935
	936	fcomplex* ap = ((fcomplex*)range.Item3 + range.Item1);
	937	while (cp < cLast) {
	938
	939	cp = func(ap++, *cp);
	940	cp++;
	941
	942	}
	943	//ap = ((double*)range.Item3 + range.Item1);
	944	//for (int i2 = range.Item2; i2-- > 0; ) {
	945	// (cp + i2) = (ap + i2) - *(cp + i2);
	946	//}
	947	//int ie = range.Item1 + range.Item2-1;
	948	//double[] locRetArr = retArr;
	949	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	950	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	951	// if (i2 >= ie) break;
	952	//}
	953
	954	break;
	955	case BinOpItMode.AAN:
	956	ap = ((fcomplex*)range.Item3 + range.Item1);
	957	bp = ((fcomplex*)range.Item4 + range.Item1);
	958	while (cp < cLast) {
	959
	960	cp++ = func(ap++, *bp++);
	961	}
	962	break;
	963	case BinOpItMode.ASI:
	964	scalar = ((fcomplex)range.Item4);
	965	while (cp < cLast) {
	966
	967	cp = func(cp, scalar);
	968	cp++;
	969	}
	970	break;
	971	case BinOpItMode.ASN:
	972	ap = ((fcomplex*)range.Item3 + range.Item1);
	973	scalar = ((fcomplex)range.Item4);
	974	while (cp < cLast) {
	975
	976	cp++ = func(ap++, scalar);
	977	}
	978	break;
	979	case BinOpItMode.SAI:
	980	scalar = ((fcomplex)range.Item3);
	981	while (cp < cLast) {
	982
	983	cp = func(scalar, cp);
	984	cp++;
	985	}
	986	break;
	987	case BinOpItMode.SAN:
	988	scalar = ((fcomplex)range.Item3);
	989	bp = ((fcomplex*)range.Item4 + range.Item1);
	990	while (cp < cLast) {
	991
	992	cp++ = func(scalar, bp++);
	993	}
	994	break;
	995	default:
	996	break;
	997	}
	998	#endregion
	999	System.Threading.Interlocked.Decrement(ref workerCount);
	1000	//retStorage.PendingEvents.Signal();
	1001	};
	1002
	1003	#region do the work
	1004	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1005	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	1006	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1007	workItemLength = outLen / workItemCount;
	1008	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1009	} else {
	1010	workItemLength = outLen / 2;
	1011	workItemCount = 2;
	1012	}
	1013	} else {
	1014	workItemLength = outLen;
	1015	workItemCount = 1;
	1016	}
	1017
	1018	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	1019
	1020	fixed ( fcomplex* arrAP = arrA)
	1021	fixed ( fcomplex* arrBP = arrB)
	1022	fixed ( fcomplex* retArrP = retArr) {
	1023
	1024	for (; i < workItemCount - 1; i++) {
	1025	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	1026	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1027	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	1028	System.Threading.Interlocked.Increment(ref workerCount);
	1029	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1030	}
	1031	// the last (or may the only) chunk is done right here
	1032	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1033	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1034	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	1035
	1036	System.Threading.SpinWait.SpinUntil(() => {
	1037	return workerCount <= 0;
	1038	});
	1039	//while (workerCount > 0) ;
	1040	}
	1041
	1042	#endregion
	1043	return new ILRetArray<fcomplex>(retStorage);
	1044	}
	1045	}
	1046
	1047	private static unsafe ILRetArray<fcomplex> applyEx(Func<fcomplex, fcomplex, fcomplex> applyFunc, ILInArray<fcomplex> A, ILInArray<fcomplex> B) {
	1048	#region parameter checking
	1049	if (isnull(A) \|\| isnull(B))
	1050	return empty<fcomplex>(ILSize.Empty00);
	1051	if (A.IsEmpty) {
	1052	return empty<fcomplex>(B.S);
	1053	} else if (B.IsEmpty) {
	1054	return empty<fcomplex>(A.S);
	1055	}
	1056	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	1057	// return add(A,B);
	1058	int dim = -1;
	1059	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	1060	if (A.S[_L] != B.S[_L]) {
	1061	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	1062	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	1063	}
	1064	dim = _L;
	1065	}
	1066	}
	1067	if (dim > 1)
	1068	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	1069	#endregion
	1070
	1071	#region parameter preparation
	1072
	1073
	1074	fcomplex[] retArr;
	1075
	1076
	1077	fcomplex[] arrA = A.GetArrayForRead();
	1078
	1079
	1080	fcomplex[] arrB = B.GetArrayForRead();
	1081	ILSize outDims;
	1082	BinOptItExMode mode;
	1083	int arrInc = 0;
	1084	int arrStepInc = 0;
	1085	int dimLen = 0;
	1086	if (A.IsVector) {
	1087	outDims = B.S;
	1088	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	1089	retArr = ILMemoryPool.Pool.New<fcomplex>(outDims.NumberOfElements);
	1090	mode = BinOptItExMode.VAN;
	1091	} else {
	1092	mode = BinOptItExMode.VAI;
	1093	}
	1094	dimLen = A.Length;
	1095	} else if (B.IsVector) {
	1096	outDims = A.S;
	1097	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	1098	retArr = ILMemoryPool.Pool.New<fcomplex>(outDims.NumberOfElements);
	1099	mode = BinOptItExMode.AVN;
	1100	} else {
	1101	mode = BinOptItExMode.AVI;
	1102	}
	1103	dimLen = B.Length;
	1104	} else {
	1105	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	1106	}
	1107	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	1108	arrStepInc = outDims.SequentialIndexDistance(dim);
	1109	#endregion
	1110
	1111	#region worker loops definition
	1112	ILDenseStorage<fcomplex> retStorage = new ILDenseStorage<fcomplex>(retArr, outDims);
	1113	int workerCount = 1;
	1114	Action<object> worker = data => {
	1115	// expects: iStart, iLen, ap, bp, cp
	1116	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	1117	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	1118
	1119	fcomplex* ap;
	1120
	1121	fcomplex* bp;
	1122
	1123	fcomplex* cp;
	1124	switch (mode) {
	1125	case BinOptItExMode.VAN:
	1126	for (int s = 0; s < range.Item2; s++) {
	1127	ap = (fcomplex*)range.Item3;
	1128	bp = (fcomplex)range.Item4 + range.Item1 + s arrStepInc; ;
	1129	cp = (fcomplex)range.Item5 + range.Item1 + s arrStepInc;
	1130	for (int l = 0; l < dimLen; l++) {
	1131
	1132	cp = applyFunc(ap, *bp);
	1133	ap++;
	1134	bp += arrInc;
	1135	cp += arrInc;
	1136	}
	1137	}
	1138	break;
	1139	case BinOptItExMode.VAI:
	1140	for (int s = 0; s < range.Item2; s++) {
	1141	ap = (fcomplex*)range.Item3;
	1142	cp = (fcomplex)range.Item5 + range.Item1 + s arrStepInc;
	1143	for (int l = 0; l < dimLen; l++) {
	1144
	1145	cp = applyFunc(ap, *cp);
	1146	ap++;
	1147	cp += arrInc;
	1148	}
	1149	}
	1150	break;
	1151	case BinOptItExMode.AVN:
	1152	for (int s = 0; s < range.Item2; s++) {
	1153	ap = (fcomplex)range.Item3 + range.Item1 + s arrStepInc;
	1154	bp = (fcomplex*)range.Item4;
	1155	cp = (fcomplex)range.Item5 + range.Item1 + s arrStepInc;
	1156	for (int l = 0; l < dimLen; l++) {
	1157
	1158	cp = applyFunc(ap, *bp);
	1159	ap += arrInc;
	1160	bp++;
	1161	cp += arrInc;
	1162	}
	1163	}
	1164	break;
	1165	case BinOptItExMode.AVI:
	1166	for (int s = 0; s < range.Item2; s++) {
	1167	bp = (fcomplex*)range.Item4;
	1168	cp = (fcomplex)range.Item5 + range.Item1 + s arrStepInc;
	1169	for (int l = 0; l < dimLen; l++) {
	1170
	1171	cp = applyFunc(cp, *bp);
	1172	bp++;
	1173	cp += arrInc;
	1174	}
	1175	}
	1176	break;
	1177	}
	1178	System.Threading.Interlocked.Decrement(ref workerCount);
	1179	};
	1180	#endregion
	1181
	1182	#region work distribution
	1183	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1184	int outLen = outDims.NumberOfElements;
	1185	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	1186	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1187	workItemLength = outLen / dimLen / workItemCount;
	1188	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1189	} else {
	1190	workItemLength = outLen / dimLen / 2;
	1191	workItemCount = 2;
	1192	}
	1193	} else {
	1194	workItemLength = outLen / dimLen;
	1195	workItemCount = 1;
	1196	}
	1197
	1198	fixed (fcomplex* arrAP = arrA)
	1199	fixed (fcomplex* arrBP = arrB)
	1200	fixed (fcomplex* retArrP = retArr) {
	1201
	1202	for (; i < workItemCount - 1; i++) {
	1203	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	1204	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1205	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	1206	System.Threading.Interlocked.Increment(ref workerCount);
	1207	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1208	}
	1209	// the last (or may the only) chunk is done right here
	1210	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1211	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1212	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	1213
	1214	ILThreadPool.Wait4Workers(ref workerCount);
	1215	}
	1216	#endregion
	1217
	1218	return new ILRetArray<fcomplex>(retStorage);
	1219	}
	1220	/// <summary>Apply an arbitrary function to two arrays</summary>
	1221	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	1222	/// <param name="A">Input array A</param>
	1223	/// <param name="B">Input array B</param>
	1224	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	1225	/// <item>
	1226	/// <term>size(A) == size(B)</term>
	1227	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	1228	/// </item>
	1229	/// <item>
	1230	/// <term>isscalar(A) \|\| isscalar(B)</term>
	1231	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	1232	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	1233	/// </item>
	1234	/// <item>
	1235	/// <term>All other cases</term>
	1236	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	1237	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	1238	/// </item>
	1239	/// </list></returns>
	1240	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	1241	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	1242	/// This feature is, however, officiallny not supported.</para></remarks>
	1243	public unsafe static ILRetArray<complex> apply(Func<complex, complex, complex> func, ILInArray<complex> A, ILInArray<complex> B) {
	1244	using (ILScope.Enter(A, B)) {
	1245	int outLen;
	1246	BinOpItMode mode;
	1247
	1248	complex[] retArr;
	1249
	1250	complex[] arrA = A.GetArrayForRead();
	1251
	1252	complex[] arrB = B.GetArrayForRead();
	1253	ILSize outDims;
	1254	#region determine operation mode
	1255	if (A.IsScalar) {
	1256	outDims = B.Size;
	1257	if (B.IsScalar) {
	1258
	1259	return new ILRetArray<complex>(new complex[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	1260	} else if (B.IsEmpty) {
	1261	return ILRetArray<complex>.empty(outDims);
	1262	} else {
	1263	outLen = outDims.NumberOfElements;
	1264	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	1265	retArr = ILMemoryPool.Pool.New<complex>(outLen);
	1266	mode = BinOpItMode.SAN;
	1267	} else {
	1268	mode = BinOpItMode.SAI;
	1269	}
	1270	}
	1271	} else {
	1272	outDims = A.Size;
	1273	if (B.IsScalar) {
	1274	if (A.IsEmpty) {
	1275	return ILRetArray<complex>.empty(A.Size);
	1276	}
	1277	outLen = A.S.NumberOfElements;
	1278	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	1279	retArr = ILMemoryPool.Pool.New<complex>(outLen);
	1280	mode = BinOpItMode.ASN;
	1281	} else {
	1282	mode = BinOpItMode.ASI;
	1283	}
	1284	} else {
	1285	// array + array
	1286	if (!A.Size.IsSameSize(B.Size)) {
	1287	return applyEx(func,A,B);
	1288	}
	1289	outLen = A.S.NumberOfElements;
	1290	if (A.TryGetStorage4InplaceOp(out retArr))
	1291	mode = BinOpItMode.AAIA;
	1292	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	1293	mode = BinOpItMode.AAIB;
	1294	} else {
	1295	retArr = ILMemoryPool.Pool.New<complex>(outLen);
	1296	mode = BinOpItMode.AAN;
	1297	}
	1298	}
	1299	}
	1300	#endregion
	1301	ILDenseStorage<complex> retStorage = new ILDenseStorage<complex>(retArr, outDims);
	1302	int i = 0, workerCount = 1;
	1303	Action<object> worker = data => {
	1304	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	1305	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	1306
	1307	complex* cLast, cp = (complex*)range.Item5 + range.Item1;
	1308
	1309	complex scalar;
	1310	cLast = cp + range.Item2;
	1311	#region loops
	1312	switch (mode) {
	1313	case BinOpItMode.AAIA:
	1314
	1315	complex* bp = ((complex*)range.Item4 + range.Item1);
	1316	while (cp < cLast) {
	1317
	1318	cp = func(cp, *bp++);
	1319	cp++;
	1320	}
	1321	break;
	1322	case BinOpItMode.AAIB:
	1323
	1324	complex* ap = ((complex*)range.Item3 + range.Item1);
	1325	while (cp < cLast) {
	1326
	1327	cp = func(ap++, *cp);
	1328	cp++;
	1329
	1330	}
	1331	//ap = ((double*)range.Item3 + range.Item1);
	1332	//for (int i2 = range.Item2; i2-- > 0; ) {
	1333	// (cp + i2) = (ap + i2) - *(cp + i2);
	1334	//}
	1335	//int ie = range.Item1 + range.Item2-1;
	1336	//double[] locRetArr = retArr;
	1337	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	1338	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	1339	// if (i2 >= ie) break;
	1340	//}
	1341
	1342	break;
	1343	case BinOpItMode.AAN:
	1344	ap = ((complex*)range.Item3 + range.Item1);
	1345	bp = ((complex*)range.Item4 + range.Item1);
	1346	while (cp < cLast) {
	1347
	1348	cp++ = func(ap++, *bp++);
	1349	}
	1350	break;
	1351	case BinOpItMode.ASI:
	1352	scalar = ((complex)range.Item4);
	1353	while (cp < cLast) {
	1354
	1355	cp = func(cp, scalar);
	1356	cp++;
	1357	}
	1358	break;
	1359	case BinOpItMode.ASN:
	1360	ap = ((complex*)range.Item3 + range.Item1);
	1361	scalar = ((complex)range.Item4);
	1362	while (cp < cLast) {
	1363
	1364	cp++ = func(ap++, scalar);
	1365	}
	1366	break;
	1367	case BinOpItMode.SAI:
	1368	scalar = ((complex)range.Item3);
	1369	while (cp < cLast) {
	1370
	1371	cp = func(scalar, cp);
	1372	cp++;
	1373	}
	1374	break;
	1375	case BinOpItMode.SAN:
	1376	scalar = ((complex)range.Item3);
	1377	bp = ((complex*)range.Item4 + range.Item1);
	1378	while (cp < cLast) {
	1379
	1380	cp++ = func(scalar, bp++);
	1381	}
	1382	break;
	1383	default:
	1384	break;
	1385	}
	1386	#endregion
	1387	System.Threading.Interlocked.Decrement(ref workerCount);
	1388	//retStorage.PendingEvents.Signal();
	1389	};
	1390
	1391	#region do the work
	1392	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1393	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	1394	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1395	workItemLength = outLen / workItemCount;
	1396	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1397	} else {
	1398	workItemLength = outLen / 2;
	1399	workItemCount = 2;
	1400	}
	1401	} else {
	1402	workItemLength = outLen;
	1403	workItemCount = 1;
	1404	}
	1405
	1406	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	1407
	1408	fixed ( complex* arrAP = arrA)
	1409	fixed ( complex* arrBP = arrB)
	1410	fixed ( complex* retArrP = retArr) {
	1411
	1412	for (; i < workItemCount - 1; i++) {
	1413	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	1414	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1415	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	1416	System.Threading.Interlocked.Increment(ref workerCount);
	1417	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1418	}
	1419	// the last (or may the only) chunk is done right here
	1420	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1421	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1422	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	1423
	1424	System.Threading.SpinWait.SpinUntil(() => {
	1425	return workerCount <= 0;
	1426	});
	1427	//while (workerCount > 0) ;
	1428	}
	1429
	1430	#endregion
	1431	return new ILRetArray<complex>(retStorage);
	1432	}
	1433	}
	1434
	1435	private static unsafe ILRetArray<complex> applyEx(Func<complex, complex, complex> applyFunc, ILInArray<complex> A, ILInArray<complex> B) {
	1436	#region parameter checking
	1437	if (isnull(A) \|\| isnull(B))
	1438	return empty<complex>(ILSize.Empty00);
	1439	if (A.IsEmpty) {
	1440	return empty<complex>(B.S);
	1441	} else if (B.IsEmpty) {
	1442	return empty<complex>(A.S);
	1443	}
	1444	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	1445	// return add(A,B);
	1446	int dim = -1;
	1447	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	1448	if (A.S[_L] != B.S[_L]) {
	1449	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	1450	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	1451	}
	1452	dim = _L;
	1453	}
	1454	}
	1455	if (dim > 1)
	1456	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	1457	#endregion
	1458
	1459	#region parameter preparation
	1460
	1461
	1462	complex[] retArr;
	1463
	1464
	1465	complex[] arrA = A.GetArrayForRead();
	1466
	1467
	1468	complex[] arrB = B.GetArrayForRead();
	1469	ILSize outDims;
	1470	BinOptItExMode mode;
	1471	int arrInc = 0;
	1472	int arrStepInc = 0;
	1473	int dimLen = 0;
	1474	if (A.IsVector) {
	1475	outDims = B.S;
	1476	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	1477	retArr = ILMemoryPool.Pool.New<complex>(outDims.NumberOfElements);
	1478	mode = BinOptItExMode.VAN;
	1479	} else {
	1480	mode = BinOptItExMode.VAI;
	1481	}
	1482	dimLen = A.Length;
	1483	} else if (B.IsVector) {
	1484	outDims = A.S;
	1485	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	1486	retArr = ILMemoryPool.Pool.New<complex>(outDims.NumberOfElements);
	1487	mode = BinOptItExMode.AVN;
	1488	} else {
	1489	mode = BinOptItExMode.AVI;
	1490	}
	1491	dimLen = B.Length;
	1492	} else {
	1493	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	1494	}
	1495	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	1496	arrStepInc = outDims.SequentialIndexDistance(dim);
	1497	#endregion
	1498
	1499	#region worker loops definition
	1500	ILDenseStorage<complex> retStorage = new ILDenseStorage<complex>(retArr, outDims);
	1501	int workerCount = 1;
	1502	Action<object> worker = data => {
	1503	// expects: iStart, iLen, ap, bp, cp
	1504	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	1505	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	1506
	1507	complex* ap;
	1508
	1509	complex* bp;
	1510
	1511	complex* cp;
	1512	switch (mode) {
	1513	case BinOptItExMode.VAN:
	1514	for (int s = 0; s < range.Item2; s++) {
	1515	ap = (complex*)range.Item3;
	1516	bp = (complex)range.Item4 + range.Item1 + s arrStepInc; ;
	1517	cp = (complex)range.Item5 + range.Item1 + s arrStepInc;
	1518	for (int l = 0; l < dimLen; l++) {
	1519
	1520	cp = applyFunc(ap, *bp);
	1521	ap++;
	1522	bp += arrInc;
	1523	cp += arrInc;
	1524	}
	1525	}
	1526	break;
	1527	case BinOptItExMode.VAI:
	1528	for (int s = 0; s < range.Item2; s++) {
	1529	ap = (complex*)range.Item3;
	1530	cp = (complex)range.Item5 + range.Item1 + s arrStepInc;
	1531	for (int l = 0; l < dimLen; l++) {
	1532
	1533	cp = applyFunc(ap, *cp);
	1534	ap++;
	1535	cp += arrInc;
	1536	}
	1537	}
	1538	break;
	1539	case BinOptItExMode.AVN:
	1540	for (int s = 0; s < range.Item2; s++) {
	1541	ap = (complex)range.Item3 + range.Item1 + s arrStepInc;
	1542	bp = (complex*)range.Item4;
	1543	cp = (complex)range.Item5 + range.Item1 + s arrStepInc;
	1544	for (int l = 0; l < dimLen; l++) {
	1545
	1546	cp = applyFunc(ap, *bp);
	1547	ap += arrInc;
	1548	bp++;
	1549	cp += arrInc;
	1550	}
	1551	}
	1552	break;
	1553	case BinOptItExMode.AVI:
	1554	for (int s = 0; s < range.Item2; s++) {
	1555	bp = (complex*)range.Item4;
	1556	cp = (complex)range.Item5 + range.Item1 + s arrStepInc;
	1557	for (int l = 0; l < dimLen; l++) {
	1558
	1559	cp = applyFunc(cp, *bp);
	1560	bp++;
	1561	cp += arrInc;
	1562	}
	1563	}
	1564	break;
	1565	}
	1566	System.Threading.Interlocked.Decrement(ref workerCount);
	1567	};
	1568	#endregion
	1569
	1570	#region work distribution
	1571	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1572	int outLen = outDims.NumberOfElements;
	1573	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	1574	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1575	workItemLength = outLen / dimLen / workItemCount;
	1576	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1577	} else {
	1578	workItemLength = outLen / dimLen / 2;
	1579	workItemCount = 2;
	1580	}
	1581	} else {
	1582	workItemLength = outLen / dimLen;
	1583	workItemCount = 1;
	1584	}
	1585
	1586	fixed (complex* arrAP = arrA)
	1587	fixed (complex* arrBP = arrB)
	1588	fixed (complex* retArrP = retArr) {
	1589
	1590	for (; i < workItemCount - 1; i++) {
	1591	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	1592	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1593	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	1594	System.Threading.Interlocked.Increment(ref workerCount);
	1595	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1596	}
	1597	// the last (or may the only) chunk is done right here
	1598	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1599	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1600	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	1601
	1602	ILThreadPool.Wait4Workers(ref workerCount);
	1603	}
	1604	#endregion
	1605
	1606	return new ILRetArray<complex>(retStorage);
	1607	}
	1608	/// <summary>Apply an arbitrary function to two arrays</summary>
	1609	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	1610	/// <param name="A">Input array A</param>
	1611	/// <param name="B">Input array B</param>
	1612	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	1613	/// <item>
	1614	/// <term>size(A) == size(B)</term>
	1615	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	1616	/// </item>
	1617	/// <item>
	1618	/// <term>isscalar(A) \|\| isscalar(B)</term>
	1619	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	1620	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	1621	/// </item>
	1622	/// <item>
	1623	/// <term>All other cases</term>
	1624	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	1625	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	1626	/// </item>
	1627	/// </list></returns>
	1628	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	1629	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	1630	/// This feature is, however, officiallny not supported.</para></remarks>
	1631	public unsafe static ILRetArray<Int64> apply(Func<Int64, Int64, Int64> func, ILInArray<Int64> A, ILInArray<Int64> B) {
	1632	using (ILScope.Enter(A, B)) {
	1633	int outLen;
	1634	BinOpItMode mode;
	1635
	1636	Int64[] retArr;
	1637
	1638	Int64[] arrA = A.GetArrayForRead();
	1639
	1640	Int64[] arrB = B.GetArrayForRead();
	1641	ILSize outDims;
	1642	#region determine operation mode
	1643	if (A.IsScalar) {
	1644	outDims = B.Size;
	1645	if (B.IsScalar) {
	1646
	1647	return new ILRetArray<Int64>(new Int64[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	1648	} else if (B.IsEmpty) {
	1649	return ILRetArray<Int64>.empty(outDims);
	1650	} else {
	1651	outLen = outDims.NumberOfElements;
	1652	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	1653	retArr = ILMemoryPool.Pool.New<Int64>(outLen);
	1654	mode = BinOpItMode.SAN;
	1655	} else {
	1656	mode = BinOpItMode.SAI;
	1657	}
	1658	}
	1659	} else {
	1660	outDims = A.Size;
	1661	if (B.IsScalar) {
	1662	if (A.IsEmpty) {
	1663	return ILRetArray<Int64>.empty(A.Size);
	1664	}
	1665	outLen = A.S.NumberOfElements;
	1666	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	1667	retArr = ILMemoryPool.Pool.New<Int64>(outLen);
	1668	mode = BinOpItMode.ASN;
	1669	} else {
	1670	mode = BinOpItMode.ASI;
	1671	}
	1672	} else {
	1673	// array + array
	1674	if (!A.Size.IsSameSize(B.Size)) {
	1675	return applyEx(func,A,B);
	1676	}
	1677	outLen = A.S.NumberOfElements;
	1678	if (A.TryGetStorage4InplaceOp(out retArr))
	1679	mode = BinOpItMode.AAIA;
	1680	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	1681	mode = BinOpItMode.AAIB;
	1682	} else {
	1683	retArr = ILMemoryPool.Pool.New<Int64>(outLen);
	1684	mode = BinOpItMode.AAN;
	1685	}
	1686	}
	1687	}
	1688	#endregion
	1689	ILDenseStorage<Int64> retStorage = new ILDenseStorage<Int64>(retArr, outDims);
	1690	int i = 0, workerCount = 1;
	1691	Action<object> worker = data => {
	1692	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	1693	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	1694
	1695	Int64* cLast, cp = (Int64*)range.Item5 + range.Item1;
	1696
	1697	Int64 scalar;
	1698	cLast = cp + range.Item2;
	1699	#region loops
	1700	switch (mode) {
	1701	case BinOpItMode.AAIA:
	1702
	1703	Int64* bp = ((Int64*)range.Item4 + range.Item1);
	1704	while (cp < cLast) {
	1705
	1706	cp = func(cp, *bp++);
	1707	cp++;
	1708	}
	1709	break;
	1710	case BinOpItMode.AAIB:
	1711
	1712	Int64* ap = ((Int64*)range.Item3 + range.Item1);
	1713	while (cp < cLast) {
	1714
	1715	cp = func(ap++, *cp);
	1716	cp++;
	1717
	1718	}
	1719	//ap = ((double*)range.Item3 + range.Item1);
	1720	//for (int i2 = range.Item2; i2-- > 0; ) {
	1721	// (cp + i2) = (ap + i2) - *(cp + i2);
	1722	//}
	1723	//int ie = range.Item1 + range.Item2-1;
	1724	//double[] locRetArr = retArr;
	1725	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	1726	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	1727	// if (i2 >= ie) break;
	1728	//}
	1729
	1730	break;
	1731	case BinOpItMode.AAN:
	1732	ap = ((Int64*)range.Item3 + range.Item1);
	1733	bp = ((Int64*)range.Item4 + range.Item1);
	1734	while (cp < cLast) {
	1735
	1736	cp++ = func(ap++, *bp++);
	1737	}
	1738	break;
	1739	case BinOpItMode.ASI:
	1740	scalar = ((Int64)range.Item4);
	1741	while (cp < cLast) {
	1742
	1743	cp = func(cp, scalar);
	1744	cp++;
	1745	}
	1746	break;
	1747	case BinOpItMode.ASN:
	1748	ap = ((Int64*)range.Item3 + range.Item1);
	1749	scalar = ((Int64)range.Item4);
	1750	while (cp < cLast) {
	1751
	1752	cp++ = func(ap++, scalar);
	1753	}
	1754	break;
	1755	case BinOpItMode.SAI:
	1756	scalar = ((Int64)range.Item3);
	1757	while (cp < cLast) {
	1758
	1759	cp = func(scalar, cp);
	1760	cp++;
	1761	}
	1762	break;
	1763	case BinOpItMode.SAN:
	1764	scalar = ((Int64)range.Item3);
	1765	bp = ((Int64*)range.Item4 + range.Item1);
	1766	while (cp < cLast) {
	1767
	1768	cp++ = func(scalar, bp++);
	1769	}
	1770	break;
	1771	default:
	1772	break;
	1773	}
	1774	#endregion
	1775	System.Threading.Interlocked.Decrement(ref workerCount);
	1776	//retStorage.PendingEvents.Signal();
	1777	};
	1778
	1779	#region do the work
	1780	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1781	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	1782	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1783	workItemLength = outLen / workItemCount;
	1784	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1785	} else {
	1786	workItemLength = outLen / 2;
	1787	workItemCount = 2;
	1788	}
	1789	} else {
	1790	workItemLength = outLen;
	1791	workItemCount = 1;
	1792	}
	1793
	1794	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	1795
	1796	fixed ( Int64* arrAP = arrA)
	1797	fixed ( Int64* arrBP = arrB)
	1798	fixed ( Int64* retArrP = retArr) {
	1799
	1800	for (; i < workItemCount - 1; i++) {
	1801	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	1802	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1803	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	1804	System.Threading.Interlocked.Increment(ref workerCount);
	1805	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1806	}
	1807	// the last (or may the only) chunk is done right here
	1808	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1809	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	1810	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	1811
	1812	System.Threading.SpinWait.SpinUntil(() => {
	1813	return workerCount <= 0;
	1814	});
	1815	//while (workerCount > 0) ;
	1816	}
	1817
	1818	#endregion
	1819	return new ILRetArray<Int64>(retStorage);
	1820	}
	1821	}
	1822
	1823	private static unsafe ILRetArray<Int64> applyEx(Func<Int64, Int64, Int64> applyFunc, ILInArray<Int64> A, ILInArray<Int64> B) {
	1824	#region parameter checking
	1825	if (isnull(A) \|\| isnull(B))
	1826	return empty<Int64>(ILSize.Empty00);
	1827	if (A.IsEmpty) {
	1828	return empty<Int64>(B.S);
	1829	} else if (B.IsEmpty) {
	1830	return empty<Int64>(A.S);
	1831	}
	1832	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	1833	// return add(A,B);
	1834	int dim = -1;
	1835	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	1836	if (A.S[_L] != B.S[_L]) {
	1837	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	1838	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	1839	}
	1840	dim = _L;
	1841	}
	1842	}
	1843	if (dim > 1)
	1844	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	1845	#endregion
	1846
	1847	#region parameter preparation
	1848
	1849
	1850	Int64[] retArr;
	1851
	1852
	1853	Int64[] arrA = A.GetArrayForRead();
	1854
	1855
	1856	Int64[] arrB = B.GetArrayForRead();
	1857	ILSize outDims;
	1858	BinOptItExMode mode;
	1859	int arrInc = 0;
	1860	int arrStepInc = 0;
	1861	int dimLen = 0;
	1862	if (A.IsVector) {
	1863	outDims = B.S;
	1864	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	1865	retArr = ILMemoryPool.Pool.New<Int64>(outDims.NumberOfElements);
	1866	mode = BinOptItExMode.VAN;
	1867	} else {
	1868	mode = BinOptItExMode.VAI;
	1869	}
	1870	dimLen = A.Length;
	1871	} else if (B.IsVector) {
	1872	outDims = A.S;
	1873	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	1874	retArr = ILMemoryPool.Pool.New<Int64>(outDims.NumberOfElements);
	1875	mode = BinOptItExMode.AVN;
	1876	} else {
	1877	mode = BinOptItExMode.AVI;
	1878	}
	1879	dimLen = B.Length;
	1880	} else {
	1881	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	1882	}
	1883	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	1884	arrStepInc = outDims.SequentialIndexDistance(dim);
	1885	#endregion
	1886
	1887	#region worker loops definition
	1888	ILDenseStorage<Int64> retStorage = new ILDenseStorage<Int64>(retArr, outDims);
	1889	int workerCount = 1;
	1890	Action<object> worker = data => {
	1891	// expects: iStart, iLen, ap, bp, cp
	1892	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	1893	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	1894
	1895	Int64* ap;
	1896
	1897	Int64* bp;
	1898
	1899	Int64* cp;
	1900	switch (mode) {
	1901	case BinOptItExMode.VAN:
	1902	for (int s = 0; s < range.Item2; s++) {
	1903	ap = (Int64*)range.Item3;
	1904	bp = (Int64)range.Item4 + range.Item1 + s arrStepInc; ;
	1905	cp = (Int64)range.Item5 + range.Item1 + s arrStepInc;
	1906	for (int l = 0; l < dimLen; l++) {
	1907
	1908	cp = applyFunc(ap, *bp);
	1909	ap++;
	1910	bp += arrInc;
	1911	cp += arrInc;
	1912	}
	1913	}
	1914	break;
	1915	case BinOptItExMode.VAI:
	1916	for (int s = 0; s < range.Item2; s++) {
	1917	ap = (Int64*)range.Item3;
	1918	cp = (Int64)range.Item5 + range.Item1 + s arrStepInc;
	1919	for (int l = 0; l < dimLen; l++) {
	1920
	1921	cp = applyFunc(ap, *cp);
	1922	ap++;
	1923	cp += arrInc;
	1924	}
	1925	}
	1926	break;
	1927	case BinOptItExMode.AVN:
	1928	for (int s = 0; s < range.Item2; s++) {
	1929	ap = (Int64)range.Item3 + range.Item1 + s arrStepInc;
	1930	bp = (Int64*)range.Item4;
	1931	cp = (Int64)range.Item5 + range.Item1 + s arrStepInc;
	1932	for (int l = 0; l < dimLen; l++) {
	1933
	1934	cp = applyFunc(ap, *bp);
	1935	ap += arrInc;
	1936	bp++;
	1937	cp += arrInc;
	1938	}
	1939	}
	1940	break;
	1941	case BinOptItExMode.AVI:
	1942	for (int s = 0; s < range.Item2; s++) {
	1943	bp = (Int64*)range.Item4;
	1944	cp = (Int64)range.Item5 + range.Item1 + s arrStepInc;
	1945	for (int l = 0; l < dimLen; l++) {
	1946
	1947	cp = applyFunc(cp, *bp);
	1948	bp++;
	1949	cp += arrInc;
	1950	}
	1951	}
	1952	break;
	1953	}
	1954	System.Threading.Interlocked.Decrement(ref workerCount);
	1955	};
	1956	#endregion
	1957
	1958	#region work distribution
	1959	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	1960	int outLen = outDims.NumberOfElements;
	1961	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	1962	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	1963	workItemLength = outLen / dimLen / workItemCount;
	1964	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	1965	} else {
	1966	workItemLength = outLen / dimLen / 2;
	1967	workItemCount = 2;
	1968	}
	1969	} else {
	1970	workItemLength = outLen / dimLen;
	1971	workItemCount = 1;
	1972	}
	1973
	1974	fixed (Int64* arrAP = arrA)
	1975	fixed (Int64* arrBP = arrB)
	1976	fixed (Int64* retArrP = retArr) {
	1977
	1978	for (; i < workItemCount - 1; i++) {
	1979	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	1980	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1981	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	1982	System.Threading.Interlocked.Increment(ref workerCount);
	1983	ILThreadPool.QueueUserWorkItem(i, worker, range);
	1984	}
	1985	// the last (or may the only) chunk is done right here
	1986	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	1987	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	1988	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	1989
	1990	ILThreadPool.Wait4Workers(ref workerCount);
	1991	}
	1992	#endregion
	1993
	1994	return new ILRetArray<Int64>(retStorage);
	1995	}
	1996	/// <summary>Apply an arbitrary function to two arrays</summary>
	1997	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	1998	/// <param name="A">Input array A</param>
	1999	/// <param name="B">Input array B</param>
	2000	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	2001	/// <item>
	2002	/// <term>size(A) == size(B)</term>
	2003	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	2004	/// </item>
	2005	/// <item>
	2006	/// <term>isscalar(A) \|\| isscalar(B)</term>
	2007	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	2008	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	2009	/// </item>
	2010	/// <item>
	2011	/// <term>All other cases</term>
	2012	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	2013	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	2014	/// </item>
	2015	/// </list></returns>
	2016	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	2017	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	2018	/// This feature is, however, officiallny not supported.</para></remarks>
	2019	public unsafe static ILRetArray<Int32> apply(Func<Int32, Int32, Int32> func, ILInArray<Int32> A, ILInArray<Int32> B) {
	2020	using (ILScope.Enter(A, B)) {
	2021	int outLen;
	2022	BinOpItMode mode;
	2023
	2024	Int32[] retArr;
	2025
	2026	Int32[] arrA = A.GetArrayForRead();
	2027
	2028	Int32[] arrB = B.GetArrayForRead();
	2029	ILSize outDims;
	2030	#region determine operation mode
	2031	if (A.IsScalar) {
	2032	outDims = B.Size;
	2033	if (B.IsScalar) {
	2034
	2035	return new ILRetArray<Int32>(new Int32[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	2036	} else if (B.IsEmpty) {
	2037	return ILRetArray<Int32>.empty(outDims);
	2038	} else {
	2039	outLen = outDims.NumberOfElements;
	2040	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	2041	retArr = ILMemoryPool.Pool.New<Int32>(outLen);
	2042	mode = BinOpItMode.SAN;
	2043	} else {
	2044	mode = BinOpItMode.SAI;
	2045	}
	2046	}
	2047	} else {
	2048	outDims = A.Size;
	2049	if (B.IsScalar) {
	2050	if (A.IsEmpty) {
	2051	return ILRetArray<Int32>.empty(A.Size);
	2052	}
	2053	outLen = A.S.NumberOfElements;
	2054	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	2055	retArr = ILMemoryPool.Pool.New<Int32>(outLen);
	2056	mode = BinOpItMode.ASN;
	2057	} else {
	2058	mode = BinOpItMode.ASI;
	2059	}
	2060	} else {
	2061	// array + array
	2062	if (!A.Size.IsSameSize(B.Size)) {
	2063	return applyEx(func,A,B);
	2064	}
	2065	outLen = A.S.NumberOfElements;
	2066	if (A.TryGetStorage4InplaceOp(out retArr))
	2067	mode = BinOpItMode.AAIA;
	2068	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	2069	mode = BinOpItMode.AAIB;
	2070	} else {
	2071	retArr = ILMemoryPool.Pool.New<Int32>(outLen);
	2072	mode = BinOpItMode.AAN;
	2073	}
	2074	}
	2075	}
	2076	#endregion
	2077	ILDenseStorage<Int32> retStorage = new ILDenseStorage<Int32>(retArr, outDims);
	2078	int i = 0, workerCount = 1;
	2079	Action<object> worker = data => {
	2080	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	2081	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	2082
	2083	Int32* cLast, cp = (Int32*)range.Item5 + range.Item1;
	2084
	2085	Int32 scalar;
	2086	cLast = cp + range.Item2;
	2087	#region loops
	2088	switch (mode) {
	2089	case BinOpItMode.AAIA:
	2090
	2091	Int32* bp = ((Int32*)range.Item4 + range.Item1);
	2092	while (cp < cLast) {
	2093
	2094	cp = func(cp, *bp++);
	2095	cp++;
	2096	}
	2097	break;
	2098	case BinOpItMode.AAIB:
	2099
	2100	Int32* ap = ((Int32*)range.Item3 + range.Item1);
	2101	while (cp < cLast) {
	2102
	2103	cp = func(ap++, *cp);
	2104	cp++;
	2105
	2106	}
	2107	//ap = ((double*)range.Item3 + range.Item1);
	2108	//for (int i2 = range.Item2; i2-- > 0; ) {
	2109	// (cp + i2) = (ap + i2) - *(cp + i2);
	2110	//}
	2111	//int ie = range.Item1 + range.Item2-1;
	2112	//double[] locRetArr = retArr;
	2113	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	2114	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	2115	// if (i2 >= ie) break;
	2116	//}
	2117
	2118	break;
	2119	case BinOpItMode.AAN:
	2120	ap = ((Int32*)range.Item3 + range.Item1);
	2121	bp = ((Int32*)range.Item4 + range.Item1);
	2122	while (cp < cLast) {
	2123
	2124	cp++ = func(ap++, *bp++);
	2125	}
	2126	break;
	2127	case BinOpItMode.ASI:
	2128	scalar = ((Int32)range.Item4);
	2129	while (cp < cLast) {
	2130
	2131	cp = func(cp, scalar);
	2132	cp++;
	2133	}
	2134	break;
	2135	case BinOpItMode.ASN:
	2136	ap = ((Int32*)range.Item3 + range.Item1);
	2137	scalar = ((Int32)range.Item4);
	2138	while (cp < cLast) {
	2139
	2140	cp++ = func(ap++, scalar);
	2141	}
	2142	break;
	2143	case BinOpItMode.SAI:
	2144	scalar = ((Int32)range.Item3);
	2145	while (cp < cLast) {
	2146
	2147	cp = func(scalar, cp);
	2148	cp++;
	2149	}
	2150	break;
	2151	case BinOpItMode.SAN:
	2152	scalar = ((Int32)range.Item3);
	2153	bp = ((Int32*)range.Item4 + range.Item1);
	2154	while (cp < cLast) {
	2155
	2156	cp++ = func(scalar, bp++);
	2157	}
	2158	break;
	2159	default:
	2160	break;
	2161	}
	2162	#endregion
	2163	System.Threading.Interlocked.Decrement(ref workerCount);
	2164	//retStorage.PendingEvents.Signal();
	2165	};
	2166
	2167	#region do the work
	2168	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	2169	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	2170	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	2171	workItemLength = outLen / workItemCount;
	2172	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	2173	} else {
	2174	workItemLength = outLen / 2;
	2175	workItemCount = 2;
	2176	}
	2177	} else {
	2178	workItemLength = outLen;
	2179	workItemCount = 1;
	2180	}
	2181
	2182	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	2183
	2184	fixed ( Int32* arrAP = arrA)
	2185	fixed ( Int32* arrBP = arrB)
	2186	fixed ( Int32* retArrP = retArr) {
	2187
	2188	for (; i < workItemCount - 1; i++) {
	2189	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	2190	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	2191	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	2192	System.Threading.Interlocked.Increment(ref workerCount);
	2193	ILThreadPool.QueueUserWorkItem(i, worker, range);
	2194	}
	2195	// the last (or may the only) chunk is done right here
	2196	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	2197	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	2198	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	2199
	2200	System.Threading.SpinWait.SpinUntil(() => {
	2201	return workerCount <= 0;
	2202	});
	2203	//while (workerCount > 0) ;
	2204	}
	2205
	2206	#endregion
	2207	return new ILRetArray<Int32>(retStorage);
	2208	}
	2209	}
	2210
	2211	private static unsafe ILRetArray<Int32> applyEx(Func<Int32, Int32, Int32> applyFunc, ILInArray<Int32> A, ILInArray<Int32> B) {
	2212	#region parameter checking
	2213	if (isnull(A) \|\| isnull(B))
	2214	return empty<Int32>(ILSize.Empty00);
	2215	if (A.IsEmpty) {
	2216	return empty<Int32>(B.S);
	2217	} else if (B.IsEmpty) {
	2218	return empty<Int32>(A.S);
	2219	}
	2220	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	2221	// return add(A,B);
	2222	int dim = -1;
	2223	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	2224	if (A.S[_L] != B.S[_L]) {
	2225	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	2226	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	2227	}
	2228	dim = _L;
	2229	}
	2230	}
	2231	if (dim > 1)
	2232	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	2233	#endregion
	2234
	2235	#region parameter preparation
	2236
	2237
	2238	Int32[] retArr;
	2239
	2240
	2241	Int32[] arrA = A.GetArrayForRead();
	2242
	2243
	2244	Int32[] arrB = B.GetArrayForRead();
	2245	ILSize outDims;
	2246	BinOptItExMode mode;
	2247	int arrInc = 0;
	2248	int arrStepInc = 0;
	2249	int dimLen = 0;
	2250	if (A.IsVector) {
	2251	outDims = B.S;
	2252	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	2253	retArr = ILMemoryPool.Pool.New<Int32>(outDims.NumberOfElements);
	2254	mode = BinOptItExMode.VAN;
	2255	} else {
	2256	mode = BinOptItExMode.VAI;
	2257	}
	2258	dimLen = A.Length;
	2259	} else if (B.IsVector) {
	2260	outDims = A.S;
	2261	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	2262	retArr = ILMemoryPool.Pool.New<Int32>(outDims.NumberOfElements);
	2263	mode = BinOptItExMode.AVN;
	2264	} else {
	2265	mode = BinOptItExMode.AVI;
	2266	}
	2267	dimLen = B.Length;
	2268	} else {
	2269	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	2270	}
	2271	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	2272	arrStepInc = outDims.SequentialIndexDistance(dim);
	2273	#endregion
	2274
	2275	#region worker loops definition
	2276	ILDenseStorage<Int32> retStorage = new ILDenseStorage<Int32>(retArr, outDims);
	2277	int workerCount = 1;
	2278	Action<object> worker = data => {
	2279	// expects: iStart, iLen, ap, bp, cp
	2280	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	2281	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	2282
	2283	Int32* ap;
	2284
	2285	Int32* bp;
	2286
	2287	Int32* cp;
	2288	switch (mode) {
	2289	case BinOptItExMode.VAN:
	2290	for (int s = 0; s < range.Item2; s++) {
	2291	ap = (Int32*)range.Item3;
	2292	bp = (Int32)range.Item4 + range.Item1 + s arrStepInc; ;
	2293	cp = (Int32)range.Item5 + range.Item1 + s arrStepInc;
	2294	for (int l = 0; l < dimLen; l++) {
	2295
	2296	cp = applyFunc(ap, *bp);
	2297	ap++;
	2298	bp += arrInc;
	2299	cp += arrInc;
	2300	}
	2301	}
	2302	break;
	2303	case BinOptItExMode.VAI:
	2304	for (int s = 0; s < range.Item2; s++) {
	2305	ap = (Int32*)range.Item3;
	2306	cp = (Int32)range.Item5 + range.Item1 + s arrStepInc;
	2307	for (int l = 0; l < dimLen; l++) {
	2308
	2309	cp = applyFunc(ap, *cp);
	2310	ap++;
	2311	cp += arrInc;
	2312	}
	2313	}
	2314	break;
	2315	case BinOptItExMode.AVN:
	2316	for (int s = 0; s < range.Item2; s++) {
	2317	ap = (Int32)range.Item3 + range.Item1 + s arrStepInc;
	2318	bp = (Int32*)range.Item4;
	2319	cp = (Int32)range.Item5 + range.Item1 + s arrStepInc;
	2320	for (int l = 0; l < dimLen; l++) {
	2321
	2322	cp = applyFunc(ap, *bp);
	2323	ap += arrInc;
	2324	bp++;
	2325	cp += arrInc;
	2326	}
	2327	}
	2328	break;
	2329	case BinOptItExMode.AVI:
	2330	for (int s = 0; s < range.Item2; s++) {
	2331	bp = (Int32*)range.Item4;
	2332	cp = (Int32)range.Item5 + range.Item1 + s arrStepInc;
	2333	for (int l = 0; l < dimLen; l++) {
	2334
	2335	cp = applyFunc(cp, *bp);
	2336	bp++;
	2337	cp += arrInc;
	2338	}
	2339	}
	2340	break;
	2341	}
	2342	System.Threading.Interlocked.Decrement(ref workerCount);
	2343	};
	2344	#endregion
	2345
	2346	#region work distribution
	2347	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	2348	int outLen = outDims.NumberOfElements;
	2349	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	2350	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	2351	workItemLength = outLen / dimLen / workItemCount;
	2352	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	2353	} else {
	2354	workItemLength = outLen / dimLen / 2;
	2355	workItemCount = 2;
	2356	}
	2357	} else {
	2358	workItemLength = outLen / dimLen;
	2359	workItemCount = 1;
	2360	}
	2361
	2362	fixed (Int32* arrAP = arrA)
	2363	fixed (Int32* arrBP = arrB)
	2364	fixed (Int32* retArrP = retArr) {
	2365
	2366	for (; i < workItemCount - 1; i++) {
	2367	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	2368	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	2369	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	2370	System.Threading.Interlocked.Increment(ref workerCount);
	2371	ILThreadPool.QueueUserWorkItem(i, worker, range);
	2372	}
	2373	// the last (or may the only) chunk is done right here
	2374	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	2375	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	2376	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	2377
	2378	ILThreadPool.Wait4Workers(ref workerCount);
	2379	}
	2380	#endregion
	2381
	2382	return new ILRetArray<Int32>(retStorage);
	2383	}
	2384	/// <summary>Apply an arbitrary function to two arrays</summary>
	2385	/// <param name="func">A function c = f(a,b), which will be applied to elements in A and B</param>
	2386	/// <param name="A">Input array A</param>
	2387	/// <param name="B">Input array B</param>
	2388	/// <returns>The combination of A and B. The result and size depends on the inputs:<list type="table">
	2389	/// <item>
	2390	/// <term>size(A) == size(B)</term>
	2391	/// <description>Same size as A/B, elementwise combination of A and B.</description>
	2392	/// </item>
	2393	/// <item>
	2394	/// <term>isscalar(A) \|\| isscalar(B)</term>
	2395	/// <description>Same size as A or B, whichever is not a scalar, the scalar value being applied to each element
	2396	/// (i.e. if the non-scalar input is empty, the result is empty).</description>
	2397	/// </item>
	2398	/// <item>
	2399	/// <term>All other cases</term>
	2400	/// <description>If A or B is a colum vector and the other parameter is an array with a matching column length, the vector is used to operate on all columns of the array.
	2401	/// Similarly, if one parameter is a row vector, it is used to operate along the rows of the other array if its number of columns matches the vector length.</description>
	2402	/// </item>
	2403	/// </list></returns>
	2404	/// <remarks><para>The <c>apply</c> function is also implemented for input if e.g. sizes (mxn) and (mx1).
	2405	/// In this case the vector argument will be combined to each column, resulting in an (mxn) array.
	2406	/// This feature is, however, officiallny not supported.</para></remarks>
	2407	public unsafe static ILRetArray<byte> apply(Func<byte, byte, byte> func, ILInArray<byte> A, ILInArray<byte> B) {
	2408	using (ILScope.Enter(A, B)) {
	2409	int outLen;
	2410	BinOpItMode mode;
	2411
	2412	byte[] retArr;
	2413
	2414	byte[] arrA = A.GetArrayForRead();
	2415
	2416	byte[] arrB = B.GetArrayForRead();
	2417	ILSize outDims;
	2418	#region determine operation mode
	2419	if (A.IsScalar) {
	2420	outDims = B.Size;
	2421	if (B.IsScalar) {
	2422
	2423	return new ILRetArray<byte>(new byte[1] { func(A.GetValue(0), B.GetValue(0)) }, A.Size);
	2424	} else if (B.IsEmpty) {
	2425	return ILRetArray<byte>.empty(outDims);
	2426	} else {
	2427	outLen = outDims.NumberOfElements;
	2428	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	2429	retArr = ILMemoryPool.Pool.New<byte>(outLen);
	2430	mode = BinOpItMode.SAN;
	2431	} else {
	2432	mode = BinOpItMode.SAI;
	2433	}
	2434	}
	2435	} else {
	2436	outDims = A.Size;
	2437	if (B.IsScalar) {
	2438	if (A.IsEmpty) {
	2439	return ILRetArray<byte>.empty(A.Size);
	2440	}
	2441	outLen = A.S.NumberOfElements;
	2442	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	2443	retArr = ILMemoryPool.Pool.New<byte>(outLen);
	2444	mode = BinOpItMode.ASN;
	2445	} else {
	2446	mode = BinOpItMode.ASI;
	2447	}
	2448	} else {
	2449	// array + array
	2450	if (!A.Size.IsSameSize(B.Size)) {
	2451	return applyEx(func,A,B);
	2452	}
	2453	outLen = A.S.NumberOfElements;
	2454	if (A.TryGetStorage4InplaceOp(out retArr))
	2455	mode = BinOpItMode.AAIA;
	2456	else if (B.TryGetStorage4InplaceOp(out retArr)) {
	2457	mode = BinOpItMode.AAIB;
	2458	} else {
	2459	retArr = ILMemoryPool.Pool.New<byte>(outLen);
	2460	mode = BinOpItMode.AAN;
	2461	}
	2462	}
	2463	}
	2464	#endregion
	2465	ILDenseStorage<byte> retStorage = new ILDenseStorage<byte>(retArr, outDims);
	2466	int i = 0, workerCount = 1;
	2467	Action<object> worker = data => {
	2468	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	2469	= (Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>)data;
	2470
	2471	byte* cLast, cp = (byte*)range.Item5 + range.Item1;
	2472
	2473	byte scalar;
	2474	cLast = cp + range.Item2;
	2475	#region loops
	2476	switch (mode) {
	2477	case BinOpItMode.AAIA:
	2478
	2479	byte* bp = ((byte*)range.Item4 + range.Item1);
	2480	while (cp < cLast) {
	2481
	2482	cp = func(cp, *bp++);
	2483	cp++;
	2484	}
	2485	break;
	2486	case BinOpItMode.AAIB:
	2487
	2488	byte* ap = ((byte*)range.Item3 + range.Item1);
	2489	while (cp < cLast) {
	2490
	2491	cp = func(ap++, *cp);
	2492	cp++;
	2493
	2494	}
	2495	//ap = ((double*)range.Item3 + range.Item1);
	2496	//for (int i2 = range.Item2; i2-- > 0; ) {
	2497	// (cp + i2) = (ap + i2) - *(cp + i2);
	2498	//}
	2499	//int ie = range.Item1 + range.Item2-1;
	2500	//double[] locRetArr = retArr;
	2501	//for (int i2 = range.Item1; i2 < locRetArr.Length; i2++) {
	2502	// locRetArr[i2] = arrA[i2] - locRetArr[i2];
	2503	// if (i2 >= ie) break;
	2504	//}
	2505
	2506	break;
	2507	case BinOpItMode.AAN:
	2508	ap = ((byte*)range.Item3 + range.Item1);
	2509	bp = ((byte*)range.Item4 + range.Item1);
	2510	while (cp < cLast) {
	2511
	2512	cp++ = func(ap++, *bp++);
	2513	}
	2514	break;
	2515	case BinOpItMode.ASI:
	2516	scalar = ((byte)range.Item4);
	2517	while (cp < cLast) {
	2518
	2519	cp = func(cp, scalar);
	2520	cp++;
	2521	}
	2522	break;
	2523	case BinOpItMode.ASN:
	2524	ap = ((byte*)range.Item3 + range.Item1);
	2525	scalar = ((byte)range.Item4);
	2526	while (cp < cLast) {
	2527
	2528	cp++ = func(ap++, scalar);
	2529	}
	2530	break;
	2531	case BinOpItMode.SAI:
	2532	scalar = ((byte)range.Item3);
	2533	while (cp < cLast) {
	2534
	2535	cp = func(scalar, cp);
	2536	cp++;
	2537	}
	2538	break;
	2539	case BinOpItMode.SAN:
	2540	scalar = ((byte)range.Item3);
	2541	bp = ((byte*)range.Item4 + range.Item1);
	2542	while (cp < cLast) {
	2543
	2544	cp++ = func(scalar, bp++);
	2545	}
	2546	break;
	2547	default:
	2548	break;
	2549	}
	2550	#endregion
	2551	System.Threading.Interlocked.Decrement(ref workerCount);
	2552	//retStorage.PendingEvents.Signal();
	2553	};
	2554
	2555	#region do the work
	2556	int workItemCount = Settings.s_maxNumberThreads, workItemLength;
	2557	if (Settings.s_maxNumberThreads > 1 && outLen / 2 > Settings.s_minParallelElement1Count) {
	2558	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	2559	workItemLength = outLen / workItemCount;
	2560	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	2561	} else {
	2562	workItemLength = outLen / 2;
	2563	workItemCount = 2;
	2564	}
	2565	} else {
	2566	workItemLength = outLen;
	2567	workItemCount = 1;
	2568	}
	2569
	2570	// retStorage.PendingEvents = new System.Threading.CountdownEvent(workItemCount);
	2571
	2572	fixed ( byte* arrAP = arrA)
	2573	fixed ( byte* arrBP = arrB)
	2574	fixed ( byte* retArrP = retArr) {
	2575
	2576	for (; i < workItemCount - 1; i++) {
	2577	Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode> range
	2578	= new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	2579	(i * workItemLength, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode);
	2580	System.Threading.Interlocked.Increment(ref workerCount);
	2581	ILThreadPool.QueueUserWorkItem(i, worker, range);
	2582	}
	2583	// the last (or may the only) chunk is done right here
	2584	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	2585	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr, BinOpItMode>
	2586	(i * workItemLength, outLen - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP, mode));
	2587
	2588	System.Threading.SpinWait.SpinUntil(() => {
	2589	return workerCount <= 0;
	2590	});
	2591	//while (workerCount > 0) ;
	2592	}
	2593
	2594	#endregion
	2595	return new ILRetArray<byte>(retStorage);
	2596	}
	2597	}
	2598
	2599	private static unsafe ILRetArray<byte> applyEx(Func<byte, byte, byte> applyFunc, ILInArray<byte> A, ILInArray<byte> B) {
	2600	#region parameter checking
	2601	if (isnull(A) \|\| isnull(B))
	2602	return empty<byte>(ILSize.Empty00);
	2603	if (A.IsEmpty) {
	2604	return empty<byte>(B.S);
	2605	} else if (B.IsEmpty) {
	2606	return empty<byte>(A.S);
	2607	}
	2608	//if (A.IsScalar \|\| B.IsScalar \|\| A.D.IsSameSize(B.D))
	2609	// return add(A,B);
	2610	int dim = -1;
	2611	for (int _L = 0; _L < Math.Max(A.S.NumberOfDimensions, B.S.NumberOfDimensions); _L++) {
	2612	if (A.S[_L] != B.S[_L]) {
	2613	if (dim >= 0 \|\| (A.S[_L] != 1 && B.S[_L] != 1)) {
	2614	throw new ILArgumentException("A and B must have the same size except for one singleton dimension in A or B");
	2615	}
	2616	dim = _L;
	2617	}
	2618	}
	2619	if (dim > 1)
	2620	throw new ILArgumentException("singleton dimension expansion currently is only supported for colum- and row vectors");
	2621	#endregion
	2622
	2623	#region parameter preparation
	2624
	2625
	2626	byte[] retArr;
	2627
	2628
	2629	byte[] arrA = A.GetArrayForRead();
	2630
	2631
	2632	byte[] arrB = B.GetArrayForRead();
	2633	ILSize outDims;
	2634	BinOptItExMode mode;
	2635	int arrInc = 0;
	2636	int arrStepInc = 0;
	2637	int dimLen = 0;
	2638	if (A.IsVector) {
	2639	outDims = B.S;
	2640	if (!B.TryGetStorage4InplaceOp(out retArr)) {
	2641	retArr = ILMemoryPool.Pool.New<byte>(outDims.NumberOfElements);
	2642	mode = BinOptItExMode.VAN;
	2643	} else {
	2644	mode = BinOptItExMode.VAI;
	2645	}
	2646	dimLen = A.Length;
	2647	} else if (B.IsVector) {
	2648	outDims = A.S;
	2649	if (!A.TryGetStorage4InplaceOp(out retArr)) {
	2650	retArr = ILMemoryPool.Pool.New<byte>(outDims.NumberOfElements);
	2651	mode = BinOptItExMode.AVN;
	2652	} else {
	2653	mode = BinOptItExMode.AVI;
	2654	}
	2655	dimLen = B.Length;
	2656	} else {
	2657	throw new ILArgumentException("A and B must have the same size except for one simgleton dimension in A or B");
	2658	}
	2659	arrInc = (dim == 0) ? outDims.SequentialIndexDistance(1) : outDims.SequentialIndexDistance(0);
	2660	arrStepInc = outDims.SequentialIndexDistance(dim);
	2661	#endregion
	2662
	2663	#region worker loops definition
	2664	ILDenseStorage<byte> retStorage = new ILDenseStorage<byte>(retArr, outDims);
	2665	int workerCount = 1;
	2666	Action<object> worker = data => {
	2667	// expects: iStart, iLen, ap, bp, cp
	2668	Tuple<int, int, IntPtr, IntPtr, IntPtr> range =
	2669	(Tuple<int, int, IntPtr, IntPtr, IntPtr>)data;
	2670
	2671	byte* ap;
	2672
	2673	byte* bp;
	2674
	2675	byte* cp;
	2676	switch (mode) {
	2677	case BinOptItExMode.VAN:
	2678	for (int s = 0; s < range.Item2; s++) {
	2679	ap = (byte*)range.Item3;
	2680	bp = (byte)range.Item4 + range.Item1 + s arrStepInc; ;
	2681	cp = (byte)range.Item5 + range.Item1 + s arrStepInc;
	2682	for (int l = 0; l < dimLen; l++) {
	2683
	2684	cp = applyFunc(ap, *bp);
	2685	ap++;
	2686	bp += arrInc;
	2687	cp += arrInc;
	2688	}
	2689	}
	2690	break;
	2691	case BinOptItExMode.VAI:
	2692	for (int s = 0; s < range.Item2; s++) {
	2693	ap = (byte*)range.Item3;
	2694	cp = (byte)range.Item5 + range.Item1 + s arrStepInc;
	2695	for (int l = 0; l < dimLen; l++) {
	2696
	2697	cp = applyFunc(ap, *cp);
	2698	ap++;
	2699	cp += arrInc;
	2700	}
	2701	}
	2702	break;
	2703	case BinOptItExMode.AVN:
	2704	for (int s = 0; s < range.Item2; s++) {
	2705	ap = (byte)range.Item3 + range.Item1 + s arrStepInc;
	2706	bp = (byte*)range.Item4;
	2707	cp = (byte)range.Item5 + range.Item1 + s arrStepInc;
	2708	for (int l = 0; l < dimLen; l++) {
	2709
	2710	cp = applyFunc(ap, *bp);
	2711	ap += arrInc;
	2712	bp++;
	2713	cp += arrInc;
	2714	}
	2715	}
	2716	break;
	2717	case BinOptItExMode.AVI:
	2718	for (int s = 0; s < range.Item2; s++) {
	2719	bp = (byte*)range.Item4;
	2720	cp = (byte)range.Item5 + range.Item1 + s arrStepInc;
	2721	for (int l = 0; l < dimLen; l++) {
	2722
	2723	cp = applyFunc(cp, *bp);
	2724	bp++;
	2725	cp += arrInc;
	2726	}
	2727	}
	2728	break;
	2729	}
	2730	System.Threading.Interlocked.Decrement(ref workerCount);
	2731	};
	2732	#endregion
	2733
	2734	#region work distribution
	2735	int i = 0, workItemCount = Settings.s_maxNumberThreads, workItemLength;
	2736	int outLen = outDims.NumberOfElements;
	2737	if (Settings.s_maxNumberThreads > 1 && outLen / 2 >= Settings.s_minParallelElement1Count) {
	2738	if (outLen / workItemCount > Settings.s_minParallelElement1Count) {
	2739	workItemLength = outLen / dimLen / workItemCount;
	2740	//workItemLength = (int)((double)outLen / workItemCount * 1.05);
	2741	} else {
	2742	workItemLength = outLen / dimLen / 2;
	2743	workItemCount = 2;
	2744	}
	2745	} else {
	2746	workItemLength = outLen / dimLen;
	2747	workItemCount = 1;
	2748	}
	2749
	2750	fixed (byte* arrAP = arrA)
	2751	fixed (byte* arrBP = arrB)
	2752	fixed (byte* retArrP = retArr) {
	2753
	2754	for (; i < workItemCount - 1; i++) {
	2755	Tuple<int, int, IntPtr, IntPtr, IntPtr> range
	2756	= new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	2757	(i * workItemLength * arrStepInc, workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP);
	2758	System.Threading.Interlocked.Increment(ref workerCount);
	2759	ILThreadPool.QueueUserWorkItem(i, worker, range);
	2760	}
	2761	// the last (or may the only) chunk is done right here
	2762	//System.Threading.Interlocked.Increment(ref retStorage.PendingTasks);
	2763	worker(new Tuple<int, int, IntPtr, IntPtr, IntPtr>
	2764	(i * workItemLength * arrStepInc, (outLen / dimLen) - i * workItemLength, (IntPtr)arrAP, (IntPtr)arrBP, (IntPtr)retArrP));
	2765
	2766	ILThreadPool.Wait4Workers(ref workerCount);
	2767	}
	2768	#endregion
	2769
	2770	return new ILRetArray<byte>(retStorage);
	2771	}
	2772
	2773	#endregion HYCALPER AUTO GENERATED CODE
	2774
	2775	}
	2776	}

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

Download in other formats:

Update cookies preferences