Context Navigation

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

Visit:

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

Line
1	///
2	/// This file is part of ILNumerics Community Edition.
3	///
4	/// ILNumerics Community Edition - high performance computing for applications.
5	/// Copyright (C) 2006 - 2012 Haymo Kutschbach, http://ilnumerics.net
6	///
7	/// ILNumerics Community Edition is free software: you can redistribute it and/or modify
8	/// it under the terms of the GNU General Public License version 3 as published by
9	/// the Free Software Foundation.
10	///
11	/// ILNumerics Community Edition is distributed in the hope that it will be useful,
12	/// but WITHOUT ANY WARRANTY; without even the implied warranty of
13	/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14	/// GNU General Public License for more details.
15	///
16	/// You should have received a copy of the GNU General Public License
17	/// along with ILNumerics Community Edition. See the file License.txt in the root
18	/// of your distribution package. If not, see <http://www.gnu.org/licenses/>.
19	///
20	/// In addition this software uses the following components and/or licenses:
21	///
22	/// =================================================================================
23	/// The Open Toolkit Library License
24	///
25	/// Copyright (c) 2006 - 2009 the Open Toolkit library.
26	///
27	/// Permission is hereby granted, free of charge, to any person obtaining a copy
28	/// of this software and associated documentation files (the "Software"), to deal
29	/// in the Software without restriction, including without limitation the rights to
30	/// use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
31	/// the Software, and to permit persons to whom the Software is furnished to do
32	/// so, subject to the following conditions:
33	///
34	/// The above copyright notice and this permission notice shall be included in all
35	/// copies or substantial portions of the Software.
36	///
37	/// =================================================================================
38	///
39
40	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