1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING.  If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if ! defined (octave_mx_op_defs_h)
27 #define octave_mx_op_defs_h 1
28 
29 #include "octave-config.h"
30 
31 #include "lo-array-errwarn.h"
32 #include "mx-op-decl.h"
33 #include "mx-inlines.cc"
34 
35 #define SNANCHK(s)                              \
36   if (octave::math::isnan (s))                  \
37     octave::err_nan_to_logical_conversion ()
38 
39 #define MNANCHK(m, MT)                          \
40   if (do_mx_check (m, mx_inline_any_nan<MT>))   \
41     octave::err_nan_to_logical_conversion ()
42 
43 // vector by scalar operations.
44 
45 #define VS_BIN_OP(R, F, OP, V, S)                                       \
46   R                                                                     \
47   F (const V& v, const S& s)                                            \
48   {                                                                     \
49     return do_ms_binary_op<R::element_type, V::element_type, S> (v, s, OP); \
50   }
51 
52 #define VS_BIN_OPS(R, V, S)                             \
53   VS_BIN_OP (R, operator +, mx_inline_add, V, S)        \
54   VS_BIN_OP (R, operator -, mx_inline_sub, V, S)        \
55   VS_BIN_OP (R, operator *, mx_inline_mul, V, S)        \
56   VS_BIN_OP (R, operator /, mx_inline_div, V, S)
57 
58 // scalar by vector by operations.
59 
60 #define SV_BIN_OP(R, F, OP, S, V)                                       \
61   R                                                                     \
62   F (const S& s, const V& v)                                            \
63   {                                                                     \
64     return do_sm_binary_op<R::element_type, S, V::element_type> (s, v, OP); \
65   }
66 
67 #define SV_BIN_OPS(R, S, V)                             \
68   SV_BIN_OP (R, operator +, mx_inline_add, S, V)        \
69   SV_BIN_OP (R, operator -, mx_inline_sub, S, V)        \
70   SV_BIN_OP (R, operator *, mx_inline_mul, S, V)        \
71   SV_BIN_OP (R, operator /, mx_inline_div, S, V)
72 
73 // vector by vector operations.
74 
75 #define VV_BIN_OP(R, F, OP, V1, V2)                                     \
76   R                                                                     \
77   F (const V1& v1, const V2& v2)                                        \
78   {                                                                     \
79     return do_mm_binary_op<R::element_type, V1::element_type, V2::element_type> (v1, v2, OP, OP, OP, #F); \
80   }
81 
82 #define VV_BIN_OPS(R, V1, V2)                           \
83   VV_BIN_OP (R, operator +, mx_inline_add, V1, V2)      \
84   VV_BIN_OP (R, operator -, mx_inline_sub, V1, V2)      \
85   VV_BIN_OP (R, product,    mx_inline_mul, V1, V2)      \
86   VV_BIN_OP (R, quotient,   mx_inline_div, V1, V2)
87 
88 // matrix by scalar operations.
89 
90 #define MS_BIN_OP(R, OP, M, S, F)                                       \
91   R                                                                     \
92   OP (const M& m, const S& s)                                           \
93   {                                                                     \
94     return do_ms_binary_op<R::element_type, M::element_type, S> (m, s, F); \
95   }
96 
97 #define MS_BIN_OPS(R, M, S)                             \
98   MS_BIN_OP (R, operator +, M, S, mx_inline_add)        \
99   MS_BIN_OP (R, operator -, M, S, mx_inline_sub)        \
100   MS_BIN_OP (R, operator *, M, S, mx_inline_mul)        \
101   MS_BIN_OP (R, operator /, M, S, mx_inline_div)
102 
103 #define MS_CMP_OP(F, OP, M, S)                                          \
104   boolMatrix                                                            \
105   F (const M& m, const S& s)                                            \
106   {                                                                     \
107     return do_ms_binary_op<bool, M::element_type, S> (m, s, OP);        \
108   }
109 
110 #define MS_CMP_OPS(M, S)                        \
111   MS_CMP_OP (mx_el_lt, mx_inline_lt, M, S)      \
112   MS_CMP_OP (mx_el_le, mx_inline_le, M, S)      \
113   MS_CMP_OP (mx_el_ge, mx_inline_ge, M, S)      \
114   MS_CMP_OP (mx_el_gt, mx_inline_gt, M, S)      \
115   MS_CMP_OP (mx_el_eq, mx_inline_eq, M, S)      \
116   MS_CMP_OP (mx_el_ne, mx_inline_ne, M, S)
117 
118 #define MS_BOOL_OP(F, OP, M, S)                                         \
119   boolMatrix                                                            \
120   F (const M& m, const S& s)                                            \
121   {                                                                     \
122     MNANCHK (m, M::element_type);                                       \
123     SNANCHK (s);                                                        \
124     return do_ms_binary_op<bool, M::element_type, S> (m, s, OP);        \
125   }
126 
127 #define MS_BOOL_OPS(M, S)                       \
128   MS_BOOL_OP (mx_el_and, mx_inline_and, M, S)   \
129   MS_BOOL_OP (mx_el_or,  mx_inline_or,  M, S)
130 
131 // scalar by matrix operations.
132 
133 #define SM_BIN_OP(R, OP, S, M, F)                                       \
134   R                                                                     \
135   OP (const S& s, const M& m)                                           \
136   {                                                                     \
137     return do_sm_binary_op<R::element_type, S, M::element_type> (s, m, F); \
138   }
139 
140 #define SM_BIN_OPS(R, S, M)                             \
141   SM_BIN_OP (R, operator +, S, M, mx_inline_add)        \
142   SM_BIN_OP (R, operator -, S, M, mx_inline_sub)        \
143   SM_BIN_OP (R, operator *, S, M, mx_inline_mul)        \
144   SM_BIN_OP (R, operator /, S, M, mx_inline_div)
145 
146 #define SM_CMP_OP(F, OP, S, M)                                          \
147   boolMatrix                                                            \
148   F (const S& s, const M& m)                                            \
149   {                                                                     \
150     return do_sm_binary_op<bool, S, M::element_type> (s, m, OP);        \
151   }
152 
153 #define SM_CMP_OPS(S, M)                        \
154   SM_CMP_OP (mx_el_lt, mx_inline_lt, S, M)      \
155   SM_CMP_OP (mx_el_le, mx_inline_le, S, M)      \
156   SM_CMP_OP (mx_el_ge, mx_inline_ge, S, M)      \
157   SM_CMP_OP (mx_el_gt, mx_inline_gt, S, M)      \
158   SM_CMP_OP (mx_el_eq, mx_inline_eq, S, M)      \
159   SM_CMP_OP (mx_el_ne, mx_inline_ne, S, M)
160 
161 #define SM_BOOL_OP(F, OP, S, M)                                         \
162   boolMatrix                                                            \
163   F (const S& s, const M& m)                                            \
164   {                                                                     \
165     SNANCHK (s);                                                        \
166     MNANCHK (m, M::element_type);                                       \
167     return do_sm_binary_op<bool, S, M::element_type> (s, m, OP);        \
168   }
169 
170 #define SM_BOOL_OPS(S, M)                       \
171   SM_BOOL_OP (mx_el_and, mx_inline_and, S, M)   \
172   SM_BOOL_OP (mx_el_or,  mx_inline_or,  S, M)
173 
174 // matrix by matrix operations.
175 
176 #define MM_BIN_OP(R, OP, M1, M2, F)                                     \
177   R                                                                     \
178   OP (const M1& m1, const M2& m2)                                       \
179   {                                                                     \
180     return do_mm_binary_op<R::element_type, M1::element_type, M2::element_type> (m1, m2, F, F, F, #OP); \
181   }
182 
183 #define MM_BIN_OPS(R, M1, M2)                           \
184   MM_BIN_OP (R, operator +, M1, M2, mx_inline_add)      \
185   MM_BIN_OP (R, operator -, M1, M2, mx_inline_sub)      \
186   MM_BIN_OP (R, product,    M1, M2, mx_inline_mul)      \
187   MM_BIN_OP (R, quotient,   M1, M2, mx_inline_div)
188 
189 #define MM_CMP_OP(F, OP, M1, M2)                                        \
190   boolMatrix                                                            \
191   F (const M1& m1, const M2& m2)                                        \
192   {                                                                     \
193     return do_mm_binary_op<bool, M1::element_type, M2::element_type> (m1, m2, OP, OP, OP, #F); \
194   }
195 
196 #define MM_CMP_OPS(M1, M2)                      \
197   MM_CMP_OP (mx_el_lt, mx_inline_lt, M1, M2)    \
198   MM_CMP_OP (mx_el_le, mx_inline_le, M1, M2)    \
199   MM_CMP_OP (mx_el_ge, mx_inline_ge, M1, M2)    \
200   MM_CMP_OP (mx_el_gt, mx_inline_gt, M1, M2)    \
201   MM_CMP_OP (mx_el_eq, mx_inline_eq, M1, M2)    \
202   MM_CMP_OP (mx_el_ne, mx_inline_ne, M1, M2)
203 
204 #define MM_BOOL_OP(F, OP, M1, M2)                                       \
205   boolMatrix                                                            \
206   F (const M1& m1, const M2& m2)                                        \
207   {                                                                     \
208     MNANCHK (m1, M1::element_type);                                     \
209     MNANCHK (m2, M2::element_type);                                     \
210     return do_mm_binary_op<bool, M1::element_type, M2::element_type> (m1, m2, OP, OP, OP, #F); \
211   }
212 
213 #define MM_BOOL_OPS(M1, M2)                     \
214   MM_BOOL_OP (mx_el_and, mx_inline_and, M1, M2) \
215   MM_BOOL_OP (mx_el_or,  mx_inline_or,  M1, M2)
216 
217 // N-D matrix by scalar operations.
218 
219 #define NDS_BIN_OP(R, OP, ND, S, F)                                     \
220   R                                                                     \
221   OP (const ND& m, const S& s)                                          \
222   {                                                                     \
223     return do_ms_binary_op<R::element_type, ND::element_type, S> (m, s, F); \
224   }
225 
226 #define NDS_BIN_OPS(R, ND, S)                           \
227   NDS_BIN_OP (R, operator +, ND, S, mx_inline_add)      \
228   NDS_BIN_OP (R, operator -, ND, S, mx_inline_sub)      \
229   NDS_BIN_OP (R, operator *, ND, S, mx_inline_mul)      \
230   NDS_BIN_OP (R, operator /, ND, S, mx_inline_div)
231 
232 #define NDS_CMP_OP(F, OP, ND, S)                                        \
233   boolNDArray                                                           \
234   F (const ND& m, const S& s)                                           \
235   {                                                                     \
236     return do_ms_binary_op<bool, ND::element_type, S> (m, s, OP);       \
237   }
238 
239 #define NDS_CMP_OPS(ND, S)                      \
240   NDS_CMP_OP (mx_el_lt, mx_inline_lt, ND, S)    \
241   NDS_CMP_OP (mx_el_le, mx_inline_le, ND, S)    \
242   NDS_CMP_OP (mx_el_ge, mx_inline_ge, ND, S)    \
243   NDS_CMP_OP (mx_el_gt, mx_inline_gt, ND, S)    \
244   NDS_CMP_OP (mx_el_eq, mx_inline_eq, ND, S)    \
245   NDS_CMP_OP (mx_el_ne, mx_inline_ne, ND, S)
246 
247 #define NDS_BOOL_OP(F, OP, ND, S)                                       \
248   boolNDArray                                                           \
249   F (const ND& m, const S& s)                                           \
250   {                                                                     \
251     MNANCHK (m, ND::element_type);                                      \
252     SNANCHK (s);                                                        \
253     return do_ms_binary_op<bool, ND::element_type, S> (m, s, OP);       \
254   }
255 
256 #define NDS_BOOL_OPS(ND, S)                             \
257   NDS_BOOL_OP (mx_el_and,     mx_inline_and,     ND, S) \
258   NDS_BOOL_OP (mx_el_or,      mx_inline_or,      ND, S) \
259   NDS_BOOL_OP (mx_el_not_and, mx_inline_not_and, ND, S) \
260   NDS_BOOL_OP (mx_el_not_or,  mx_inline_not_or,  ND, S) \
261   NDS_BOOL_OP (mx_el_and_not, mx_inline_and_not, ND, S) \
262   NDS_BOOL_OP (mx_el_or_not,  mx_inline_or_not,  ND, S)
263 
264 // scalar by N-D matrix operations.
265 
266 #define SND_BIN_OP(R, OP, S, ND, F)                                     \
267   R                                                                     \
268   OP (const S& s, const ND& m)                                          \
269   {                                                                     \
270     return do_sm_binary_op<R::element_type, S, ND::element_type> (s, m, F); \
271   }
272 
273 #define SND_BIN_OPS(R, S, ND)                           \
274   SND_BIN_OP (R, operator +, S, ND, mx_inline_add)      \
275   SND_BIN_OP (R, operator -, S, ND, mx_inline_sub)      \
276   SND_BIN_OP (R, operator *, S, ND, mx_inline_mul)      \
277   SND_BIN_OP (R, operator /, S, ND, mx_inline_div)
278 
279 #define SND_CMP_OP(F, OP, S, ND)                                        \
280   boolNDArray                                                           \
281   F (const S& s, const ND& m)                                           \
282   {                                                                     \
283     return do_sm_binary_op<bool, S, ND::element_type> (s, m, OP);       \
284   }
285 
286 #define SND_CMP_OPS(S, ND)                      \
287   SND_CMP_OP (mx_el_lt, mx_inline_lt, S, ND)    \
288   SND_CMP_OP (mx_el_le, mx_inline_le, S, ND)    \
289   SND_CMP_OP (mx_el_ge, mx_inline_ge, S, ND)    \
290   SND_CMP_OP (mx_el_gt, mx_inline_gt, S, ND)    \
291   SND_CMP_OP (mx_el_eq, mx_inline_eq, S, ND)    \
292   SND_CMP_OP (mx_el_ne, mx_inline_ne, S, ND)
293 
294 #define SND_BOOL_OP(F, OP, S, ND)                                       \
295   boolNDArray                                                           \
296   F (const S& s, const ND& m)                                           \
297   {                                                                     \
298     SNANCHK (s);                                                        \
299     MNANCHK (m, ND::element_type);                                      \
300     return do_sm_binary_op<bool, S, ND::element_type> (s, m, OP);       \
301   }
302 
303 #define SND_BOOL_OPS(S, ND)                             \
304   SND_BOOL_OP (mx_el_and,     mx_inline_and,     S, ND) \
305   SND_BOOL_OP (mx_el_or,      mx_inline_or,      S, ND) \
306   SND_BOOL_OP (mx_el_not_and, mx_inline_not_and, S, ND) \
307   SND_BOOL_OP (mx_el_not_or,  mx_inline_not_or,  S, ND) \
308   SND_BOOL_OP (mx_el_and_not, mx_inline_and_not, S, ND) \
309   SND_BOOL_OP (mx_el_or_not,  mx_inline_or_not,  S, ND)
310 
311 // N-D matrix by N-D matrix operations.
312 
313 #define NDND_BIN_OP(R, OP, ND1, ND2, F)                                 \
314   R                                                                     \
315   OP (const ND1& m1, const ND2& m2)                                     \
316   {                                                                     \
317     return do_mm_binary_op<R::element_type, ND1::element_type, ND2::element_type> (m1, m2, F, F, F, #OP); \
318   }
319 
320 #define NDND_BIN_OPS(R, ND1, ND2)                       \
321   NDND_BIN_OP (R, operator +, ND1, ND2, mx_inline_add)  \
322   NDND_BIN_OP (R, operator -, ND1, ND2, mx_inline_sub)  \
323   NDND_BIN_OP (R, product,    ND1, ND2, mx_inline_mul)  \
324   NDND_BIN_OP (R, quotient,   ND1, ND2, mx_inline_div)
325 
326 #define NDND_CMP_OP(F, OP, ND1, ND2)                                    \
327   boolNDArray                                                           \
328   F (const ND1& m1, const ND2& m2)                                      \
329   {                                                                     \
330     return do_mm_binary_op<bool, ND1::element_type, ND2::element_type> (m1, m2, OP, OP, OP, #F); \
331   }
332 
333 #define NDND_CMP_OPS(ND1, ND2)                          \
334   NDND_CMP_OP (mx_el_lt, mx_inline_lt, ND1, ND2)        \
335   NDND_CMP_OP (mx_el_le, mx_inline_le, ND1, ND2)        \
336   NDND_CMP_OP (mx_el_ge, mx_inline_ge, ND1, ND2)        \
337   NDND_CMP_OP (mx_el_gt, mx_inline_gt, ND1, ND2)        \
338   NDND_CMP_OP (mx_el_eq, mx_inline_eq, ND1, ND2)        \
339   NDND_CMP_OP (mx_el_ne, mx_inline_ne, ND1, ND2)
340 
341 #define NDND_BOOL_OP(F, OP, ND1, ND2)                                   \
342   boolNDArray                                                           \
343   F (const ND1& m1, const ND2& m2)                                      \
344   {                                                                     \
345     MNANCHK (m1, ND1::element_type);                                    \
346     MNANCHK (m2, ND2::element_type);                                    \
347     return do_mm_binary_op<bool, ND1::element_type, ND2::element_type> (m1, m2, OP, OP, OP, #F); \
348   }
349 
350 #define NDND_BOOL_OPS(ND1, ND2)                                 \
351   NDND_BOOL_OP (mx_el_and,     mx_inline_and,     ND1, ND2)     \
352   NDND_BOOL_OP (mx_el_or,      mx_inline_or,      ND1, ND2)     \
353   NDND_BOOL_OP (mx_el_not_and, mx_inline_not_and, ND1, ND2)     \
354   NDND_BOOL_OP (mx_el_not_or,  mx_inline_not_or,  ND1, ND2)     \
355   NDND_BOOL_OP (mx_el_and_not, mx_inline_and_not, ND1, ND2)     \
356   NDND_BOOL_OP (mx_el_or_not,  mx_inline_or_not,  ND1, ND2)
357 
358 // scalar by diagonal matrix operations.
359 
360 #define SDM_BIN_OP(R, OP, S, DM)                        \
361   R                                                     \
362   operator OP (const S& s, const DM& dm)                \
363   {                                                     \
364     R r (dm.rows (), dm.cols ());                       \
365                                                         \
366     for (octave_idx_type i = 0; i < dm.length (); i++)  \
367       r.dgxelem (i) = s OP dm.dgelem (i);               \
368                                                         \
369     return r;                                           \
370   }
371 
372 #define SDM_BIN_OPS(R, S, DM)                   \
373   SDM_BIN_OP (R, *, S, DM)
374 
375 // diagonal matrix by scalar operations.
376 
377 #define DMS_BIN_OP(R, OP, DM, S)                        \
378   R                                                     \
379   operator OP (const DM& dm, const S& s)                \
380   {                                                     \
381     R r (dm.rows (), dm.cols ());                       \
382                                                         \
383     for (octave_idx_type i = 0; i < dm.length (); i++)  \
384       r.dgxelem (i) = dm.dgelem (i) OP s;               \
385                                                         \
386     return r;                                           \
387   }
388 
389 #define DMS_BIN_OPS(R, DM, S)                   \
390   DMS_BIN_OP (R, *, DM, S)                      \
391   DMS_BIN_OP (R, /, DM, S)
392 
393 // matrix by diagonal matrix operations.
394 
395 #define MDM_BIN_OP(R, OP, M, DM, OPEQ)                          \
396   R                                                             \
397   OP (const M& m, const DM& dm)                                 \
398   {                                                             \
399     R r;                                                        \
400                                                                 \
401     octave_idx_type m_nr = m.rows ();                           \
402     octave_idx_type m_nc = m.cols ();                           \
403                                                                 \
404     octave_idx_type dm_nr = dm.rows ();                         \
405     octave_idx_type dm_nc = dm.cols ();                         \
406                                                                 \
407     if (m_nr != dm_nr || m_nc != dm_nc)                         \
408       octave::err_nonconformant (#OP, m_nr, m_nc, dm_nr, dm_nc);        \
409                                                                 \
410     r.resize (m_nr, m_nc);                                      \
411                                                                 \
412     if (m_nr > 0 && m_nc > 0)                                   \
413       {                                                         \
414         r = R (m);                                              \
415                                                                 \
416         octave_idx_type len = dm.length ();                     \
417                                                                 \
418         for (octave_idx_type i = 0; i < len; i++)               \
419           r.elem (i, i) OPEQ dm.elem (i, i);                    \
420       }                                                         \
421                                                                 \
422     return r;                                                   \
423   }
424 
425 #define MDM_MULTIPLY_OP(R, M, DM)                                       \
426   R                                                                     \
427   operator * (const M& m, const DM& dm)                                 \
428   {                                                                     \
429     R r;                                                                \
430                                                                         \
431     R::element_type r_zero = R::element_type ();                        \
432                                                                         \
433     octave_idx_type m_nr = m.rows ();                                   \
434     octave_idx_type m_nc = m.cols ();                                   \
435                                                                         \
436     octave_idx_type dm_nr = dm.rows ();                                 \
437     octave_idx_type dm_nc = dm.cols ();                                 \
438                                                                         \
439     if (m_nc != dm_nr)                                                  \
440       octave::err_nonconformant ("operator *", m_nr, m_nc, dm_nr, dm_nc);       \
441                                                                         \
442     r = R (m_nr, dm_nc);                                                \
443     R::element_type *rd = r.fortran_vec ();                             \
444     const M::element_type *md = m.data ();                              \
445     const DM::element_type *dd = dm.data ();                            \
446                                                                         \
447     octave_idx_type len = dm.length ();                                 \
448     for (octave_idx_type i = 0; i < len; i++)                           \
449       {                                                                 \
450         mx_inline_mul (m_nr, rd, md, dd[i]);                            \
451         rd += m_nr; md += m_nr;                                         \
452       }                                                                 \
453     mx_inline_fill (m_nr * (dm_nc - len), rd, r_zero);                  \
454                                                                         \
455     return r;                                                           \
456   }
457 
458 #define MDM_BIN_OPS(R, M, DM)                   \
459   MDM_BIN_OP (R, operator +, M, DM, +=)         \
460   MDM_BIN_OP (R, operator -, M, DM, -=)         \
461   MDM_MULTIPLY_OP (R, M, DM)
462 
463 // diagonal matrix by matrix operations.
464 
465 #define DMM_BIN_OP(R, OP, DM, M, OPEQ, PREOP)                   \
466   R                                                             \
467   OP (const DM& dm, const M& m)                                 \
468   {                                                             \
469     R r;                                                        \
470                                                                 \
471     octave_idx_type dm_nr = dm.rows ();                         \
472     octave_idx_type dm_nc = dm.cols ();                         \
473                                                                 \
474     octave_idx_type m_nr = m.rows ();                           \
475     octave_idx_type m_nc = m.cols ();                           \
476                                                                 \
477     if (dm_nr != m_nr || dm_nc != m_nc)                         \
478       octave::err_nonconformant (#OP, dm_nr, dm_nc, m_nr, m_nc);        \
479     else                                                        \
480       {                                                         \
481         if (m_nr > 0 && m_nc > 0)                               \
482           {                                                     \
483             r = R (PREOP m);                                    \
484                                                                 \
485             octave_idx_type len = dm.length ();                 \
486                                                                 \
487             for (octave_idx_type i = 0; i < len; i++)           \
488               r.elem (i, i) OPEQ dm.elem (i, i);                \
489           }                                                     \
490         else                                                    \
491           r.resize (m_nr, m_nc);                                \
492       }                                                         \
493                                                                 \
494     return r;                                                   \
495   }
496 
497 #define DMM_MULTIPLY_OP(R, DM, M)                                       \
498   R                                                                     \
499   operator * (const DM& dm, const M& m)                                 \
500   {                                                                     \
501     R r;                                                                \
502                                                                         \
503     R::element_type r_zero = R::element_type ();                        \
504                                                                         \
505     octave_idx_type dm_nr = dm.rows ();                                 \
506     octave_idx_type dm_nc = dm.cols ();                                 \
507                                                                         \
508     octave_idx_type m_nr = m.rows ();                                   \
509     octave_idx_type m_nc = m.cols ();                                   \
510                                                                         \
511     if (dm_nc != m_nr)                                                  \
512       octave::err_nonconformant ("operator *", dm_nr, dm_nc, m_nr, m_nc);       \
513                                                                         \
514     r = R (dm_nr, m_nc);                                                \
515     R::element_type *rd = r.fortran_vec ();                             \
516     const M::element_type *md = m.data ();                              \
517     const DM::element_type *dd = dm.data ();                            \
518                                                                         \
519     octave_idx_type len = dm.length ();                                 \
520     for (octave_idx_type i = 0; i < m_nc; i++)                          \
521       {                                                                 \
522         mx_inline_mul (len, rd, md, dd);                                \
523         rd += len; md += m_nr;                                          \
524         mx_inline_fill (dm_nr - len, rd, r_zero);                       \
525         rd += dm_nr - len;                                              \
526       }                                                                 \
527                                                                         \
528     return r;                                                           \
529   }
530 
531 #define DMM_BIN_OPS(R, DM, M)                   \
532   DMM_BIN_OP (R, operator +, DM, M, +=, )       \
533   DMM_BIN_OP (R, operator -, DM, M, +=, -)      \
534   DMM_MULTIPLY_OP (R, DM, M)
535 
536 // diagonal matrix by diagonal matrix operations.
537 
538 #define DMDM_BIN_OP(R, OP, DM1, DM2, F)                                 \
539   R                                                                     \
540   OP (const DM1& dm1, const DM2& dm2)                                   \
541   {                                                                     \
542     R r;                                                                \
543                                                                         \
544     octave_idx_type dm1_nr = dm1.rows ();                               \
545     octave_idx_type dm1_nc = dm1.cols ();                               \
546                                                                         \
547     octave_idx_type dm2_nr = dm2.rows ();                               \
548     octave_idx_type dm2_nc = dm2.cols ();                               \
549                                                                         \
550     if (dm1_nr != dm2_nr || dm1_nc != dm2_nc)                           \
551       octave::err_nonconformant (#OP, dm1_nr, dm1_nc, dm2_nr, dm2_nc);          \
552                                                                         \
553     r.resize (dm1_nr, dm1_nc);                                          \
554                                                                         \
555     if (dm1_nr > 0 && dm1_nc > 0)                                       \
556       F (dm1.length (), r.fortran_vec (), dm1.data (), dm2.data ());    \
557                                                                         \
558     return r;                                                           \
559   }
560 
561 #define DMDM_BIN_OPS(R, DM1, DM2)                       \
562   DMDM_BIN_OP (R, operator +, DM1, DM2, mx_inline_add)  \
563   DMDM_BIN_OP (R, operator -, DM1, DM2, mx_inline_sub)  \
564   DMDM_BIN_OP (R, product,    DM1, DM2, mx_inline_mul)
565 
566 // scalar by N-D array min/max ops
567 
568 #define SND_MINMAX_FCN(FCN, OP, T, S)                                   \
569   T                                                                     \
570   FCN (S d, const T& m)                                                 \
571   {                                                                     \
572     return do_sm_binary_op<T::element_type, S, T::element_type> (d, m, mx_inline_x##FCN); \
573   }
574 
575 #define NDS_MINMAX_FCN(FCN, OP, T, S)                                   \
576   T                                                                     \
577   FCN (const T& m, S d)                                                 \
578   {                                                                     \
579     return do_ms_binary_op<T::element_type, T::element_type, S> (m, d, mx_inline_x##FCN); \
580   }
581 
582 #define NDND_MINMAX_FCN(FCN, OP, T, S)                                  \
583   T                                                                     \
584   FCN (const T& a, const T& b)                                          \
585   {                                                                     \
586     return do_mm_binary_op<T::element_type, T::element_type, T::element_type> (a, b, mx_inline_x##FCN, mx_inline_x##FCN, mx_inline_x##FCN, #FCN); \
587   }
588 
589 #define MINMAX_FCNS(T, S)                       \
590   SND_MINMAX_FCN (min, <, T, S)                 \
591   NDS_MINMAX_FCN (min, <, T, S)                 \
592   NDND_MINMAX_FCN (min, <, T, S)                \
593   SND_MINMAX_FCN (max, >, T, S)                 \
594   NDS_MINMAX_FCN (max, >, T, S)                 \
595   NDND_MINMAX_FCN (max, >, T, S)
596 
597 // permutation matrix by matrix ops and vice versa
598 
599 #define PMM_MULTIPLY_OP(PM, M)                                          \
600   M operator * (const PM& p, const M& x)                                \
601   {                                                                     \
602     octave_idx_type nr = x.rows ();                                     \
603     octave_idx_type nc = x.columns ();                                  \
604     M result;                                                           \
605     if (p.columns () != nr)                                             \
606       octave::err_nonconformant ("operator *", p.rows (), p.columns (), nr, nc); \
607     else                                                                \
608       {                                                                 \
609         result = M (nr, nc);                                            \
610         result.assign (p.col_perm_vec (), idx_vector::colon, x);        \
611       }                                                                 \
612                                                                         \
613     return result;                                                      \
614   }
615 
616 #define MPM_MULTIPLY_OP(M, PM)                                          \
617   M operator * (const M& x, const PM& p)                                \
618   {                                                                     \
619     octave_idx_type nr = x.rows ();                                     \
620     octave_idx_type nc = x.columns ();                                  \
621     M result;                                                           \
622     if (p.rows () != nc)                                                \
623       octave::err_nonconformant ("operator *", nr, nc, p.rows (), p.columns ()); \
624                                                                         \
625     result = x.index (idx_vector::colon, p.col_perm_vec ());            \
626                                                                         \
627     return result;                                                      \
628   }
629 
630 #define PMM_BIN_OPS(R, PM, M)                   \
631   PMM_MULTIPLY_OP(PM, M);
632 
633 #define MPM_BIN_OPS(R, M, PM)                   \
634   MPM_MULTIPLY_OP(M, PM);
635 
636 #define NDND_MAPPER_BODY(R, NAME)               \
637   R retval (dims ());                           \
638   octave_idx_type n = numel ();                 \
639   for (octave_idx_type i = 0; i < n; i++)       \
640     retval.xelem (i) = NAME (elem (i));         \
641   return retval;
642 
643 #endif
644