Mercurial > hg > camir-aes2014

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