Mercurial > hg > camir-aes2014

comparison toolboxes/FullBNT-1.0.7/bnt/potentials/genops.c @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e9a9cd732c1e
1 /*
2
3 GENOPS.C
4 Generalized arithmetic operators overloading built-in functions.
5
6 written by Douglas M. Schwarz
7 schwarz@servtech.com
8 26 December 1998
9 Last modified: 2 April 1999
10
11 Copyright 1998-1999 by Douglas M. Schwarz. All rights reserved.
12
13 */
14
15
16 /*
17
18 Build MEX file by entering the appropriate command at the MATLAB prompt
19 (-D<name> option is equivalent to #define <name> in source file):
20
21 mex genops.c -DPLUS_MEX -output plus
22 mex genops.c -DMINUS_MEX -output minus
23 mex genops.c -DTIMES_MEX -output times
24 mex genops.c -DRDIVIDE_MEX -output rdivide
25 mex genops.c -DLDIVIDE_MEX -output ldivide
26 mex genops.c -DPOWER_MEX -output power
27 mex genops.c -DEQ_MEX -output eq
28 mex genops.c -DNE_MEX -output ne
29 mex genops.c -DLT_MEX -output lt
30 mex genops.c -DGT_MEX -output gt
31 mex genops.c -DLE_MEX -output le
32 mex genops.c -DGE_MEX -output ge
33
34 */
35
36 /* This file has been formatted for a tab equal to 4 spaces. */
37
38 #if defined(EQ_MEX) || defined(NE_MEX) || defined(LT_MEX) || defined(GT_MEX) \
39 || defined(LE_MEX) || defined(GE_MEX)
40 #define RELOP_MEX
41 #endif
42
43 #include "mex.h"
44 #include "matrix.h"
45 #ifdef POWER_MEX
46 #include <math.h>
47 #define PI 3.141592653589793
48 #endif
49
50 bool allequal(int, const int *, const int *);
51 void removeZeroImag(double *, double *, int, const int *, int, mxArray **);
52
53 #define xMat prhs[0]
54 #define yMat prhs[1]
55 #define zMat plhs[0]
56
57 #define min(A,B) ((A) < (B) ? (A) : (B))
58 #define max(A,B) ((A) > (B) ? (A) : (B))
59
60 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
61 {
62 double *xrp, *xip, *yrp, *yip;
63 #ifndef RELOP_MEX
64 double *zr, *zi, *zip;
65 #endif
66 double *zrp, *zrend;
67 int xnd, ynd, numElements = 1;
68 const int *xdim, *ydim;
69 bool xcmplx, ycmplx;
70 mxClassID yclass;
71 int *s, ndim, *sx, *sy, i, *cpsx, *cpsy;
72 int *subs, *s1, *cpsx2, *cpsy2;
73 int ix = 0, iy = 0;
74 mxArray *args[3], *result[1];
75 #if defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
76 double denom;
77 #endif
78 #ifdef POWER_MEX
79 double mag, theta, phi, magx;
80 int flops = 0;
81 #endif
82
83
84 if (nrhs != 2)
85 mexErrMsgTxt("Incorrect number of inputs.");
86
87 if (nlhs > 1)
88 mexErrMsgTxt("Too many output arguments.");
89
90 xnd = mxGetNumberOfDimensions(xMat);
91 ynd = mxGetNumberOfDimensions(yMat);
92 xdim = mxGetDimensions(xMat);
93 ydim = mxGetDimensions(yMat);
94
95 yclass = mxGetClassID(yMat);
96
97 /* If the built-in function in MATLAB can handle the arguments
98 then use that. */
99 if (yclass != mxDOUBLE_CLASS ||
100 (xnd == 2 && xdim[0] == 1 && xdim[1] == 1) ||
101 (ynd == 2 && ydim[0] == 1 && ydim[1] == 1) ||
102 (xnd == ynd && allequal(xnd,xdim,ydim)))
103 {
104 #ifdef PLUS_MEX
105 args[0] = mxCreateString("plus");
106 #elif defined(MINUS_MEX)
107 args[0] = mxCreateString("minus");
108 #elif defined(TIMES_MEX)
109 args[0] = mxCreateString("times");
110 #elif defined(RDIVIDE_MEX)
111 args[0] = mxCreateString("rdivide");
112 #elif defined(LDIVIDE_MEX)
113 args[0] = mxCreateString("ldivide");
114 #elif defined(POWER_MEX)
115 args[0] = mxCreateString("power");
116 #elif defined(EQ_MEX)
117 args[0] = mxCreateString("eq");
118 #elif defined(NE_MEX)
119 args[0] = mxCreateString("ne");
120 #elif defined(LT_MEX)
121 args[0] = mxCreateString("lt");
122 #elif defined(GT_MEX)
123 args[0] = mxCreateString("gt");
124 #elif defined(LE_MEX)
125 args[0] = mxCreateString("le");
126 #elif defined(GE_MEX)
127 args[0] = mxCreateString("ge");
128 #endif
129 args[1] = (mxArray *)xMat;
130 args[2] = (mxArray *)yMat;
131 mexCallMATLAB(1, result, 3, args, "builtin");
132 mxDestroyArray(args[0]);
133 zMat = result[0];
134 }
135 else /* X and Y are both N-D and different dimensionality. */
136 {
137 ndim = max(xnd,ynd);
138 sx = (int *)mxMalloc(sizeof(int)*ndim);
139 sy = (int *)mxMalloc(sizeof(int)*ndim);
140 s = (int *)mxMalloc(sizeof(int)*ndim);
141 s1 = (int *)mxMalloc(sizeof(int)*ndim);
142 *(cpsx = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
143 *(cpsy = (int *)mxMalloc(sizeof(int)*ndim)) = 1;
144 subs = (int *)mxMalloc(sizeof(int)*ndim);
145 cpsx2 = (int *)mxMalloc(sizeof(int)*ndim);
146 cpsy2 = (int *)mxMalloc(sizeof(int)*ndim);
147 for (i = 0; i < ndim; i++)
148 {
149 subs[i] = 0;
150 sx[i] = (i < xnd) ? xdim[i] : 1;
151 sy[i] = (i < ynd) ? ydim[i] : 1;
152 if (sx[i] == sy[i])
153 s[i] = sx[i];
154 else if (sx[i] == 1)
155 s[i] = sy[i];
156 else if (sy[i] == 1)
157 s[i] = sx[i];
158 else
159 {
160 mxFree(sx);
161 mxFree(sy);
162 mxFree(s);
163 mxFree(s1);
164 mxFree(cpsx);
165 mxFree(cpsy);
166 mxFree(subs);
167 mxFree(cpsx2);
168 mxFree(cpsy2);
169 mexErrMsgTxt("Array dimensions are not appropriate.");
170 }
171 s1[i] = s[i] - 1;
172 numElements *= s[i];
173 }
174
175 for (i = 0; i < ndim-1; i++)
176 {
177 cpsx[i+1] = cpsx[i]*sx[i]--;
178 cpsy[i+1] = cpsy[i]*sy[i]--;
179 cpsx2[i] = cpsx[i]*sx[i];
180 cpsy2[i] = cpsy[i]*sy[i];
181 }
182 cpsx2[ndim-1] = cpsx[ndim-1]*(--sx[ndim-1]);
183 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
184
185 xcmplx = mxIsComplex(xMat);
186 ycmplx = mxIsComplex(yMat);
187
188 if (!xcmplx && !ycmplx) /* X and Y both N-D, both real. */
189 {
190 #ifdef POWER_MEX
191 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
192 zrp = zr = mxGetPr(zMat);
193 zip = zi = mxGetPi(zMat);
194 #elif defined(RELOP_MEX)
195 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
196 mxSetLogical(zMat);
197 zrp = mxGetPr(zMat);
198 #else
199 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
200 zrp = mxGetPr(zMat);
201 #endif
202 xrp = mxGetPr(xMat);
203 yrp = mxGetPr(yMat);
204 zrend = zrp + numElements;
205 while (zrp < zrend)
206 {
207 #ifdef PLUS_MEX
208 *zrp++ = *xrp + *yrp;
209 #elif defined(MINUS_MEX)
210 *zrp++ = *xrp - *yrp;
211 #elif defined(TIMES_MEX)
212 *zrp++ = *xrp * *yrp;
213 #elif defined(RDIVIDE_MEX)
214 *zrp++ = *xrp / *yrp;
215 #elif defined(LDIVIDE_MEX)
216 *zrp++ = *yrp / *xrp;
217 #elif defined(POWER_MEX)
218 if (*xrp < 0.0 && *yrp != floor(*yrp))
219 {
220 mag = pow(-*xrp,*yrp);
221 theta = PI * *yrp;
222 *zrp++ = mag*cos(theta);
223 *zip++ = mag*sin(theta);
224 flops += 18;
225 }
226 else
227 {
228 *zrp++ = pow(*xrp,*yrp);
229 *zip++ = 0.0;
230 flops++;
231 }
232 #elif defined(EQ_MEX)
233 *zrp++ = (*xrp == *yrp);
234 #elif defined(NE_MEX)
235 *zrp++ = (*xrp != *yrp);
236 #elif defined(LT_MEX)
237 *zrp++ = (*xrp < *yrp);
238 #elif defined(GT_MEX)
239 *zrp++ = (*xrp > *yrp);
240 #elif defined(LE_MEX)
241 *zrp++ = (*xrp <= *yrp);
242 #elif defined(GE_MEX)
243 *zrp++ = (*xrp >= *yrp);
244 #endif
245 for (i = 0; i < ndim; i++)
246 {
247 if (subs[i] == s1[i])
248 {
249 subs[i] = 0;
250 if (sx[i])
251 xrp -= cpsx2[i];
252 if (sy[i])
253 yrp -= cpsy2[i];
254 }
255 else
256 {
257 subs[i]++;
258 if (sx[i])
259 xrp += cpsx[i];
260 if (sy[i])
261 yrp += cpsy[i];
262 break;
263 }
264 }
265 }
266 #ifdef POWER_MEX
267 mexAddFlops(flops);
268 removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
269 #elif !defined(RELOP_MEX)
270 mexAddFlops(numElements);
271 #endif
272 }
273 else if (!xcmplx && ycmplx) /* X and Y both N-D, X real, Y complex. */
274 {
275 #ifdef POWER_MEX
276 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
277 zrp = zr = mxGetPr(zMat);
278 zip = zi = mxGetPi(zMat);
279 #elif defined(RELOP_MEX)
280 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
281 mxSetLogical(zMat);
282 zrp = mxGetPr(zMat);
283 #else
284 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
285 zrp = mxGetPr(zMat);
286 zip = mxGetPi(zMat);
287 #endif
288 xrp = mxGetPr(xMat);
289 yrp = mxGetPr(yMat);
290 yip = mxGetPi(yMat);
291 zrend = zrp + numElements;
292 while (zrp < zrend)
293 {
294 #ifdef PLUS_MEX
295 *zrp++ = *xrp + *yrp;
296 *zip++ = *yip;
297 #elif defined(MINUS_MEX)
298 *zrp++ = *xrp - *yrp;
299 *zip++ = -*yip;
300 #elif defined(TIMES_MEX)
301 *zrp++ = *xrp * *yrp;
302 *zip++ = *xrp * *yip;
303 #elif defined(RDIVIDE_MEX)
304 denom = *yrp * *yrp + *yip * *yip;
305 *zrp++ = (*xrp * *yrp)/denom;
306 *zip++ = (-*xrp * *yip)/denom;
307 #elif defined(LDIVIDE_MEX)
308 *zrp++ = *yrp / *xrp;
309 *zip++ = *yip / *xrp;
310 #elif defined(POWER_MEX)
311 if (*yip == 0.0)
312 {
313 if (*xrp < 0.0 && *yrp != floor(*yrp))
314 {
315 mag = pow(-*xrp,*yrp);
316 theta = PI * *yrp;
317 *zrp++ = mag*cos(theta);
318 *zip++ = mag*sin(theta);
319 flops += 18;
320 }
321 else
322 {
323 *zrp++ = pow(*xrp,*yrp);
324 *zip++ = 0.0;
325 flops++;
326 }
327 }
328 else
329 {
330 if (*xrp < 0.0)
331 {
332 mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
333 theta = *yip * log(-*xrp) + PI * *yrp;
334 *zrp++ = mag*cos(theta);
335 *zip++ = mag*sin(theta);
336 flops += 18;
337 }
338 else
339 {
340 mag = pow(*xrp,*yrp);
341 theta = *yip * log(*xrp);
342 *zrp++ = mag*cos(theta);
343 *zip++ = mag*sin(theta);
344 flops += 13;
345 }
346 }
347 #elif defined(EQ_MEX)
348 *zrp++ = (*xrp == *yrp) && (*yip == 0.0);
349 #elif defined(NE_MEX)
350 *zrp++ = (*xrp != *yrp) || (*yip != 0.0);
351 #elif defined(LT_MEX)
352 *zrp++ = (*xrp < *yrp);
353 #elif defined(GT_MEX)
354 *zrp++ = (*xrp > *yrp);
355 #elif defined(LE_MEX)
356 *zrp++ = (*xrp <= *yrp);
357 #elif defined(GE_MEX)
358 *zrp++ = (*xrp >= *yrp);
359 #endif
360 for (i = 0; i < ndim; i++)
361 {
362 if (subs[i] == s1[i])
363 {
364 subs[i] = 0;
365 if (sx[i])
366 xrp -= cpsx2[i];
367 if (sy[i])
368 {
369 yrp -= cpsy2[i];
370 yip -= cpsy2[i];
371 }
372 }
373 else
374 {
375 subs[i]++;
376 if (sx[i])
377 xrp += cpsx[i];
378 if (sy[i])
379 {
380 yrp += cpsy[i];
381 yip += cpsy[i];
382 }
383 break;
384 }
385 }
386 }
387 #if defined(PLUS_MEX) || defined(MINUS_MEX)
388 mexAddFlops(2*numElements);
389 #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
390 mexAddFlops(6*numElements);
391 #elif defined(POWER_MEX)
392 mexAddFlops(flops);
393 removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
394 #endif
395 }
396 else if (xcmplx && !ycmplx) /* X and Y both N-D, X complex, Y real. */
397 {
398 #ifdef POWER_MEX
399 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
400 zrp = zr = mxGetPr(zMat);
401 zip = zi = mxGetPi(zMat);
402 #elif defined(RELOP_MEX)
403 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
404 mxSetLogical(zMat);
405 zrp = mxGetPr(zMat);
406 #else
407 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
408 zrp = mxGetPr(zMat);
409 zip = mxGetPi(zMat);
410 #endif
411 xrp = mxGetPr(xMat);
412 xip = mxGetPi(xMat);
413 yrp = mxGetPr(yMat);
414 zrend = zrp + numElements;
415 while (zrp < zrend)
416 {
417 #ifdef PLUS_MEX
418 *zrp++ = *xrp + *yrp;
419 *zip++ = *xip;
420 #elif defined(MINUS_MEX)
421 *zrp++ = *xrp - *yrp;
422 *zip++ = *xip;
423 #elif defined(TIMES_MEX)
424 *zrp++ = *xrp * *yrp;
425 *zip++ = *xip * *yrp;
426 #elif defined(RDIVIDE_MEX)
427 *zrp++ = *xrp / *yrp;
428 *zip++ = *xip / *yrp;
429 #elif defined(LDIVIDE_MEX)
430 denom = *xrp * *xrp + *xip * *xip;
431 *zrp++ = (*xrp * *yrp)/denom;
432 *zip++ = (-*xip * *yrp)/denom;
433 #elif defined(POWER_MEX)
434 if (*xip == 0.0)
435 {
436 if (*xrp < 0.0 && *yrp != floor(*yrp))
437 {
438 mag = pow(-*xrp,*yrp);
439 theta = PI * *yrp;
440 *zrp++ = mag*cos(theta);
441 *zip++ = mag*sin(theta);
442 flops += 18;
443 }
444 else
445 {
446 *zrp++ = pow(*xrp,*yrp);
447 *zip++ = 0.0;
448 flops++;
449 }
450 }
451 else
452 {
453 mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
454 theta = *yrp*atan2(*xip,*xrp);
455 *zrp++ = mag*cos(theta);
456 *zip++ = mag*sin(theta);
457 flops += 18;
458 }
459 #elif defined(EQ_MEX)
460 *zrp++ = (*xrp == *yrp) && (*xip == 0.0);
461 #elif defined(NE_MEX)
462 *zrp++ = (*xrp != *yrp) || (*xip != 0.0);
463 #elif defined(LT_MEX)
464 *zrp++ = (*xrp < *yrp);
465 #elif defined(GT_MEX)
466 *zrp++ = (*xrp > *yrp);
467 #elif defined(LE_MEX)
468 *zrp++ = (*xrp <= *yrp);
469 #elif defined(GE_MEX)
470 *zrp++ = (*xrp >= *yrp);
471 #endif
472 for (i = 0; i < ndim; i++)
473 {
474 if (subs[i] == s1[i])
475 {
476 subs[i] = 0;
477 if (sx[i])
478 {
479 xrp -= cpsx2[i];
480 xip -= cpsx2[i];
481 }
482 if (sy[i])
483 yrp -= cpsy2[i];
484 }
485 else
486 {
487 subs[i]++;
488 if (sx[i])
489 {
490 xrp += cpsx[i];
491 xip += cpsx[i];
492 }
493 if (sy[i])
494 yrp += cpsy[i];
495 break;
496 }
497 }
498 }
499 #if defined(PLUS_MEX) || defined(MINUS_MEX)
500 mexAddFlops(2*numElements);
501 #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
502 mexAddFlops(6*numElements);
503 #elif defined(POWER_MEX)
504 mexAddFlops(flops);
505 removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
506 #endif
507 }
508 else if (xcmplx && ycmplx) /* X and Y both N-D, both complex. */
509 {
510 #if defined(RELOP_MEX)
511 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
512 mxSetLogical(zMat);
513 zrp = mxGetPr(zMat);
514 #else
515 zMat = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxCOMPLEX);
516 zrp = zr = mxGetPr(zMat);
517 zip = zi = mxGetPi(zMat);
518 #endif
519 xrp = mxGetPr(xMat);
520 xip = mxGetPi(xMat);
521 yrp = mxGetPr(yMat);
522 yip = mxGetPi(yMat);
523 zrend = zrp + numElements;
524 while (zrp < zrend)
525 {
526 #ifdef PLUS_MEX
527 *zrp++ = *xrp + *yrp;
528 *zip++ = *xip + *yip;
529 #elif defined(MINUS_MEX)
530 *zrp++ = *xrp - *yrp;
531 *zip++ = *xip - *yip;
532 #elif defined(TIMES_MEX)
533 *zrp++ = *xrp * *yrp - *xip * *yip;
534 *zip++ = *xip * *yrp + *xrp * *yip;
535 #elif defined(RDIVIDE_MEX)
536 denom = *yrp * *yrp + *yip * *yip;
537 *zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
538 *zip++ = (*xip * *yrp - *xrp * *yip)/denom;
539 #elif defined(LDIVIDE_MEX)
540 denom = *xrp * *xrp + *xip * *xip;
541 *zrp++ = (*xrp * *yrp + *xip * *yip)/denom;
542 *zip++ = (*xrp * *yip - *xip * *yrp)/denom;
543 #elif defined(POWER_MEX)
544 if (*xip == 0.0 && *yip == 0.0)
545 {
546 if (*xrp < 0.0 && *yrp != floor(*yrp))
547 {
548 mag = pow(-*xrp,*yrp);
549 theta = PI * *yrp;
550 *zrp++ = mag*cos(theta);
551 *zip++ = mag*sin(theta);
552 flops += 18;
553 }
554 else
555 {
556 *zrp++ = pow(*xrp,*yrp);
557 *zip++ = 0.0;
558 flops++;
559 }
560 }
561 else if (*xip == 0.0)
562 {
563 if (*xrp < 0.0)
564 {
565 mag = pow(-*xrp,*yrp)*exp(-PI * *yip);
566 theta = *yip * log(-*xrp) + PI * *yrp;
567 *zrp++ = mag*cos(theta);
568 *zip++ = mag*sin(theta);
569 flops += 18;
570 }
571 else
572 {
573 mag = pow(*xrp,*yrp);
574 theta = *yip * log(*xrp);
575 *zrp++ = mag*cos(theta);
576 *zip++ = mag*sin(theta);
577 flops += 13;
578 }
579 }
580 else if (*yip == 0.0)
581 {
582 mag = pow(*xrp * *xrp + *xip * *xip,0.5 * *yrp);
583 theta = *yrp * atan2(*xip,*xrp);
584 *zrp++ = mag*cos(theta);
585 *zip++ = mag*sin(theta);
586 flops += 18;
587 }
588 else
589 {
590 magx = sqrt(*xrp * *xrp + *xip * *xip);
591 phi = atan2(*xip,*xrp);
592 mag = pow(magx,*yrp)*exp(-*yip * phi);
593 theta = *yip * log(magx) + *yrp * phi;
594 *zrp++ = mag*cos(theta);
595 *zip++ = mag*sin(theta);
596 flops += 18;
597 }
598 #elif defined(EQ_MEX)
599 *zrp++ = (*xrp == *yrp) && (*xip == *yip);
600 #elif defined(NE_MEX)
601 *zrp++ = (*xrp != *yrp) || (*xip != *yip);
602 #elif defined(LT_MEX)
603 *zrp++ = (*xrp < *yrp);
604 #elif defined(GT_MEX)
605 *zrp++ = (*xrp > *yrp);
606 #elif defined(LE_MEX)
607 *zrp++ = (*xrp <= *yrp);
608 #elif defined(GE_MEX)
609 *zrp++ = (*xrp >= *yrp);
610 #endif
611 for (i = 0; i < ndim; i++)
612 {
613 if (subs[i] == s1[i])
614 {
615 subs[i] = 0;
616 if (sx[i])
617 {
618 xrp -= cpsx2[i];
619 xip -= cpsx2[i];
620 }
621 if (sy[i])
622 {
623 yrp -= cpsy2[i];
624 yip -= cpsy2[i];
625 }
626 }
627 else
628 {
629 subs[i]++;
630 if (sx[i])
631 {
632 xrp += cpsx[i];
633 xip += cpsx[i];
634 }
635 if (sy[i])
636 {
637 yrp += cpsy[i];
638 yip += cpsy[i];
639 }
640 break;
641 }
642 }
643 }
644 #if defined(PLUS_MEX) || defined(MINUS_MEX)
645 mexAddFlops(2*numElements);
646 #elif defined(TIMES_MEX) || defined(RDIVIDE_MEX) || defined(LDIVIDE_MEX)
647 mexAddFlops(6*numElements);
648 #elif defined(POWER_MEX)
649 mexAddFlops(flops);
650 #endif
651 #ifndef RELOP_MEX
652 removeZeroImag(zr, zi, ndim, (const int *)s, numElements, &zMat);
653 #endif
654 }
655 }
656 }
657
658
659 /***********************************************************
660 * *
661 * Tests to see if the vectors xdim and ydim are equal. *
662 * *
663 ***********************************************************/
664 bool allequal(int ndim, const int *xdim, const int *ydim)
665 {
666 int i;
667 bool result = true;
668
669 for (i = 0; i < ndim; i++)
670 result = result && (xdim[i] == ydim[i]);
671
672 return(result);
673 }
674
675
676 /******************************************************************************
677 * *
678 * Tests to see if every imaginary element is identically zero and, if so, *
679 * creates a new array which is real and copies the real elements to it. *
680 * *
681 ******************************************************************************/
682 void removeZeroImag(double *zr, double *zi, int ndim, const int *s,
683 int numElements, mxArray *plhs[])
684 {
685 double *zrend, *ziend, *zip, *z1p, *z2p;
686 bool allImZero = true;
687 mxArray *temp;
688
689 zip = zi;
690 ziend = zi + numElements;
691 while (zip < ziend)
692 {
693 allImZero = allImZero && (*zip++ == 0.0);
694 if (!allImZero)
695 return;
696 }
697
698 temp = mxCreateNumericArray(ndim, s, mxDOUBLE_CLASS, mxREAL);
699 z1p = zr;
700 z2p = mxGetPr(temp);
701 zrend = z1p + numElements;
702 while (z1p < zrend)
703 *z2p++ = *z1p++;
704 mxDestroyArray(plhs[0]);
705 plhs[0] = temp;
706 return;
707 }