00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if !defined (octave_mx_inlines_h)
00024 #define octave_mx_inlines_h 1
00025
00026 #include <cstddef>
00027
00028 #include "oct-cmplx.h"
00029
00030 #define VS_OP_FCN(F, OP) \
00031 template <class R, class V, class S> \
00032 inline void \
00033 F ## _vs (R *r, const V *v, size_t n, S s) \
00034 { \
00035 for (size_t i = 0; i < n; i++) \
00036 r[i] = v[i] OP s; \
00037 }
00038
00039 VS_OP_FCN (mx_inline_add, +)
00040 VS_OP_FCN (mx_inline_subtract, -)
00041 VS_OP_FCN (mx_inline_multiply, *)
00042 VS_OP_FCN (mx_inline_divide, /)
00043
00044 #define VS_OP(F, OP, R, V, S) \
00045 static inline R * \
00046 F (const V *v, size_t n, S s) \
00047 { \
00048 R *r = 0; \
00049 if (n > 0) \
00050 { \
00051 r = new R [n]; \
00052 F ## _vs (r, v, n, s); \
00053 } \
00054 return r; \
00055 }
00056
00057 #define VS_OPS(R, V, S) \
00058 VS_OP (mx_inline_add, +, R, V, S) \
00059 VS_OP (mx_inline_subtract, -, R, V, S) \
00060 VS_OP (mx_inline_multiply, *, R, V, S) \
00061 VS_OP (mx_inline_divide, /, R, V, S)
00062
00063 VS_OPS (double, double, double)
00064 VS_OPS (Complex, double, Complex)
00065 VS_OPS (Complex, Complex, double)
00066 VS_OPS (Complex, Complex, Complex)
00067
00068 #define SV_OP_FCN(F, OP) \
00069 template <class R, class S, class V> \
00070 inline void \
00071 F ## _sv (R *r, S s, const V *v, size_t n) \
00072 { \
00073 for (size_t i = 0; i < n; i++) \
00074 r[i] = s OP v[i]; \
00075 } \
00076
00077 SV_OP_FCN (mx_inline_add, +)
00078 SV_OP_FCN (mx_inline_subtract, -)
00079 SV_OP_FCN (mx_inline_multiply, *)
00080 SV_OP_FCN (mx_inline_divide, /)
00081
00082 #define SV_OP(F, OP, R, S, V) \
00083 static inline R * \
00084 F (S s, const V *v, size_t n) \
00085 { \
00086 R *r = 0; \
00087 if (n > 0) \
00088 { \
00089 r = new R [n]; \
00090 F ## _sv (r, s, v, n); \
00091 } \
00092 return r; \
00093 }
00094
00095 #define SV_OPS(R, S, V) \
00096 SV_OP (mx_inline_add, +, R, S, V) \
00097 SV_OP (mx_inline_subtract, -, R, S, V) \
00098 SV_OP (mx_inline_multiply, *, R, S, V) \
00099 SV_OP (mx_inline_divide, /, R, S, V)
00100
00101 SV_OPS (double, double, double)
00102 SV_OPS (Complex, double, Complex)
00103 SV_OPS (Complex, Complex, double)
00104 SV_OPS (Complex, Complex, Complex)
00105
00106 #define VV_OP_FCN(F, OP) \
00107 template <class R, class T1, class T2> \
00108 inline void \
00109 F ## _vv (R *r, const T1 *v1, const T2 *v2, size_t n) \
00110 { \
00111 for (size_t i = 0; i < n; i++) \
00112 r[i] = v1[i] OP v2[i]; \
00113 } \
00114
00115 VV_OP_FCN (mx_inline_add, +)
00116 VV_OP_FCN (mx_inline_subtract, -)
00117 VV_OP_FCN (mx_inline_multiply, *)
00118 VV_OP_FCN (mx_inline_divide, /)
00119
00120 #define VV_OP(F, OP, R, T1, T2) \
00121 static inline R * \
00122 F (const T1 *v1, const T2 *v2, size_t n) \
00123 { \
00124 R *r = 0; \
00125 if (n > 0) \
00126 { \
00127 r = new R [n]; \
00128 F ## _vv (r, v1, v2, n); \
00129 } \
00130 return r; \
00131 }
00132
00133 #define VV_OPS(R, T1, T2) \
00134 VV_OP (mx_inline_add, +, R, T1, T2) \
00135 VV_OP (mx_inline_subtract, -, R, T1, T2) \
00136 VV_OP (mx_inline_multiply, *, R, T1, T2) \
00137 VV_OP (mx_inline_divide, /, R, T1, T2)
00138
00139 VV_OPS (double, double, double)
00140 VV_OPS (Complex, double, Complex)
00141 VV_OPS (Complex, Complex, double)
00142 VV_OPS (Complex, Complex, Complex)
00143
00144 #define VS_OP2(F, OP, V, S) \
00145 static inline V * \
00146 F (V *v, size_t n, S s) \
00147 { \
00148 for (size_t i = 0; i < n; i++) \
00149 v[i] OP s; \
00150 return v; \
00151 }
00152
00153 #define VS_OP2S(V, S) \
00154 VS_OP2 (mx_inline_add2, +=, V, S) \
00155 VS_OP2 (mx_inline_subtract2, -=, V, S) \
00156 VS_OP2 (mx_inline_multiply2, *=, V, S) \
00157 VS_OP2 (mx_inline_divide2, /=, V, S) \
00158 VS_OP2 (mx_inline_copy, =, V, S)
00159
00160 VS_OP2S (double, double)
00161 VS_OP2S (Complex, double)
00162 VS_OP2S (Complex, Complex)
00163
00164 #define VV_OP2(F, OP, T1, T2) \
00165 static inline T1 * \
00166 F (T1 *v1, const T2 *v2, size_t n) \
00167 { \
00168 for (size_t i = 0; i < n; i++) \
00169 v1[i] OP v2[i]; \
00170 return v1; \
00171 }
00172
00173 #define VV_OP2S(T1, T2) \
00174 VV_OP2 (mx_inline_add2, +=, T1, T2) \
00175 VV_OP2 (mx_inline_subtract2, -=, T1, T2) \
00176 VV_OP2 (mx_inline_multiply2, *=, T1, T2) \
00177 VV_OP2 (mx_inline_divide2, /=, T1, T2) \
00178 VV_OP2 (mx_inline_copy, =, T1, T2)
00179
00180 VV_OP2S (double, double)
00181 VV_OP2S (Complex, double)
00182 VV_OP2S (Complex, Complex)
00183
00184 #define OP_EQ_FCN(T1, T2) \
00185 static inline bool \
00186 mx_inline_equal (const T1 *x, const T2 *y, size_t n) \
00187 { \
00188 for (size_t i = 0; i < n; i++) \
00189 if (x[i] != y[i]) \
00190 return false; \
00191 return true; \
00192 }
00193
00194 OP_EQ_FCN (bool, bool)
00195 OP_EQ_FCN (char, char)
00196 OP_EQ_FCN (double, double)
00197 OP_EQ_FCN (Complex, Complex)
00198
00199 #define OP_DUP_FCN(OP, F, R, T) \
00200 static inline R * \
00201 F (const T *x, size_t n) \
00202 { \
00203 R *r = 0; \
00204 if (n > 0) \
00205 { \
00206 r = new R [n]; \
00207 for (size_t i = 0; i < n; i++) \
00208 r[i] = OP (x[i]); \
00209 } \
00210 return r; \
00211 }
00212
00213 OP_DUP_FCN (, mx_inline_dup, double, double)
00214 OP_DUP_FCN (, mx_inline_dup, Complex, Complex)
00215
00216
00217
00218 OP_DUP_FCN (0.0 ==, mx_inline_not, double, double)
00219 OP_DUP_FCN (0.0 ==, mx_inline_not, double, Complex)
00220
00221 OP_DUP_FCN (, mx_inline_make_complex, Complex, double)
00222
00223 OP_DUP_FCN (-, mx_inline_change_sign, double, double)
00224 OP_DUP_FCN (-, mx_inline_change_sign, Complex, Complex)
00225
00226 OP_DUP_FCN (real, mx_inline_real_dup, double, Complex)
00227 OP_DUP_FCN (imag, mx_inline_imag_dup, double, Complex)
00228 OP_DUP_FCN (conj, mx_inline_conj_dup, Complex, Complex)
00229
00230
00231
00232 #define MX_CUMULATIVE_OP(RET_TYPE, ELT_TYPE, OP) \
00233 \
00234 int nr = rows (); \
00235 int nc = cols (); \
00236 \
00237 RET_TYPE retval (nr, nc); \
00238 \
00239 if (nr > 0 && nc > 0) \
00240 { \
00241 if ((nr == 1 && dim == -1) || dim == 1) \
00242 { \
00243 for (int i = 0; i < nr; i++) \
00244 { \
00245 ELT_TYPE t = elem (i, 0); \
00246 for (int j = 0; j < nc; j++) \
00247 { \
00248 retval.elem (i, j) = t; \
00249 if (j < nc - 1) \
00250 t OP elem (i, j+1); \
00251 } \
00252 } \
00253 } \
00254 else \
00255 { \
00256 for (int j = 0; j < nc; j++) \
00257 { \
00258 ELT_TYPE t = elem (0, j); \
00259 for (int i = 0; i < nr; i++) \
00260 { \
00261 retval.elem (i, j) = t; \
00262 if (i < nr - 1) \
00263 t OP elem (i+1, j); \
00264 } \
00265 } \
00266 } \
00267 } \
00268 \
00269 return retval
00270
00271 #define MX_BASE_REDUCTION_OP(RET_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, \
00272 MT_RESULT) \
00273 \
00274 int nr = rows (); \
00275 int nc = cols (); \
00276 \
00277 RET_TYPE retval; \
00278 \
00279 if (nr > 0 && nc > 0) \
00280 { \
00281 if ((nr == 1 && dim == -1) || dim == 1) \
00282 { \
00283 retval.resize (nr, 1); \
00284 for (int i = 0; i < nr; i++) \
00285 { \
00286 retval.elem (i, 0) = INIT_VAL; \
00287 for (int j = 0; j < nc; j++) \
00288 { \
00289 ROW_EXPR; \
00290 } \
00291 } \
00292 } \
00293 else \
00294 { \
00295 retval.resize (1, nc); \
00296 for (int j = 0; j < nc; j++) \
00297 { \
00298 retval.elem (0, j) = INIT_VAL; \
00299 for (int i = 0; i < nr; i++) \
00300 { \
00301 COL_EXPR; \
00302 } \
00303 } \
00304 } \
00305 } \
00306 else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \
00307 retval.resize (1, 1, MT_RESULT); \
00308 else if (nr == 0 && (dim == 0 || dim == -1)) \
00309 retval.resize (1, nc, MT_RESULT); \
00310 else if (nc == 0 && dim == 1) \
00311 retval.resize (nr, 1, MT_RESULT); \
00312 else \
00313 retval.resize (nr > 0, nc > 0); \
00314 \
00315 return retval
00316
00317 #define MX_REDUCTION_OP_ROW_EXPR(OP) \
00318 retval.elem (i, 0) OP elem (i, j)
00319
00320 #define MX_REDUCTION_OP_COL_EXPR(OP) \
00321 retval.elem (0, j) OP elem (i, j)
00322
00323 #define MX_REDUCTION_OP(RET_TYPE, OP, INIT_VAL, MT_RESULT) \
00324 MX_BASE_REDUCTION_OP (RET_TYPE, \
00325 MX_REDUCTION_OP_ROW_EXPR (OP), \
00326 MX_REDUCTION_OP_COL_EXPR (OP), \
00327 INIT_VAL, MT_RESULT)
00328
00329 #define MX_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \
00330 if (elem (i, j) TEST_OP 0.0) \
00331 { \
00332 retval.elem (i, 0) = TEST_TRUE_VAL; \
00333 break; \
00334 }
00335
00336 #define MX_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \
00337 if (elem (i, j) TEST_OP 0.0) \
00338 { \
00339 retval.elem (0, j) = TEST_TRUE_VAL; \
00340 break; \
00341 }
00342
00343 #define MX_ANY_ALL_OP(DIM, INIT_VAL, TEST_OP, TEST_TRUE_VAL) \
00344 MX_BASE_REDUCTION_OP (boolMatrix, \
00345 MX_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \
00346 MX_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \
00347 INIT_VAL, INIT_VAL)
00348
00349 #define MX_ALL_OP(DIM) MX_ANY_ALL_OP (DIM, true, ==, false)
00350
00351 #define MX_ANY_OP(DIM) MX_ANY_ALL_OP (DIM, false, !=, true)
00352
00353 #define MX_ND_ALL_EXPR elem (iter_idx) == 0
00354
00355 #define MX_ND_ANY_EXPR elem (iter_idx) != 0
00356
00357 #define MX_ND_ALL_EVAL(TEST_EXPR) \
00358 if (TEST_EXPR) \
00359 { \
00360 if (dim > -1) \
00361 iter_idx (dim) = 0; \
00362 retval (iter_idx) = 0; \
00363 break; \
00364 } \
00365 else \
00366 { \
00367 if (dim > -1) \
00368 iter_idx (dim) = 0; \
00369 retval (iter_idx) = 1; \
00370 } \
00371
00372 #define MX_ND_ANY_EVAL(TEST_EXPR) \
00373 if (TEST_EXPR) \
00374 { \
00375 if (dim > -1) \
00376 iter_idx (dim) = 0; \
00377 retval (iter_idx) = 1; \
00378 break; \
00379 }
00380
00381 #define MX_ND_REDUCTION(EVAL_EXPR, END_EXPR, VAL, ACC_DECL, \
00382 RET_TYPE, RET_ELT_TYPE) \
00383 \
00384 RET_TYPE retval; \
00385 \
00386 dim_vector dv = this->dims (); \
00387 \
00388 int empty = true; \
00389 \
00390 for (int i = 0; i < dv.length (); i++) \
00391 { \
00392 if (dv(i) > 0) \
00393 { \
00394 empty = false; \
00395 break; \
00396 } \
00397 } \
00398 \
00399 if (empty) \
00400 { \
00401 dim_vector retval_dv (1, 1); \
00402 retval.resize (retval_dv, VAL); \
00403 return retval; \
00404 } \
00405 \
00406 if (dim == -1) \
00407 { \
00408 for (int i = 0; i < dv.length (); i++) \
00409 { \
00410 if (dv (i) != 1) \
00411 { \
00412 dim = i; \
00413 break; \
00414 } \
00415 } \
00416 \
00417 if (dim == -1) \
00418 dim = 0; \
00419 } \
00420 \
00421 int squeezed = 0; \
00422 \
00423 for (int i = 0; i < dv.length (); i++) \
00424 { \
00425 if (dv(i) == 0) \
00426 { \
00427 squeezed = 1;\
00428 break;\
00429 } \
00430 } \
00431 \
00432 if (squeezed) \
00433 { \
00434 dv(dim) = 1; \
00435 \
00436 retval.resize (dv, VAL); \
00437 \
00438 return retval; \
00439 } \
00440 \
00441 \
00442 int dim_length = 1; \
00443 \
00444 \
00445 \
00446 \
00447 \
00448 if (dim >= dv.length ()) \
00449 dim = -1; \
00450 else \
00451 dim_length = dv(dim); \
00452 \
00453 if (dim > -1) \
00454 dv(dim) = 1; \
00455 \
00456 \
00457 int num_iter = (this->numel () / dim_length); \
00458 \
00459 \
00460 retval.resize (dv, RET_ELT_TYPE ()); \
00461 \
00462 Array<int> iter_idx (dv.length (), 0); \
00463 \
00464 \
00465 \
00466 \
00467 for (int j = 0; j < num_iter; j++) \
00468 { \
00469 ACC_DECL;\
00470 for (int i = 0; i < dim_length; i++) \
00471 { \
00472 if (dim > -1) \
00473 iter_idx(dim) = i; \
00474 \
00475 EVAL_EXPR; \
00476 } \
00477 \
00478 if (dim > -1) \
00479 iter_idx (dim) = 0; \
00480 \
00481 END_EXPR;\
00482 \
00483 increment_index (iter_idx, dv); \
00484 } \
00485 \
00486 retval.chop_trailing_singletons (); \
00487 \
00488 return retval
00489
00490 #define MX_ND_REAL_OP_REDUCTION(ASN_EXPR, INIT_VAL) \
00491 MX_ND_REDUCTION (acc ASN_EXPR, retval.elem (iter_idx) = acc, \
00492 INIT_VAL, double acc = INIT_VAL, NDArray, double)
00493
00494 #define MX_ND_COMPLEX_OP_REDUCTION(ASN_EXPR, INIT_VAL) \
00495 MX_ND_REDUCTION (acc ASN_EXPR, retval.elem (iter_idx) = acc, \
00496 INIT_VAL, Complex acc = INIT_VAL, ComplexNDArray, Complex)
00497
00498 #define MX_ND_ANY_ALL_REDUCTION(EVAL_EXPR, VAL) \
00499 MX_ND_REDUCTION (EVAL_EXPR, , VAL, , boolNDArray, bool)
00500
00501 #define MX_ND_CUMULATIVE_OP(RET_TYPE, ACC_TYPE, VAL, OP) \
00502 RET_TYPE retval; \
00503 \
00504 dim_vector dv = this->dims (); \
00505 \
00506 int empty = true; \
00507 \
00508 \
00509 if (dim >= dv.length ()) \
00510 { \
00511 retval = RET_TYPE (*this); \
00512 return retval; \
00513 } \
00514 \
00515 \
00516 for (int i = 0; i < dv.length (); i++) \
00517 { \
00518 if (dv(i) > 0) \
00519 { \
00520 empty = false; \
00521 break; \
00522 } \
00523 } \
00524 \
00525 if (empty) \
00526 { \
00527 retval.resize (dv); \
00528 return retval; \
00529 } \
00530 \
00531 \
00532 if (dim == -1) \
00533 { \
00534 for (int i = 0; i < dv.length (); i++) \
00535 { \
00536 if (dv (i) != 1) \
00537 { \
00538 dim = i; \
00539 break; \
00540 } \
00541 } \
00542 \
00543 if (dim == -1) \
00544 dim = 0; \
00545 } \
00546 \
00547 \
00548 \
00549 int squeezed = 0; \
00550 \
00551 for (int i = 0; i < dv.length (); i++) \
00552 { \
00553 if (dv(i) == 0) \
00554 { \
00555 squeezed = 1; \
00556 break; \
00557 } \
00558 } \
00559 \
00560 if (squeezed) \
00561 { \
00562 retval.resize (dv); \
00563 return retval; \
00564 } \
00565 \
00566 \
00567 retval.resize (dv, VAL); \
00568 \
00569 \
00570 int dim_length = 1; \
00571 \
00572 dim_length = dv (dim); \
00573 \
00574 dv (dim) = 1; \
00575 \
00576 \
00577 int num_iter = (this->numel () / dim_length); \
00578 \
00579 Array<int> iter_idx (dv.length (), 0); \
00580 \
00581 \
00582 \
00583 \
00584 for (int j = 0; j < num_iter; j++) \
00585 { \
00586 for (int i = 0; i < dim_length; i++) \
00587 { \
00588 if (i > 0) \
00589 { \
00590 iter_idx (dim) = i - 1; \
00591 \
00592 ACC_TYPE prev_sum = retval (iter_idx); \
00593 \
00594 iter_idx (dim) = i; \
00595 \
00596 retval (iter_idx) = elem (iter_idx) OP prev_sum; \
00597 } \
00598 else \
00599 retval (iter_idx) = elem (iter_idx); \
00600 } \
00601 \
00602 if (dim > -1) \
00603 iter_idx (dim) = 0; \
00604 \
00605 increment_index (iter_idx, dv); \
00606 } \
00607 \
00608 return retval
00609
00610 #endif
00611
00612
00613
00614
00615
00616