bug 13918: Unmanaged operations on hypermatrix did not call overload functions.
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_subtraction.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2008-2008 - DIGITEO - Antoine ELIAS
4 *  Copyright (C) 2010-2010 - DIGITEO - Bruno JOFRET
5 *
6 *  This file must be used under the terms of the CeCILL.
7 *  This source file is licensed as described in the file COPYING, which
8 *  you should have received as part of this distribution.  The terms
9 *  are also available at
10 *  http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 *
12 */
13
14 #include <algorithm>
15
16 #include "types_subtraction.hxx"
17 #include "types_opposite.hxx"
18 #include "double.hxx"
19 #include "polynom.hxx"
20 #include "int.hxx"
21 #include "sparse.hxx"
22 #include "generic_operations.hxx"
23
24 extern "C"
25 {
26 #include "elem_common.h" //dset
27 #include "Sciwarning.h" //Sciwarning
28 }
29
30 using namespace types;
31
32 //define arrays on operation functions
33 static sub_function pSubfunction[types::InternalType::IdLast][types::InternalType::IdLast] = {NULL};
34
35 void fillSubtractFunction()
36 {
37 #define scilab_fill_sub(id1, id2, func, typeIn1, typeIn2, typeOut) \
38     pSubfunction[types::InternalType::Id ## id1][types::InternalType::Id ## id2] = (sub_function)&sub_##func<typeIn1, typeIn2, typeOut>
39
40     //Double
41     //Matrix - Matrix
42     scilab_fill_sub(Double, Double, M_M, Double, Double, Double);
43     scilab_fill_sub(Double, Int8, M_M, Double, Int8, Int8);
44     scilab_fill_sub(Double, UInt8, M_M, Double, UInt8, UInt8);
45     scilab_fill_sub(Double, Int16, M_M, Double, Int16, Int16);
46     scilab_fill_sub(Double, UInt16, M_M, Double, UInt16, UInt16);
47     scilab_fill_sub(Double, Int32, M_M, Double, Int32, Int32);
48     scilab_fill_sub(Double, UInt32, M_M, Double, UInt32, UInt32);
49     scilab_fill_sub(Double, Int64, M_M, Double, Int64, Int64);
50     scilab_fill_sub(Double, UInt64, M_M, Double, UInt64, UInt64);
51     scilab_fill_sub(Double, Bool, M_M, Double, Bool, Double);
52     scilab_fill_sub(Double, Polynom, M_M, Double, Polynom, Polynom);
53     scilab_fill_sub(Double, Sparse, M_M, Double, Sparse, Double);
54
55     //Matrix - Matrix Complex
56     scilab_fill_sub(Double, DoubleComplex, M_MC, Double, Double, Double);
57     scilab_fill_sub(Double, PolynomComplex, M_M, Double, Polynom, Polynom);
58     scilab_fill_sub(Double, SparseComplex, M_M, Double, Sparse, Double);
59
60     //Matrix - Scalar
61     scilab_fill_sub(Double, ScalarDouble, M_S, Double, Double, Double);
62     scilab_fill_sub(Double, ScalarInt8, M_S, Double, Int8, Int8);
63     scilab_fill_sub(Double, ScalarUInt8, M_S, Double, UInt8, UInt8);
64     scilab_fill_sub(Double, ScalarInt16, M_S, Double, Int16, Int16);
65     scilab_fill_sub(Double, ScalarUInt16, M_S, Double, UInt16, UInt16);
66     scilab_fill_sub(Double, ScalarInt32, M_S, Double, Int32, Int32);
67     scilab_fill_sub(Double, ScalarUInt32, M_S, Double, UInt32, UInt32);
68     scilab_fill_sub(Double, ScalarInt64, M_S, Double, Int64, Int64);
69     scilab_fill_sub(Double, ScalarUInt64, M_S, Double, UInt64, UInt64);
70     scilab_fill_sub(Double, ScalarBool, M_S, Double, Bool, Double);
71     scilab_fill_sub(Double, ScalarPolynom, M_M, Double, Polynom, Polynom);
72
73     //Matrix - Scalar Complex
74     scilab_fill_sub(Double, ScalarDoubleComplex, M_SC, Double, Double, Double);
75     scilab_fill_sub(Double, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
76     //Matrix - Empty
77     scilab_fill_sub(Double, Empty, M_E, Double, Double, Double);
78
79
80     //Matrix Complex - Matrix
81     scilab_fill_sub(DoubleComplex, Double, MC_M, Double, Double, Double);
82     scilab_fill_sub(DoubleComplex, DoubleComplex, MC_MC, Double, Double, Double);
83     scilab_fill_sub(DoubleComplex, ScalarDouble, MC_S, Double, Double, Double);
84     scilab_fill_sub(DoubleComplex, ScalarDoubleComplex, MC_SC, Double, Double, Double);
85     scilab_fill_sub(DoubleComplex, Empty, MC_E, Double, Double, Double);
86     scilab_fill_sub(DoubleComplex, Polynom, M_M, Double, Polynom, Polynom);
87     scilab_fill_sub(DoubleComplex, PolynomComplex, M_M, Double, Polynom, Polynom);
88     scilab_fill_sub(DoubleComplex, ScalarPolynom, M_M, Double, Polynom, Polynom);
89     scilab_fill_sub(DoubleComplex, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
90     scilab_fill_sub(DoubleComplex, Sparse, M_M, Double, Sparse, Double);
91     scilab_fill_sub(DoubleComplex, SparseComplex, M_M, Double, Sparse, Double);
92
93     //Scalar - Matrix
94     scilab_fill_sub(ScalarDouble, Double, S_M, Double, Double, Double);
95     scilab_fill_sub(ScalarDouble, Int8, S_M, Double, Int8, Int8);
96     scilab_fill_sub(ScalarDouble, UInt8, S_M, Double, UInt8, UInt8);
97     scilab_fill_sub(ScalarDouble, Int16, S_M, Double, Int16, Int16);
98     scilab_fill_sub(ScalarDouble, UInt16, S_M, Double, UInt16, UInt16);
99     scilab_fill_sub(ScalarDouble, Int32, S_M, Double, Int32, Int32);
100     scilab_fill_sub(ScalarDouble, UInt32, S_M, Double, UInt32, UInt32);
101     scilab_fill_sub(ScalarDouble, Int64, S_M, Double, Int64, Int64);
102     scilab_fill_sub(ScalarDouble, UInt64, S_M, Double, UInt64, UInt64);
103     scilab_fill_sub(ScalarDouble, Bool, S_M, Double, Bool, Double);
104     scilab_fill_sub(ScalarDouble, Polynom, M_M, Double, Polynom, Polynom);
105     scilab_fill_sub(ScalarDouble, Sparse, M_M, Double, Sparse, Double);
106
107     //Scalar - Matrix Complex
108     scilab_fill_sub(ScalarDouble, DoubleComplex, S_MC, Double, Double, Double);
109     scilab_fill_sub(ScalarDouble, PolynomComplex, M_M, Double, Polynom, Polynom);
110     scilab_fill_sub(ScalarDouble, SparseComplex, M_M, Double, Sparse, Double);
111
112     //Scalar - Scalar
113     scilab_fill_sub(ScalarDouble, ScalarDouble, S_S, Double, Double, Double);
114     scilab_fill_sub(ScalarDouble, ScalarInt8, S_S, Double, Int8, Int8);
115     scilab_fill_sub(ScalarDouble, ScalarUInt8, S_S, Double, UInt8, UInt8);
116     scilab_fill_sub(ScalarDouble, ScalarInt16, S_S, Double, Int16, Int16);
117     scilab_fill_sub(ScalarDouble, ScalarUInt16, S_S, Double, UInt16, UInt16);
118     scilab_fill_sub(ScalarDouble, ScalarInt32, S_S, Double, Int32, Int32);
119     scilab_fill_sub(ScalarDouble, ScalarUInt32, S_S, Double, UInt32, UInt32);
120     scilab_fill_sub(ScalarDouble, ScalarInt64, S_S, Double, Int64, Int64);
121     scilab_fill_sub(ScalarDouble, ScalarUInt64, S_S, Double, UInt64, UInt64);
122     scilab_fill_sub(ScalarDouble, ScalarBool, S_S, Double, Bool, Double);
123     scilab_fill_sub(ScalarDouble, ScalarPolynom, M_M, Double, Polynom, Polynom);
124
125     //Scalar - Scalar Complex
126     scilab_fill_sub(ScalarDouble, ScalarDoubleComplex, S_SC, Double, Double, Double);
127     scilab_fill_sub(ScalarDouble, PolynomComplex, M_M, Double, Polynom, Polynom);
128     scilab_fill_sub(ScalarDouble, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
129
130     //Scalar - Empty
131     scilab_fill_sub(ScalarDouble, Empty, S_E, Double, Double, Double);
132
133     //Scalar Complex - Matrix
134     scilab_fill_sub(ScalarDoubleComplex, Double, SC_M, Double, Double, Double);
135     scilab_fill_sub(ScalarDoubleComplex, Polynom, M_M, Double, Polynom, Polynom);
136     scilab_fill_sub(ScalarDoubleComplex, Sparse, M_M, Double, Sparse, Double);
137     //Scalar Complex - Matrix Complex
138     scilab_fill_sub(ScalarDoubleComplex, DoubleComplex, SC_MC, Double, Double, Double);
139     scilab_fill_sub(ScalarDoubleComplex, PolynomComplex, M_M, Double, Polynom, Polynom);
140     scilab_fill_sub(ScalarDoubleComplex, SparseComplex, M_M, Double, Sparse, Double);
141     //Scalar Complex - Scalar
142     scilab_fill_sub(ScalarDoubleComplex, ScalarDouble, SC_S, Double, Double, Double);
143     scilab_fill_sub(ScalarDoubleComplex, ScalarPolynom, M_M, Double, Polynom, Polynom);
144     //Scalar Complex - Scalar Complex
145     scilab_fill_sub(ScalarDoubleComplex, ScalarDoubleComplex, SC_SC, Double, Double, Double);
146     scilab_fill_sub(ScalarDoubleComplex, ScalarPolynomComplex, M_M, Double, Polynom, Polynom);
147     //Scalar Complex - Empty
148     scilab_fill_sub(ScalarDoubleComplex, Empty, SC_E, Double, Double, Double);
149
150     //Empty - Matrix
151     scilab_fill_sub(Empty, Double, E_M, Double, Double, Double);
152     scilab_fill_sub(Empty, Int8, E_M, Double, Int8, Int8);
153     scilab_fill_sub(Empty, UInt8, E_M, Double, UInt8, UInt8);
154     scilab_fill_sub(Empty, Int16, E_M, Double, Int16, Int16);
155     scilab_fill_sub(Empty, UInt16, E_M, Double, UInt16, UInt16);
156     scilab_fill_sub(Empty, Int32, E_M, Double, Int32, Int32);
157     scilab_fill_sub(Empty, UInt32, E_M, Double, UInt32, UInt32);
158     scilab_fill_sub(Empty, Int64, E_M, Double, Int64, Int64);
159     scilab_fill_sub(Empty, UInt64, E_M, Double, UInt64, UInt64);
160     //scilab_fill_sub(Empty, Bool, E_M, Double, Bool, Double);
161     scilab_fill_sub(Empty, Polynom, E_M, Double, Polynom, Polynom);
162     scilab_fill_sub(Empty, PolynomComplex, E_MC, Double, Polynom, Polynom);
163     scilab_fill_sub(Empty, Sparse, E_M, Double, Sparse, Sparse);
164     scilab_fill_sub(Empty, SparseComplex, E_MC, Double, Sparse, Sparse);
165
166     //Empty - Matrix Complex
167     scilab_fill_sub(Empty, DoubleComplex, E_MC, Double, Double, Double);
168     //Empty - Scalar
169     scilab_fill_sub(Empty, ScalarDouble, E_M, Double, Double, Double);
170     scilab_fill_sub(Empty, ScalarInt8, E_M, Double, Int8, Int8);
171     scilab_fill_sub(Empty, ScalarUInt8, E_M, Double, UInt8, UInt8);
172     scilab_fill_sub(Empty, ScalarInt16, E_M, Double, Int16, Int16);
173     scilab_fill_sub(Empty, ScalarUInt16, E_M, Double, UInt16, UInt16);
174     scilab_fill_sub(Empty, ScalarInt32, E_M, Double, Int32, Int32);
175     scilab_fill_sub(Empty, ScalarUInt32, E_M, Double, UInt32, UInt32);
176     scilab_fill_sub(Empty, ScalarInt64, E_M, Double, Int64, Int64);
177     scilab_fill_sub(Empty, ScalarUInt64, E_M, Double, UInt64, UInt64);
178     scilab_fill_sub(Empty, ScalarBool, E_M, Double, Bool, Double);
179     scilab_fill_sub(Empty, ScalarPolynom, E_M, Double, Polynom, Polynom);
180
181     //Empty - Scalar Complex
182     scilab_fill_sub(Empty, ScalarDoubleComplex, E_MC, Double, Double, Double);
183     scilab_fill_sub(Empty, ScalarPolynomComplex, E_MC, Double, Polynom, Polynom);
184     //Empty - Empty
185     scilab_fill_sub(Empty, Empty, E_M, Double, Double, Double);
186     //Empty - eye
187     scilab_fill_sub(Empty, Identity, E_I, Double, Double, Double);
188     scilab_fill_sub(Empty, IdentityComplex, E_IC, Double, Double, Double);
189
190     //Matrix - Identity
191     scilab_fill_sub(Double, Identity, M_I, Double, Double, Double);
192     scilab_fill_sub(Double, IdentityComplex, M_IC, Double, Double, Double);
193     scilab_fill_sub(DoubleComplex, Identity, MC_I, Double, Double, Double);
194     scilab_fill_sub(DoubleComplex, IdentityComplex, MC_IC, Double, Double, Double);
195     scilab_fill_sub(ScalarDouble, Identity, S_I, Double, Double, Double);
196     scilab_fill_sub(ScalarDouble, IdentityComplex, S_IC, Double, Double, Double);
197     scilab_fill_sub(ScalarDoubleComplex, Identity, SC_I, Double, Double, Double);
198     scilab_fill_sub(ScalarDoubleComplex, IdentityComplex, SC_IC, Double, Double, Double);
199
200     //Int8
201     //Matrix - Matrix
202     scilab_fill_sub(Int8, Double, M_M, Int8, Double, Int8);
203     scilab_fill_sub(Int8, Int8, M_M, Int8, Int8, Int8);
204     scilab_fill_sub(Int8, UInt8, M_M, Int8, UInt8, UInt8);
205     scilab_fill_sub(Int8, Int16, M_M, Int8, Int16, Int16);
206     scilab_fill_sub(Int8, UInt16, M_M, Int8, UInt16, UInt16);
207     scilab_fill_sub(Int8, Int32, M_M, Int8, Int32, Int32);
208     scilab_fill_sub(Int8, UInt32, M_M, Int8, UInt32, UInt32);
209     scilab_fill_sub(Int8, Int64, M_M, Int8, Int64, Int64);
210     scilab_fill_sub(Int8, UInt64, M_M, Int8, UInt64, UInt64);
211     scilab_fill_sub(Int8, Bool, M_M, Int8, Bool, Int8);
212     scilab_fill_sub(Int8, Empty, M_E, Int8, Double, Int8);
213
214     //Matrix - Scalar
215     scilab_fill_sub(Int8, ScalarDouble, M_S, Int8, Double, Int8);
216     scilab_fill_sub(Int8, ScalarInt8, M_S, Int8, Int8, Int8);
217     scilab_fill_sub(Int8, ScalarUInt8, M_S, Int8, UInt8, UInt8);
218     scilab_fill_sub(Int8, ScalarInt16, M_S, Int8, Int16, Int16);
219     scilab_fill_sub(Int8, ScalarUInt16, M_S, Int8, UInt16, UInt16);
220     scilab_fill_sub(Int8, ScalarInt32, M_S, Int8, Int32, Int32);
221     scilab_fill_sub(Int8, ScalarUInt32, M_S, Int8, UInt32, UInt32);
222     scilab_fill_sub(Int8, ScalarInt64, M_S, Int8, Int64, Int64);
223     scilab_fill_sub(Int8, ScalarUInt64, M_S, Int8, UInt64, UInt64);
224     scilab_fill_sub(Int8, ScalarBool, M_S, Int8, Bool, Int8);
225
226     //Scalar - Matrix
227     scilab_fill_sub(ScalarInt8, Double, S_M, Int8, Double, Int8);
228     scilab_fill_sub(ScalarInt8, Int8, S_M, Int8, Int8, Int8);
229     scilab_fill_sub(ScalarInt8, UInt8, S_M, Int8, UInt8, UInt8);
230     scilab_fill_sub(ScalarInt8, Int16, S_M, Int8, Int16, Int16);
231     scilab_fill_sub(ScalarInt8, UInt16, S_M, Int8, UInt16, UInt16);
232     scilab_fill_sub(ScalarInt8, Int32, S_M, Int8, Int32, Int32);
233     scilab_fill_sub(ScalarInt8, UInt32, S_M, Int8, UInt32, UInt32);
234     scilab_fill_sub(ScalarInt8, Int64, S_M, Int8, Int64, Int64);
235     scilab_fill_sub(ScalarInt8, UInt64, S_M, Int8, UInt64, UInt64);
236     scilab_fill_sub(ScalarInt8, Bool, S_M, Int8, Bool, Int8);
237     scilab_fill_sub(ScalarInt8, Empty, S_E, Int8, Double, Int8);
238
239     //Scalar - Scalar
240     scilab_fill_sub(ScalarInt8, ScalarDouble, S_S, Int8, Double, Int8);
241     scilab_fill_sub(ScalarInt8, ScalarInt8, S_S, Int8, Int8, Int8);
242     scilab_fill_sub(ScalarInt8, ScalarUInt8, S_S, Int8, UInt8, UInt8);
243     scilab_fill_sub(ScalarInt8, ScalarInt16, S_S, Int8, Int16, Int16);
244     scilab_fill_sub(ScalarInt8, ScalarUInt16, S_S, Int8, UInt16, UInt16);
245     scilab_fill_sub(ScalarInt8, ScalarInt32, S_S, Int8, Int32, Int32);
246     scilab_fill_sub(ScalarInt8, ScalarUInt32, S_S, Int8, UInt32, UInt32);
247     scilab_fill_sub(ScalarInt8, ScalarInt64, S_S, Int8, Int64, Int64);
248     scilab_fill_sub(ScalarInt8, ScalarUInt64, S_S, Int8, UInt64, UInt64);
249     scilab_fill_sub(ScalarInt8, ScalarBool, S_S, Int8, Bool, Int8);
250
251     //UInt8
252     //Matrix - Matrix
253     scilab_fill_sub(UInt8, Double, M_M, UInt8, Double, UInt8);
254     scilab_fill_sub(UInt8, Int8, M_M, UInt8, Int8, UInt8);
255     scilab_fill_sub(UInt8, UInt8, M_M, UInt8, UInt8, UInt8);
256     scilab_fill_sub(UInt8, Int16, M_M, UInt8, Int16, UInt16);
257     scilab_fill_sub(UInt8, UInt16, M_M, UInt8, UInt16, UInt16);
258     scilab_fill_sub(UInt8, Int32, M_M, UInt8, Int32, UInt32);
259     scilab_fill_sub(UInt8, UInt32, M_M, UInt8, UInt32, UInt32);
260     scilab_fill_sub(UInt8, Int64, M_M, UInt8, Int64, UInt64);
261     scilab_fill_sub(UInt8, UInt64, M_M, UInt8, UInt64, UInt64);
262     scilab_fill_sub(UInt8, Bool, M_M, UInt8, Bool, UInt8);
263     scilab_fill_sub(UInt8, Empty, M_E, UInt8, Double, UInt8);
264
265     //Matrix - Scalar
266     scilab_fill_sub(UInt8, ScalarDouble, M_S, UInt8, Double, UInt8);
267     scilab_fill_sub(UInt8, ScalarInt8, M_S, UInt8, Int8, UInt8);
268     scilab_fill_sub(UInt8, ScalarUInt8, M_S, UInt8, UInt8, UInt8);
269     scilab_fill_sub(UInt8, ScalarInt16, M_S, UInt8, Int16, UInt16);
270     scilab_fill_sub(UInt8, ScalarUInt16, M_S, UInt8, UInt16, UInt16);
271     scilab_fill_sub(UInt8, ScalarInt32, M_S, UInt8, Int32, UInt32);
272     scilab_fill_sub(UInt8, ScalarUInt32, M_S, UInt8, UInt32, UInt32);
273     scilab_fill_sub(UInt8, ScalarInt64, M_S, UInt8, Int64, UInt64);
274     scilab_fill_sub(UInt8, ScalarUInt64, M_S, UInt8, UInt64, UInt64);
275     scilab_fill_sub(UInt8, ScalarBool, M_S, UInt8, Bool, UInt8);
276
277     //Scalar - Matrix
278     scilab_fill_sub(ScalarUInt8, Double, S_M, UInt8, Double, UInt8);
279     scilab_fill_sub(ScalarUInt8, Int8, S_M, UInt8, Int8, UInt8);
280     scilab_fill_sub(ScalarUInt8, UInt8, S_M, UInt8, UInt8, UInt8);
281     scilab_fill_sub(ScalarUInt8, Int16, S_M, UInt8, Int16, UInt16);
282     scilab_fill_sub(ScalarUInt8, UInt16, S_M, UInt8, UInt16, UInt16);
283     scilab_fill_sub(ScalarUInt8, Int32, S_M, UInt8, Int32, UInt32);
284     scilab_fill_sub(ScalarUInt8, UInt32, S_M, UInt8, UInt32, UInt32);
285     scilab_fill_sub(ScalarUInt8, Int64, S_M, UInt8, Int64, UInt64);
286     scilab_fill_sub(ScalarUInt8, UInt64, S_M, UInt8, UInt64, UInt64);
287     scilab_fill_sub(ScalarUInt8, Bool, S_M, UInt8, Bool, UInt8);
288     scilab_fill_sub(ScalarUInt8, Empty, S_E, UInt8, Double, UInt8);
289
290     //Scalar - Scalar
291     scilab_fill_sub(ScalarUInt8, ScalarDouble, S_S, UInt8, Double, UInt8);
292     scilab_fill_sub(ScalarUInt8, ScalarInt8, S_S, UInt8, Int8, UInt8);
293     scilab_fill_sub(ScalarUInt8, ScalarUInt8, S_S, UInt8, UInt8, UInt8);
294     scilab_fill_sub(ScalarUInt8, ScalarInt16, S_S, UInt8, Int16, UInt16);
295     scilab_fill_sub(ScalarUInt8, ScalarUInt16, S_S, UInt8, UInt16, UInt16);
296     scilab_fill_sub(ScalarUInt8, ScalarInt32, S_S, UInt8, Int32, UInt32);
297     scilab_fill_sub(ScalarUInt8, ScalarUInt32, S_S, UInt8, UInt32, UInt32);
298     scilab_fill_sub(ScalarUInt8, ScalarInt64, S_S, UInt8, Int64, UInt64);
299     scilab_fill_sub(ScalarUInt8, ScalarUInt64, S_S, UInt8, UInt64, UInt64);
300     scilab_fill_sub(ScalarUInt8, ScalarBool, S_S, UInt8, Bool, UInt8);
301
302     //Int16
303     //Matrix - Matrix
304     scilab_fill_sub(Int16, Double, M_M, Int16, Double, Int16);
305     scilab_fill_sub(Int16, Int8, M_M, Int16, Int8, Int16);
306     scilab_fill_sub(Int16, UInt8, M_M, Int16, UInt8, UInt16);
307     scilab_fill_sub(Int16, Int16, M_M, Int16, Int16, Int16);
308     scilab_fill_sub(Int16, UInt16, M_M, Int16, UInt16, UInt16);
309     scilab_fill_sub(Int16, Int32, M_M, Int16, Int32, Int32);
310     scilab_fill_sub(Int16, UInt32, M_M, Int16, UInt32, UInt32);
311     scilab_fill_sub(Int16, Int64, M_M, Int16, Int64, Int64);
312     scilab_fill_sub(Int16, UInt64, M_M, Int16, UInt64, UInt64);
313     scilab_fill_sub(Int16, Bool, M_M, Int16, Bool, Int16);
314     scilab_fill_sub(Int16, Empty, M_E, Int16, Double, Int16);
315
316     //Matrix - Scalar
317     scilab_fill_sub(Int16, ScalarDouble, M_S, Int16, Double, Int16);
318     scilab_fill_sub(Int16, ScalarInt8, M_S, Int16, Int8, Int16);
319     scilab_fill_sub(Int16, ScalarUInt8, M_S, Int16, UInt8, UInt16);
320     scilab_fill_sub(Int16, ScalarInt16, M_S, Int16, Int16, Int16);
321     scilab_fill_sub(Int16, ScalarUInt16, M_S, Int16, UInt16, UInt16);
322     scilab_fill_sub(Int16, ScalarInt32, M_S, Int16, Int32, Int32);
323     scilab_fill_sub(Int16, ScalarUInt32, M_S, Int16, UInt32, UInt32);
324     scilab_fill_sub(Int16, ScalarInt64, M_S, Int16, Int64, Int64);
325     scilab_fill_sub(Int16, ScalarUInt64, M_S, Int16, UInt64, UInt64);
326     scilab_fill_sub(Int16, ScalarBool, M_S, Int16, Bool, Int16);
327
328     //Scalar - Matrix
329     scilab_fill_sub(ScalarInt16, Double, S_M, Int16, Double, Int16);
330     scilab_fill_sub(ScalarInt16, Int8, S_M, Int16, Int8, Int16);
331     scilab_fill_sub(ScalarInt16, UInt8, S_M, Int16, UInt8, UInt16);
332     scilab_fill_sub(ScalarInt16, Int16, S_M, Int16, Int16, Int16);
333     scilab_fill_sub(ScalarInt16, UInt16, S_M, Int16, UInt16, UInt16);
334     scilab_fill_sub(ScalarInt16, Int32, S_M, Int16, Int32, Int32);
335     scilab_fill_sub(ScalarInt16, UInt32, S_M, Int16, UInt32, UInt32);
336     scilab_fill_sub(ScalarInt16, Int64, S_M, Int16, Int64, Int64);
337     scilab_fill_sub(ScalarInt16, UInt64, S_M, Int16, UInt64, UInt64);
338     scilab_fill_sub(ScalarInt16, Bool, S_M, Int16, Bool, Int16);
339     scilab_fill_sub(ScalarInt16, Empty, S_E, Int16, Double, Int16);
340
341     //Scalar - Scalar
342     scilab_fill_sub(ScalarInt16, ScalarDouble, S_S, Int16, Double, Int16);
343     scilab_fill_sub(ScalarInt16, ScalarInt8, S_S, Int16, Int8, Int16);
344     scilab_fill_sub(ScalarInt16, ScalarUInt8, S_S, Int16, UInt8, UInt16);
345     scilab_fill_sub(ScalarInt16, ScalarInt16, S_S, Int16, Int16, Int16);
346     scilab_fill_sub(ScalarInt16, ScalarUInt16, S_S, Int16, UInt16, UInt16);
347     scilab_fill_sub(ScalarInt16, ScalarInt32, S_S, Int16, Int32, Int32);
348     scilab_fill_sub(ScalarInt16, ScalarUInt32, S_S, Int16, UInt32, UInt32);
349     scilab_fill_sub(ScalarInt16, ScalarInt64, S_S, Int16, Int64, Int64);
350     scilab_fill_sub(ScalarInt16, ScalarUInt64, S_S, Int16, UInt64, UInt64);
351     scilab_fill_sub(ScalarInt16, ScalarBool, S_S, Int16, Bool, Int16);
352
353     //UInt16
354     //Matrix - Matrix
355     scilab_fill_sub(UInt16, Double, M_M, UInt16, Double, UInt16);
356     scilab_fill_sub(UInt16, Int8, M_M, UInt16, Int8, UInt16);
357     scilab_fill_sub(UInt16, UInt8, M_M, UInt16, UInt8, UInt16);
358     scilab_fill_sub(UInt16, Int16, M_M, UInt16, Int16, UInt16);
359     scilab_fill_sub(UInt16, UInt16, M_M, UInt16, UInt16, UInt16);
360     scilab_fill_sub(UInt16, Int32, M_M, UInt16, Int32, UInt32);
361     scilab_fill_sub(UInt16, UInt32, M_M, UInt16, UInt32, UInt32);
362     scilab_fill_sub(UInt16, Int64, M_M, UInt16, Int64, UInt64);
363     scilab_fill_sub(UInt16, UInt64, M_M, UInt16, UInt64, UInt64);
364     scilab_fill_sub(UInt16, Bool, M_M, UInt16, Bool, UInt16);
365     scilab_fill_sub(UInt16, Empty, M_E, UInt16, Double, UInt16);
366
367     //Matrix - Scalar
368     scilab_fill_sub(UInt16, ScalarDouble, M_S, UInt16, Double, UInt16);
369     scilab_fill_sub(UInt16, ScalarInt8, M_S, UInt16, Int8, UInt16);
370     scilab_fill_sub(UInt16, ScalarUInt8, M_S, UInt16, UInt8, UInt16);
371     scilab_fill_sub(UInt16, ScalarInt16, M_S, UInt16, Int16, UInt16);
372     scilab_fill_sub(UInt16, ScalarUInt16, M_S, UInt16, UInt16, UInt16);
373     scilab_fill_sub(UInt16, ScalarInt32, M_S, UInt16, Int32, UInt32);
374     scilab_fill_sub(UInt16, ScalarUInt32, M_S, UInt16, UInt32, UInt32);
375     scilab_fill_sub(UInt16, ScalarInt64, M_S, UInt16, Int64, UInt64);
376     scilab_fill_sub(UInt16, ScalarUInt64, M_S, UInt16, UInt64, UInt64);
377     scilab_fill_sub(UInt16, ScalarBool, M_S, UInt16, Bool, UInt16);
378
379     //Scalar - Matrix
380     scilab_fill_sub(ScalarUInt16, Double, S_M, UInt16, Double, UInt16);
381     scilab_fill_sub(ScalarUInt16, Int8, S_M, UInt16, Int8, UInt16);
382     scilab_fill_sub(ScalarUInt16, UInt8, S_M, UInt16, UInt8, UInt16);
383     scilab_fill_sub(ScalarUInt16, Int16, S_M, UInt16, Int16, UInt16);
384     scilab_fill_sub(ScalarUInt16, UInt16, S_M, UInt16, UInt16, UInt16);
385     scilab_fill_sub(ScalarUInt16, Int32, S_M, UInt16, Int32, UInt32);
386     scilab_fill_sub(ScalarUInt16, UInt32, S_M, UInt16, UInt32, UInt32);
387     scilab_fill_sub(ScalarUInt16, Int64, S_M, UInt16, Int64, UInt64);
388     scilab_fill_sub(ScalarUInt16, UInt64, S_M, UInt16, UInt64, UInt64);
389     scilab_fill_sub(ScalarUInt16, Bool, S_M, UInt16, Bool, UInt16);
390     scilab_fill_sub(ScalarUInt16, Empty, S_E, UInt16, Double, UInt16);
391
392     //Scalar - Scalar
393     scilab_fill_sub(ScalarUInt16, ScalarDouble, S_S, UInt16, Double, UInt16);
394     scilab_fill_sub(ScalarUInt16, ScalarInt8, S_S, UInt16, Int8, UInt16);
395     scilab_fill_sub(ScalarUInt16, ScalarUInt8, S_S, UInt16, UInt8, UInt16);
396     scilab_fill_sub(ScalarUInt16, ScalarInt16, S_S, UInt16, Int16, UInt16);
397     scilab_fill_sub(ScalarUInt16, ScalarUInt16, S_S, UInt16, UInt16, UInt16);
398     scilab_fill_sub(ScalarUInt16, ScalarInt32, S_S, UInt16, Int32, UInt32);
399     scilab_fill_sub(ScalarUInt16, ScalarUInt32, S_S, UInt16, UInt32, UInt32);
400     scilab_fill_sub(ScalarUInt16, ScalarInt64, S_S, UInt16, Int64, UInt64);
401     scilab_fill_sub(ScalarUInt16, ScalarUInt64, S_S, UInt16, UInt64, UInt64);
402     scilab_fill_sub(ScalarUInt16, ScalarBool, S_S, UInt16, Bool, UInt16);
403
404     //Int32
405     //Matrix - Matrix
406     scilab_fill_sub(Int32, Double, M_M, Int32, Double, Int32);
407     scilab_fill_sub(Int32, Int8, M_M, Int32, Int8, Int32);
408     scilab_fill_sub(Int32, UInt8, M_M, Int32, UInt8, UInt32);
409     scilab_fill_sub(Int32, Int16, M_M, Int32, Int16, Int32);
410     scilab_fill_sub(Int32, UInt16, M_M, Int32, UInt16, UInt32);
411     scilab_fill_sub(Int32, Int32, M_M, Int32, Int32, Int32);
412     scilab_fill_sub(Int32, UInt32, M_M, Int32, UInt32, UInt32);
413     scilab_fill_sub(Int32, Int64, M_M, Int32, Int64, Int64);
414     scilab_fill_sub(Int32, UInt64, M_M, Int32, UInt64, UInt64);
415     scilab_fill_sub(Int32, Bool, M_M, Int32, Bool, Int32);
416     scilab_fill_sub(Int32, Empty, M_E, Int32, Double, Int32);
417
418     //Matrix - Scalar
419     scilab_fill_sub(Int32, ScalarDouble, M_S, Int32, Double, Int32);
420     scilab_fill_sub(Int32, ScalarInt8, M_S, Int32, Int8, Int32);
421     scilab_fill_sub(Int32, ScalarUInt8, M_S, Int32, UInt8, UInt32);
422     scilab_fill_sub(Int32, ScalarInt16, M_S, Int32, Int16, Int32);
423     scilab_fill_sub(Int32, ScalarUInt16, M_S, Int32, UInt16, UInt32);
424     scilab_fill_sub(Int32, ScalarInt32, M_S, Int32, Int32, Int32);
425     scilab_fill_sub(Int32, ScalarUInt32, M_S, Int32, UInt32, UInt32);
426     scilab_fill_sub(Int32, ScalarInt64, M_S, Int32, Int64, Int64);
427     scilab_fill_sub(Int32, ScalarUInt64, M_S, Int32, UInt64, UInt64);
428     scilab_fill_sub(Int32, ScalarBool, M_S, Int32, Bool, Int32);
429
430     //Scalar - Matrix
431     scilab_fill_sub(ScalarInt32, Double, S_M, Int32, Double, Int32);
432     scilab_fill_sub(ScalarInt32, Int8, S_M, Int32, Int8, Int32);
433     scilab_fill_sub(ScalarInt32, UInt8, S_M, Int32, UInt8, UInt32);
434     scilab_fill_sub(ScalarInt32, Int16, S_M, Int32, Int16, Int32);
435     scilab_fill_sub(ScalarInt32, UInt16, S_M, Int32, UInt16, UInt32);
436     scilab_fill_sub(ScalarInt32, Int32, S_M, Int32, Int32, Int32);
437     scilab_fill_sub(ScalarInt32, UInt32, S_M, Int32, UInt32, UInt32);
438     scilab_fill_sub(ScalarInt32, Int64, S_M, Int32, Int64, Int64);
439     scilab_fill_sub(ScalarInt32, UInt64, S_M, Int32, UInt64, UInt64);
440     scilab_fill_sub(ScalarInt32, Bool, S_M, Int32, Bool, Int32);
441     scilab_fill_sub(ScalarInt32, Empty, S_E, Int32, Double, Int32);
442
443     //Scalar - Scalar
444     scilab_fill_sub(ScalarInt32, ScalarDouble, S_S, Int32, Double, Int32);
445     scilab_fill_sub(ScalarInt32, ScalarInt8, S_S, Int32, Int8, Int32);
446     scilab_fill_sub(ScalarInt32, ScalarUInt8, S_S, Int32, UInt8, UInt32);
447     scilab_fill_sub(ScalarInt32, ScalarInt16, S_S, Int32, Int16, Int32);
448     scilab_fill_sub(ScalarInt32, ScalarUInt16, S_S, Int32, UInt16, UInt32);
449     scilab_fill_sub(ScalarInt32, ScalarInt32, S_S, Int32, Int32, Int32);
450     scilab_fill_sub(ScalarInt32, ScalarUInt32, S_S, Int32, UInt32, UInt32);
451     scilab_fill_sub(ScalarInt32, ScalarInt64, S_S, Int32, Int64, Int64);
452     scilab_fill_sub(ScalarInt32, ScalarUInt64, S_S, Int32, UInt64, UInt64);
453     scilab_fill_sub(ScalarInt32, ScalarBool, S_S, Int32, Bool, Int32);
454
455     //UInt32
456     //Matrix - Matrix
457     scilab_fill_sub(UInt32, Double, M_M, UInt32, Double, UInt32);
458     scilab_fill_sub(UInt32, Int8, M_M, UInt32, Int8, UInt32);
459     scilab_fill_sub(UInt32, UInt8, M_M, UInt32, UInt8, UInt32);
460     scilab_fill_sub(UInt32, Int16, M_M, UInt32, Int16, UInt32);
461     scilab_fill_sub(UInt32, UInt16, M_M, UInt32, UInt16, UInt32);
462     scilab_fill_sub(UInt32, Int32, M_M, UInt32, Int32, UInt32);
463     scilab_fill_sub(UInt32, UInt32, M_M, UInt32, UInt32, UInt32);
464     scilab_fill_sub(UInt32, Int64, M_M, UInt32, Int64, UInt64);
465     scilab_fill_sub(UInt32, UInt64, M_M, UInt32, UInt64, UInt64);
466     scilab_fill_sub(UInt32, Bool, M_M, UInt32, Bool, UInt32);
467     scilab_fill_sub(UInt32, Empty, M_E, UInt32, Double, UInt32);
468
469     //Matrix - Scalar
470     scilab_fill_sub(UInt32, ScalarDouble, M_S, UInt32, Double, UInt32);
471     scilab_fill_sub(UInt32, ScalarInt8, M_S, UInt32, Int8, UInt32);
472     scilab_fill_sub(UInt32, ScalarUInt8, M_S, UInt32, UInt8, UInt32);
473     scilab_fill_sub(UInt32, ScalarInt16, M_S, UInt32, Int16, UInt32);
474     scilab_fill_sub(UInt32, ScalarUInt16, M_S, UInt32, UInt16, UInt32);
475     scilab_fill_sub(UInt32, ScalarInt32, M_S, UInt32, Int32, UInt32);
476     scilab_fill_sub(UInt32, ScalarUInt32, M_S, UInt32, UInt32, UInt32);
477     scilab_fill_sub(UInt32, ScalarInt64, M_S, UInt32, Int64, UInt64);
478     scilab_fill_sub(UInt32, ScalarUInt64, M_S, UInt32, UInt64, UInt64);
479     scilab_fill_sub(UInt32, ScalarBool, M_S, UInt32, Bool, UInt32);
480
481     //Scalar - Matrix
482     scilab_fill_sub(ScalarUInt32, Double, S_M, UInt32, Double, UInt32);
483     scilab_fill_sub(ScalarUInt32, Int8, S_M, UInt32, Int8, UInt32);
484     scilab_fill_sub(ScalarUInt32, UInt8, S_M, UInt32, UInt8, UInt32);
485     scilab_fill_sub(ScalarUInt32, Int16, S_M, UInt32, Int16, UInt32);
486     scilab_fill_sub(ScalarUInt32, UInt16, S_M, UInt32, UInt16, UInt32);
487     scilab_fill_sub(ScalarUInt32, Int32, S_M, UInt32, Int32, UInt32);
488     scilab_fill_sub(ScalarUInt32, UInt32, S_M, UInt32, UInt32, UInt32);
489     scilab_fill_sub(ScalarUInt32, Int64, S_M, UInt32, Int64, UInt64);
490     scilab_fill_sub(ScalarUInt32, UInt64, S_M, UInt32, UInt64, UInt64);
491     scilab_fill_sub(ScalarUInt32, Bool, S_M, UInt32, Bool, UInt32);
492     scilab_fill_sub(ScalarUInt32, Empty, S_E, UInt32, Double, UInt32);
493
494     //Scalar - Scalar
495     scilab_fill_sub(ScalarUInt32, ScalarDouble, S_S, UInt32, Double, UInt32);
496     scilab_fill_sub(ScalarUInt32, ScalarInt8, S_S, UInt32, Int8, UInt32);
497     scilab_fill_sub(ScalarUInt32, ScalarUInt8, S_S, UInt32, UInt8, UInt32);
498     scilab_fill_sub(ScalarUInt32, ScalarInt16, S_S, UInt32, Int16, UInt32);
499     scilab_fill_sub(ScalarUInt32, ScalarUInt16, S_S, UInt32, UInt16, UInt32);
500     scilab_fill_sub(ScalarUInt32, ScalarInt32, S_S, UInt32, Int32, UInt32);
501     scilab_fill_sub(ScalarUInt32, ScalarUInt32, S_S, UInt32, UInt32, UInt32);
502     scilab_fill_sub(ScalarUInt32, ScalarInt64, S_S, UInt32, Int64, UInt64);
503     scilab_fill_sub(ScalarUInt32, ScalarUInt64, S_S, UInt32, UInt64, UInt64);
504     scilab_fill_sub(ScalarUInt32, ScalarBool, S_S, UInt32, Bool, UInt32);
505
506     //Int64
507     //Matrix - Matrix
508     scilab_fill_sub(Int64, Double, M_M, Int64, Double, Int64);
509     scilab_fill_sub(Int64, Int8, M_M, Int64, Int8, Int64);
510     scilab_fill_sub(Int64, UInt8, M_M, Int64, UInt8, UInt64);
511     scilab_fill_sub(Int64, Int16, M_M, Int64, Int16, Int64);
512     scilab_fill_sub(Int64, UInt16, M_M, Int64, UInt16, UInt64);
513     scilab_fill_sub(Int64, Int32, M_M, Int64, Int32, Int64);
514     scilab_fill_sub(Int64, UInt32, M_M, Int64, UInt32, UInt64);
515     scilab_fill_sub(Int64, Int64, M_M, Int64, Int64, Int64);
516     scilab_fill_sub(Int64, UInt64, M_M, Int64, UInt64, UInt64);
517     scilab_fill_sub(Int64, Bool, M_M, Int64, Bool, Int64);
518     scilab_fill_sub(Int64, Empty, M_E, Int64, Double, Int64);
519
520     //Matrix - Scalar
521     scilab_fill_sub(Int64, ScalarDouble, M_S, Int64, Double, Int64);
522     scilab_fill_sub(Int64, ScalarInt8, M_S, Int64, Int8, Int64);
523     scilab_fill_sub(Int64, ScalarUInt8, M_S, Int64, UInt8, UInt64);
524     scilab_fill_sub(Int64, ScalarInt16, M_S, Int64, Int16, Int64);
525     scilab_fill_sub(Int64, ScalarUInt16, M_S, Int64, UInt16, UInt64);
526     scilab_fill_sub(Int64, ScalarInt32, M_S, Int64, Int32, Int64);
527     scilab_fill_sub(Int64, ScalarUInt32, M_S, Int64, UInt32, UInt64);
528     scilab_fill_sub(Int64, ScalarInt64, M_S, Int64, Int64, Int64);
529     scilab_fill_sub(Int64, ScalarUInt64, M_S, Int64, UInt64, UInt64);
530     scilab_fill_sub(Int64, ScalarBool, M_S, Int64, Bool, Int64);
531
532     //Scalar - Matrix
533     scilab_fill_sub(ScalarInt64, Double, S_M, Int64, Double, Int64);
534     scilab_fill_sub(ScalarInt64, Int8, S_M, Int64, Int8, Int64);
535     scilab_fill_sub(ScalarInt64, UInt8, S_M, Int64, UInt8, UInt64);
536     scilab_fill_sub(ScalarInt64, Int16, S_M, Int64, Int16, Int64);
537     scilab_fill_sub(ScalarInt64, UInt16, S_M, Int64, UInt16, UInt64);
538     scilab_fill_sub(ScalarInt64, Int32, S_M, Int64, Int32, Int64);
539     scilab_fill_sub(ScalarInt64, UInt32, S_M, Int64, UInt32, UInt64);
540     scilab_fill_sub(ScalarInt64, Int64, S_M, Int64, Int64, Int64);
541     scilab_fill_sub(ScalarInt64, UInt64, S_M, Int64, UInt64, UInt64);
542     scilab_fill_sub(ScalarInt64, Bool, S_M, Int64, Bool, Int64);
543     scilab_fill_sub(ScalarInt64, Empty, S_E, Int64, Double, Int64);
544
545     //Scalar - Scalar
546     scilab_fill_sub(ScalarInt64, ScalarDouble, S_S, Int64, Double, Int64);
547     scilab_fill_sub(ScalarInt64, ScalarInt8, S_S, Int64, Int8, Int64);
548     scilab_fill_sub(ScalarInt64, ScalarUInt8, S_S, Int64, UInt8, UInt64);
549     scilab_fill_sub(ScalarInt64, ScalarInt16, S_S, Int64, Int16, Int64);
550     scilab_fill_sub(ScalarInt64, ScalarUInt16, S_S, Int64, UInt16, UInt64);
551     scilab_fill_sub(ScalarInt64, ScalarInt32, S_S, Int64, Int32, Int64);
552     scilab_fill_sub(ScalarInt64, ScalarUInt32, S_S, Int64, UInt32, UInt64);
553     scilab_fill_sub(ScalarInt64, ScalarInt64, S_S, Int64, Int64, Int64);
554     scilab_fill_sub(ScalarInt64, ScalarUInt64, S_S, Int64, UInt64, UInt64);
555     scilab_fill_sub(ScalarInt64, ScalarBool, S_S, Int64, Bool, Int64);
556
557     //UInt64
558     //Matrix - Matrix
559     scilab_fill_sub(UInt64, Double, M_M, UInt64, Double, UInt64);
560     scilab_fill_sub(UInt64, Int8, M_M, UInt64, Int8, UInt64);
561     scilab_fill_sub(UInt64, UInt8, M_M, UInt64, UInt8, UInt64);
562     scilab_fill_sub(UInt64, Int16, M_M, UInt64, Int16, UInt64);
563     scilab_fill_sub(UInt64, UInt16, M_M, UInt64, UInt16, UInt64);
564     scilab_fill_sub(UInt64, Int32, M_M, UInt64, Int32, UInt64);
565     scilab_fill_sub(UInt64, UInt32, M_M, UInt64, UInt32, UInt64);
566     scilab_fill_sub(UInt64, Int64, M_M, UInt64, Int64, UInt64);
567     scilab_fill_sub(UInt64, UInt64, M_M, UInt64, UInt64, UInt64);
568     scilab_fill_sub(UInt64, Bool, M_M, UInt64, Bool, UInt64);
569     scilab_fill_sub(UInt64, Empty, M_E, UInt64, Double, UInt64);
570
571     //Matrix - Scalar
572     scilab_fill_sub(UInt64, ScalarDouble, M_S, UInt64, Double, UInt64);
573     scilab_fill_sub(UInt64, ScalarInt8, M_S, UInt64, Int8, UInt64);
574     scilab_fill_sub(UInt64, ScalarUInt8, M_S, UInt64, UInt8, UInt64);
575     scilab_fill_sub(UInt64, ScalarInt16, M_S, UInt64, Int16, UInt64);
576     scilab_fill_sub(UInt64, ScalarUInt16, M_S, UInt64, UInt16, UInt64);
577     scilab_fill_sub(UInt64, ScalarInt32, M_S, UInt64, Int32, UInt64);
578     scilab_fill_sub(UInt64, ScalarUInt32, M_S, UInt64, UInt32, UInt64);
579     scilab_fill_sub(UInt64, ScalarInt64, M_S, UInt64, Int64, UInt64);
580     scilab_fill_sub(UInt64, ScalarUInt64, M_S, UInt64, UInt64, UInt64);
581     scilab_fill_sub(UInt64, ScalarBool, M_S, UInt64, Bool, UInt64);
582
583     //Scalar - Matrix
584     scilab_fill_sub(ScalarUInt64, Double, S_M, UInt64, Double, UInt64);
585     scilab_fill_sub(ScalarUInt64, Int8, S_M, UInt64, Int8, UInt64);
586     scilab_fill_sub(ScalarUInt64, UInt8, S_M, UInt64, UInt8, UInt64);
587     scilab_fill_sub(ScalarUInt64, Int16, S_M, UInt64, Int16, UInt64);
588     scilab_fill_sub(ScalarUInt64, UInt16, S_M, UInt64, UInt16, UInt64);
589     scilab_fill_sub(ScalarUInt64, Int32, S_M, UInt64, Int32, UInt64);
590     scilab_fill_sub(ScalarUInt64, UInt32, S_M, UInt64, UInt32, UInt64);
591     scilab_fill_sub(ScalarUInt64, Int64, S_M, UInt64, Int64, UInt64);
592     scilab_fill_sub(ScalarUInt64, UInt64, S_M, UInt64, UInt64, UInt64);
593     scilab_fill_sub(ScalarUInt64, Bool, S_M, UInt64, Bool, UInt64);
594     scilab_fill_sub(ScalarUInt64, Empty, S_E, UInt64, Double, UInt64);
595
596     //Scalar - Scalar
597     scilab_fill_sub(ScalarUInt64, ScalarDouble, S_S, UInt64, Double, UInt64);
598     scilab_fill_sub(ScalarUInt64, ScalarInt8, S_S, UInt64, Int8, UInt64);
599     scilab_fill_sub(ScalarUInt64, ScalarUInt8, S_S, UInt64, UInt8, UInt64);
600     scilab_fill_sub(ScalarUInt64, ScalarInt16, S_S, UInt64, Int16, UInt64);
601     scilab_fill_sub(ScalarUInt64, ScalarUInt16, S_S, UInt64, UInt16, UInt64);
602     scilab_fill_sub(ScalarUInt64, ScalarInt32, S_S, UInt64, Int32, UInt64);
603     scilab_fill_sub(ScalarUInt64, ScalarUInt32, S_S, UInt64, UInt32, UInt64);
604     scilab_fill_sub(ScalarUInt64, ScalarInt64, S_S, UInt64, Int64, UInt64);
605     scilab_fill_sub(ScalarUInt64, ScalarUInt64, S_S, UInt64, UInt64, UInt64);
606     scilab_fill_sub(ScalarUInt64, ScalarBool, S_S, UInt64, Bool, UInt64);
607
608     //Bool
609     //Matrix - Matrix
610     scilab_fill_sub(Bool, Double, M_M, Bool, Double, Double);
611     scilab_fill_sub(Bool, Int8, M_M, Bool, Int8, Int8);
612     scilab_fill_sub(Bool, UInt8, M_M, Bool, UInt8, UInt8);
613     scilab_fill_sub(Bool, Int16, M_M, Bool, Int16, Int16);
614     scilab_fill_sub(Bool, UInt16, M_M, Bool, UInt16, UInt16);
615     scilab_fill_sub(Bool, Int32, M_M, Bool, Int32, Int32);
616     scilab_fill_sub(Bool, UInt32, M_M, Bool, UInt32, UInt32);
617     scilab_fill_sub(Bool, Int64, M_M, Bool, Int64, Int64);
618     scilab_fill_sub(Bool, UInt64, M_M, Bool, UInt64, UInt64);
619     scilab_fill_sub(Bool, Bool, M_M, Bool, Bool, Bool);
620     scilab_fill_sub(Bool, Empty, M_E, Bool, Double, Double);
621
622     //Matrix - Scalar
623     scilab_fill_sub(Bool, ScalarDouble, M_S, Bool, Double, Double);
624     scilab_fill_sub(Bool, ScalarInt8, M_S, Bool, Int8, Int8);
625     scilab_fill_sub(Bool, ScalarUInt8, M_S, Bool, UInt8, UInt8);
626     scilab_fill_sub(Bool, ScalarInt16, M_S, Bool, Int16, Int16);
627     scilab_fill_sub(Bool, ScalarUInt16, M_S, Bool, UInt16, UInt16);
628     scilab_fill_sub(Bool, ScalarInt32, M_S, Bool, Int32, Int32);
629     scilab_fill_sub(Bool, ScalarUInt32, M_S, Bool, UInt32, UInt32);
630     scilab_fill_sub(Bool, ScalarInt64, M_S, Bool, Int64, Int64);
631     scilab_fill_sub(Bool, ScalarUInt64, M_S, Bool, UInt64, UInt64);
632     scilab_fill_sub(Bool, ScalarBool, M_S, Bool, Bool, Bool);
633
634     //Scalar - Matrix
635     scilab_fill_sub(ScalarBool, Double, S_M, Bool, Double, Double);
636     scilab_fill_sub(ScalarBool, Int8, S_M, Bool, Int8, Int8);
637     scilab_fill_sub(ScalarBool, UInt8, S_M, Bool, UInt8, UInt8);
638     scilab_fill_sub(ScalarBool, Int16, S_M, Bool, Int16, Int16);
639     scilab_fill_sub(ScalarBool, UInt16, S_M, Bool, UInt16, UInt16);
640     scilab_fill_sub(ScalarBool, Int32, S_M, Bool, Int32, Int32);
641     scilab_fill_sub(ScalarBool, UInt32, S_M, Bool, UInt32, UInt32);
642     scilab_fill_sub(ScalarBool, Int64, S_M, Bool, Int64, Int64);
643     scilab_fill_sub(ScalarBool, UInt64, S_M, Bool, UInt64, UInt64);
644     scilab_fill_sub(ScalarBool, Bool, S_M, Bool, Bool, Bool);
645     scilab_fill_sub(ScalarBool, Empty, S_E, Bool, Double, Double);
646
647     //Scalar - Scalar
648     scilab_fill_sub(ScalarBool, ScalarDouble, S_S, Bool, Double, Double);
649     scilab_fill_sub(ScalarBool, ScalarInt8, S_S, Bool, Int8, Int8);
650     scilab_fill_sub(ScalarBool, ScalarUInt8, S_S, Bool, UInt8, UInt8);
651     scilab_fill_sub(ScalarBool, ScalarInt16, S_S, Bool, Int16, Int16);
652     scilab_fill_sub(ScalarBool, ScalarUInt16, S_S, Bool, UInt16, UInt16);
653     scilab_fill_sub(ScalarBool, ScalarInt32, S_S, Bool, Int32, Int32);
654     scilab_fill_sub(ScalarBool, ScalarUInt32, S_S, Bool, UInt32, UInt32);
655     scilab_fill_sub(ScalarBool, ScalarInt64, S_S, Bool, Int64, Int64);
656     scilab_fill_sub(ScalarBool, ScalarUInt64, S_S, Bool, UInt64, UInt64);
657     scilab_fill_sub(ScalarBool, ScalarBool, S_S, Bool, Bool, Bool);
658
659     //Identity
660     scilab_fill_sub(Identity, Double, I_M, Double, Double, Double);
661     scilab_fill_sub(Identity, DoubleComplex, I_MC, Double, Double, Double);
662     scilab_fill_sub(Identity, ScalarDouble, I_S, Double, Double, Double);
663     scilab_fill_sub(Identity, ScalarDoubleComplex, I_SC, Double, Double, Double);
664     scilab_fill_sub(Identity, Identity, I_I, Double, Double, Double);
665     scilab_fill_sub(Identity, IdentityComplex, I_IC, Double, Double, Double);
666     scilab_fill_sub(Identity, Empty, I_E, Double, Double, Double);
667
668     scilab_fill_sub(Identity, Polynom, I_M, Double, Polynom, Polynom);
669     scilab_fill_sub(Identity, PolynomComplex, I_MC, Double, Polynom, Polynom);
670     scilab_fill_sub(Identity, ScalarPolynom, I_M, Double, Polynom, Polynom);
671     scilab_fill_sub(Identity, ScalarPolynomComplex, I_MC, Double, Polynom, Polynom);
672     scilab_fill_sub(Identity, Sparse, M_M, Double, Sparse, Sparse);
673     scilab_fill_sub(Identity, SparseComplex, M_M, Double, Sparse, Sparse);
674
675     scilab_fill_sub(IdentityComplex, Double, IC_M, Double, Double, Double);
676     scilab_fill_sub(IdentityComplex, DoubleComplex, IC_MC, Double, Double, Double);
677     scilab_fill_sub(IdentityComplex, ScalarDouble, IC_S, Double, Double, Double);
678     scilab_fill_sub(IdentityComplex, ScalarDoubleComplex, IC_SC, Double, Double, Double);
679     scilab_fill_sub(IdentityComplex, Identity, IC_I, Double, Double, Double);
680     scilab_fill_sub(IdentityComplex, IdentityComplex, IC_IC, Double, Double, Double);
681     scilab_fill_sub(IdentityComplex, Empty, IC_E, Double, Double, Double);
682
683     scilab_fill_sub(IdentityComplex, Polynom, IC_M, Double, Polynom, Polynom);
684     scilab_fill_sub(IdentityComplex, PolynomComplex, IC_MC, Double, Polynom, Polynom);
685     scilab_fill_sub(IdentityComplex, ScalarPolynom, IC_M, Double, Polynom, Polynom);
686     scilab_fill_sub(IdentityComplex, ScalarPolynomComplex, IC_MC, Double, Polynom, Polynom);
687     scilab_fill_sub(IdentityComplex, Sparse, M_M, Double, Sparse, Sparse);
688     scilab_fill_sub(IdentityComplex, SparseComplex, M_M, Double, Sparse, Sparse);
689
690     scilab_fill_sub(Identity, ScalarInt8, I_S, Double, Int8, Int8);
691     scilab_fill_sub(Identity, ScalarUInt8, I_S, Double, UInt8, UInt8);
692     scilab_fill_sub(Identity, ScalarInt16, I_S, Double, Int16, Int16);
693     scilab_fill_sub(Identity, ScalarUInt16, I_S, Double, UInt16, UInt16);
694     scilab_fill_sub(Identity, ScalarInt32, I_S, Double, Int32, Int32);
695     scilab_fill_sub(Identity, ScalarUInt32, I_S, Double, UInt32, UInt32);
696     scilab_fill_sub(Identity, ScalarInt64, I_S, Double, Int64, Int64);
697     scilab_fill_sub(Identity, ScalarUInt64, I_S, Double, UInt64, UInt64);
698
699     scilab_fill_sub(Identity, Int8, I_M, Double, Int8, Int8);
700     scilab_fill_sub(Identity, UInt8, I_M, Double, UInt8, UInt8);
701     scilab_fill_sub(Identity, Int16, I_M, Double, Int16, Int16);
702     scilab_fill_sub(Identity, UInt16, I_M, Double, UInt16, UInt16);
703     scilab_fill_sub(Identity, Int32, I_M, Double, Int32, Int32);
704     scilab_fill_sub(Identity, UInt32, I_M, Double, UInt32, UInt32);
705     scilab_fill_sub(Identity, Int64, I_M, Double, Int64, Int64);
706     scilab_fill_sub(Identity, UInt64, I_M, Double, UInt64, UInt64);
707
708     //Polynom
709
710     //poly - poly
711     scilab_fill_sub(Polynom, Polynom, M_M, Polynom, Polynom, Polynom);
712     scilab_fill_sub(Polynom, PolynomComplex, M_M, Polynom, Polynom, Polynom);
713     scilab_fill_sub(PolynomComplex, Polynom, M_M, Polynom, Polynom, Polynom);
714     scilab_fill_sub(PolynomComplex, PolynomComplex, M_M, Polynom, Polynom, Polynom);
715
716     //poly - scalar poly
717     scilab_fill_sub(Polynom, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
718     scilab_fill_sub(Polynom, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
719     scilab_fill_sub(PolynomComplex, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
720     scilab_fill_sub(PolynomComplex, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
721
722     //poly - double
723     scilab_fill_sub(Polynom, Double, M_M, Polynom, Double, Polynom);
724     scilab_fill_sub(Polynom, DoubleComplex, M_M, Polynom, Double, Polynom);
725     scilab_fill_sub(PolynomComplex, Double, M_M, Polynom, Double, Polynom);
726     scilab_fill_sub(PolynomComplex, DoubleComplex, M_M, Polynom, Double, Polynom);
727
728     //poly - scalar double
729     scilab_fill_sub(Polynom, ScalarDouble, M_M, Polynom, Double, Polynom);
730     scilab_fill_sub(Polynom, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
731     scilab_fill_sub(PolynomComplex, ScalarDouble, M_M, Polynom, Double, Polynom);
732     scilab_fill_sub(PolynomComplex, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
733
734     //poly - []
735     scilab_fill_sub(Polynom, Empty, M_E, Polynom, Double, Polynom);
736     scilab_fill_sub(PolynomComplex, Empty, M_E, Polynom, Double, Polynom);
737
738     //poly - eye
739     scilab_fill_sub(Polynom, Identity, M_M, Polynom, Double, Polynom);
740     scilab_fill_sub(Polynom, IdentityComplex, M_M, Polynom, Double, Polynom);
741     scilab_fill_sub(PolynomComplex, Identity, M_M, Polynom, Double, Polynom);
742     scilab_fill_sub(PolynomComplex, IdentityComplex, M_M, Polynom, Double, Polynom);
743
744     //scalar poly - poly
745     scilab_fill_sub(ScalarPolynom, Polynom, M_M, Polynom, Polynom, Polynom);
746     scilab_fill_sub(ScalarPolynom, PolynomComplex, M_M, Polynom, Polynom, Polynom);
747     scilab_fill_sub(ScalarPolynomComplex, Polynom, M_M, Polynom, Polynom, Polynom);
748     scilab_fill_sub(ScalarPolynomComplex, PolynomComplex, M_M, Polynom, Polynom, Polynom);
749
750     //scalar poly - scalar poly
751     scilab_fill_sub(ScalarPolynom, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
752     scilab_fill_sub(ScalarPolynom, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
753     scilab_fill_sub(ScalarPolynomComplex, ScalarPolynom, M_M, Polynom, Polynom, Polynom);
754     scilab_fill_sub(ScalarPolynomComplex, ScalarPolynomComplex, M_M, Polynom, Polynom, Polynom);
755
756     //scalar poly - double
757     scilab_fill_sub(ScalarPolynom, Double, M_M, Polynom, Double, Polynom);
758     scilab_fill_sub(ScalarPolynom, DoubleComplex, M_M, Polynom, Double, Polynom);
759     scilab_fill_sub(ScalarPolynomComplex, Double, M_M, Polynom, Double, Polynom);
760     scilab_fill_sub(ScalarPolynomComplex, DoubleComplex, M_M, Polynom, Double, Polynom);
761
762     //scalar poly - scalar double
763     scilab_fill_sub(ScalarPolynom, ScalarDouble, M_M, Polynom, Double, Polynom);
764     scilab_fill_sub(ScalarPolynom, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
765     scilab_fill_sub(ScalarPolynomComplex, ScalarDouble, M_M, Polynom, Double, Polynom);
766     scilab_fill_sub(ScalarPolynomComplex, ScalarDoubleComplex, M_M, Polynom, Double, Polynom);
767
768     //scalar poly - []
769     scilab_fill_sub(ScalarPolynom, Empty, M_E, Polynom, Double, Polynom);
770     scilab_fill_sub(ScalarPolynomComplex, Empty, M_E, Polynom, Double, Polynom);
771
772     //scalar poly - eye
773     scilab_fill_sub(ScalarPolynom, Identity, M_M, Polynom, Double, Polynom);
774     scilab_fill_sub(ScalarPolynom, IdentityComplex, M_M, Polynom, Double, Polynom);
775     scilab_fill_sub(ScalarPolynomComplex, Identity, M_M, Polynom, Double, Polynom);
776     scilab_fill_sub(ScalarPolynomComplex, IdentityComplex, M_M, Polynom, Double, Polynom);
777
778     //Sparse
779     scilab_fill_sub(Sparse, Sparse, M_M, Sparse, Sparse, Sparse);
780     scilab_fill_sub(Sparse, SparseComplex, M_M, Sparse, Sparse, Sparse);
781     scilab_fill_sub(Sparse, Double, M_M, Sparse, Double, Double);
782     scilab_fill_sub(Sparse, DoubleComplex, M_M, Sparse, Double, Double);
783     scilab_fill_sub(Sparse, ScalarDouble, M_M, Sparse, Double, Double);
784     scilab_fill_sub(Sparse, ScalarDoubleComplex, M_M, Sparse, Double, Double);
785
786     scilab_fill_sub(Sparse, Empty, M_E, Sparse, Double, Sparse);
787     scilab_fill_sub(Sparse, Identity, M_M, Sparse, Double, Sparse);
788     scilab_fill_sub(Sparse, IdentityComplex, M_M, Sparse, Double, Sparse);
789
790     scilab_fill_sub(SparseComplex, Sparse, M_M, Sparse, Sparse, Sparse);
791     scilab_fill_sub(SparseComplex, SparseComplex, M_M, Sparse, Sparse, Sparse);
792     scilab_fill_sub(SparseComplex, Double, M_M, Sparse, Double, Double);
793     scilab_fill_sub(SparseComplex, DoubleComplex, M_M, Sparse, Double, Double);
794     scilab_fill_sub(SparseComplex, ScalarDouble, M_M, Sparse, Double, Double);
795     scilab_fill_sub(SparseComplex, ScalarDoubleComplex, M_M, Sparse, Double, Double);
796
797     scilab_fill_sub(SparseComplex, Empty, M_E, Sparse, Double, Sparse);
798     scilab_fill_sub(SparseComplex, Identity, M_M, Sparse, Double, Sparse);
799     scilab_fill_sub(SparseComplex, IdentityComplex, M_M, Sparse, Double, Sparse);
800
801 #undef scilab_fill_sub
802 }
803
804 InternalType* GenericMinus(InternalType* _pLeftOperand, InternalType* _pRightOperand)
805 {
806     InternalType *pResult = NULL;
807
808     sub_function sub = pSubfunction[_pLeftOperand->getId()][_pRightOperand->getId()];
809     if (sub)
810     {
811         pResult = sub(_pLeftOperand, _pRightOperand);
812         if (pResult)
813         {
814             return pResult;
815         }
816     }
817
818     /*
819     ** Default case : Return NULL will Call Overloading.
820     */
821     return NULL;
822 }
823
824 template<class T, class U, class O>
825 InternalType* sub_M_M(T *_pL, U *_pR)
826 {
827     int iDimsL = _pL->getDims();
828     int iDimsR = _pR->getDims();
829
830     if (iDimsL != iDimsR)
831     {
832         return nullptr;
833     }
834
835     int* piDimsL = _pL->getDimsArray();
836     int* piDimsR = _pR->getDimsArray();
837
838     for (int i = 0 ; i < iDimsL ; ++i)
839     {
840         if (piDimsL[i] != piDimsR[i])
841         {
842             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
843         }
844     }
845
846     O* pOut = new O(iDimsL, piDimsL);
847
848     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(), pOut->get());
849     return pOut;
850 }
851
852 template<class T, class U, class O>
853 InternalType* sub_M_MC(T *_pL, U *_pR)
854 {
855     int iDimsL = _pL->getDims();
856     int iDimsR = _pR->getDims();
857
858     if (iDimsL != iDimsR)
859     {
860         return nullptr;
861     }
862
863     int* piDimsL = _pL->getDimsArray();
864     int* piDimsR = _pR->getDimsArray();
865
866     for (int i = 0 ; i < iDimsL ; ++i)
867     {
868         if (piDimsL[i] != piDimsR[i])
869         {
870             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
871         }
872     }
873
874     O* pOut = new O(iDimsL, piDimsL, true);
875
876     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
877     return pOut;
878 }
879
880 template<class T, class U, class O>
881 InternalType* sub_M_S(T *_pL, U *_pR)
882 {
883     O* pOut = new O(_pL->getDims(), _pL->getDimsArray());
884     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(0), pOut->get());
885     return pOut;
886 }
887
888 template<class T, class U, class O>
889 InternalType* sub_M_SC(T *_pL, U *_pR)
890 {
891     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
892     sub(_pL->get(), (size_t)_pL->getSize(), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
893     return pOut;
894 }
895
896 template<class T, class U, class O>
897 InternalType* sub_M_E(T *_pL, U * /*_pR*/)
898 {
899     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
900     Double* pOut = Double::Empty();
901     sub();
902     return pOut;
903 }
904
905
906 template<class T, class U, class O>
907 InternalType* sub_MC_M(T *_pL, U *_pR)
908 {
909     int iDimsL = _pL->getDims();
910     int iDimsR = _pR->getDims();
911
912     if (iDimsL != iDimsR)
913     {
914         return nullptr;
915     }
916
917     int* piDimsL = _pL->getDimsArray();
918     int* piDimsR = _pR->getDimsArray();
919
920     for (int i = 0 ; i < iDimsL ; ++i)
921     {
922         if (piDimsL[i] != piDimsR[i])
923         {
924             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
925         }
926     }
927
928     O* pOut = new O(iDimsL, piDimsL, true);
929
930     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(), pOut->get(), pOut->getImg());
931     return pOut;
932 }
933
934 template<class T, class U, class O>
935 InternalType* sub_MC_MC(T *_pL, U *_pR)
936 {
937     int iDimsL = _pL->getDims();
938     int iDimsR = _pR->getDims();
939
940     if (iDimsL != iDimsR)
941     {
942         return nullptr;
943     }
944
945     int* piDimsL = _pL->getDimsArray();
946     int* piDimsR = _pR->getDimsArray();
947
948     for (int i = 0 ; i < iDimsL ; ++i)
949     {
950         if (piDimsL[i] != piDimsR[i])
951         {
952             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
953         }
954     }
955
956     O* pOut = new O(iDimsL, piDimsL, true);
957
958     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
959     return pOut;
960 }
961
962 template<class T, class U, class O>
963 InternalType* sub_MC_S(T *_pL, U *_pR)
964 {
965     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
966     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(0), pOut->get(), pOut->getImg());
967     return pOut;
968 }
969
970 template<class T, class U, class O>
971 InternalType* sub_MC_SC(T *_pL, U *_pR)
972 {
973     O* pOut = new O(_pL->getDims(), _pL->getDimsArray(), true);
974     sub(_pL->get(), _pL->getImg(), (size_t)_pL->getSize(), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
975     return pOut;
976 }
977
978 template<class T, class U, class O>
979 InternalType* sub_MC_E(T *_pL, U * /*_pR*/)
980 {
981     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
982     Double* pOut = Double::Empty();
983     sub();
984     return pOut;
985 }
986
987
988 template<class T, class U, class O>
989 InternalType* sub_S_M(T *_pL, U *_pR)
990 {
991     O* pOut = new O(_pR->getDims(), _pR->getDimsArray());
992     sub(_pL->get(0), (size_t)_pR->getSize(), _pR->get(), pOut->get());
993     return pOut;
994 }
995
996 template<class T, class U, class O>
997 InternalType* sub_S_MC(T *_pL, U *_pR)
998 {
999     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1000     sub(_pL->get(0), (size_t)_pR->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
1001     return pOut;
1002 }
1003
1004 template<class T, class U, class O>
1005 InternalType* sub_S_S(T *_pL, U *_pR)
1006 {
1007     O* pOut = new O(0);
1008     sub(_pL->get(0), _pR->get(0), pOut->get());
1009     return pOut;
1010 }
1011
1012 template<class T, class U, class O>
1013 InternalType* sub_S_SC(T *_pL, U *_pR)
1014 {
1015     O* pOut = new O(0.0, 0.0);
1016     sub(_pL->get(), (size_t)1, _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1017     return pOut;
1018 }
1019
1020 template<class T, class U, class O>
1021 InternalType* sub_S_E(T *_pL, U * /*_pR*/)
1022 {
1023     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1024     Double* pOut = Double::Empty();
1025     sub();
1026     return pOut;
1027 }
1028
1029
1030 template<class T, class U, class O>
1031 InternalType* sub_SC_M(T *_pL, U *_pR)
1032 {
1033     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1034     sub(_pL->get(0), _pL->getImg(0), (size_t)_pR->getSize(), _pR->get(), pOut->get(), pOut->getImg());
1035     return pOut;
1036 }
1037
1038 template<class T, class U, class O>
1039 InternalType* sub_SC_MC(T *_pL, U *_pR)
1040 {
1041     O* pOut = new O(_pR->getDims(), _pR->getDimsArray(), true);
1042     sub(_pL->get(0), _pL->getImg(0), (size_t)_pR->getSize(), _pR->get(), _pR->getImg(), pOut->get(), pOut->getImg());
1043     return pOut;
1044 }
1045
1046 template<class T, class U, class O>
1047 InternalType* sub_SC_S(T *_pL, U *_pR)
1048 {
1049     O* pOut = new O(0.0, 0.0);
1050     sub(_pL->get(), _pL->getImg(), (size_t)1, _pR->get(0), pOut->get(), pOut->getImg());
1051     return pOut;
1052 }
1053
1054 template<class T, class U, class O>
1055 InternalType* sub_SC_SC(T *_pL, U *_pR)
1056 {
1057     O* pOut = new O(0.0, 0.0);
1058     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1059     return pOut;
1060 }
1061
1062 template<class T, class U, class O>
1063 InternalType* sub_SC_E(T *_pL, U * /*_pR*/)
1064 {
1065     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1066     Double* pOut = Double::Empty();
1067     sub();
1068     return pOut;
1069 }
1070
1071
1072 template<class T, class U, class O>
1073 InternalType* sub_E_M(T * /*_pL*/, U *_pR)
1074 {
1075     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1076     Double* pOut = Double::Empty();
1077     sub();
1078     return pOut;
1079 }
1080
1081 template<class T, class U, class O>
1082 InternalType* sub_E_MC(T * /*_pL*/, U *_pR)
1083 {
1084     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1085     Double* pOut = Double::Empty();
1086     sub();
1087     return pOut;
1088 }
1089
1090 template<class T, class U, class O>
1091 InternalType* sub_E_E(T * /*_pL*/, U * /*_pR*/)
1092 {
1093     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1094     Double* pOut = Double::Empty();
1095     sub();
1096     return pOut;
1097 }
1098
1099 template<class T, class U, class O>
1100 InternalType* sub_I_M(T *_pL, U *_pR)
1101 {
1102     int iDims = _pR->getDims();
1103     int* piDims = _pR->getDimsArray();
1104     O* pOut = (O*)opposite_M<U, O>(_pR);
1105     double dblLeft = _pL->get(0);
1106     int iLeadDims = piDims[0];
1107     int* piIndex = new int[iDims];
1108     piIndex[0] = 0;
1109
1110     //find smaller dims
1111     for (int i = 1 ; i < iDims ; ++i)
1112     {
1113         //init
1114         piIndex[i] = 0;
1115
1116         if (iLeadDims > piDims[i])
1117         {
1118             iLeadDims = piDims[i];
1119         }
1120     }
1121
1122     for (int i = 0 ; i < iLeadDims ; ++i)
1123     {
1124         for (int j = 0 ; j < iDims ; ++j)
1125         {
1126             piIndex[j] = i;
1127         }
1128
1129         int index = _pR->getIndex(piIndex);
1130         sub(dblLeft, _pR->get(index), pOut->get() + index);
1131     }
1132
1133     delete[] piIndex;
1134     return pOut;
1135 }
1136
1137 template<class T, class U, class O>
1138 InternalType* sub_I_MC(T *_pL, U *_pR)
1139 {
1140     int iDims = _pR->getDims();
1141     int* piDims = _pR->getDimsArray();
1142     O* pOut = (O*)opposite_MC<U, O>(_pR);
1143     double* pdblOut = pOut->get();
1144     double* pdblRight = _pR->get();
1145     double dblLeft = _pL->get(0);
1146     int iLeadDims = piDims[0];
1147     int* piIndex = new int[iDims];
1148     piIndex[0] = 0;
1149
1150     //find smaller dims
1151     for (int i = 1 ; i < iDims ; ++i)
1152     {
1153         //init
1154         piIndex[i] = 0;
1155
1156         if (iLeadDims > piDims[i])
1157         {
1158             iLeadDims = piDims[i];
1159         }
1160     }
1161
1162     for (int i = 0 ; i < iLeadDims ; ++i)
1163     {
1164         for (int j = 0 ; j < iDims ; ++j)
1165         {
1166             piIndex[j] = i;
1167         }
1168
1169         int index = _pR->getIndex(piIndex);
1170         sub(dblLeft, pdblRight[index], pdblOut + index);
1171     }
1172
1173     delete[] piIndex;
1174     return pOut;
1175 }
1176
1177 template<class T, class U, class O>
1178 InternalType* sub_IC_M(T *_pL, U *_pR)
1179 {
1180     int iDims = _pR->getDims();
1181     int* piDims = _pR->getDimsArray();
1182     O* pOut = (O*)opposite_M<U, O>(_pR);
1183     pOut->setComplex(true);
1184     double* pdblOutR = pOut->get();
1185     double* pdblOutI = pOut->getImg();
1186     double* pdblRight = _pR->get();
1187     double dblLeftR = _pL->get(0);
1188     double dblLeftI = _pL->getImg(0);
1189     int iLeadDims = piDims[0];
1190     int* piIndex = new int[iDims];
1191     piIndex[0] = 0;
1192
1193     //find smaller dims
1194     for (int i = 1 ; i < iDims ; ++i)
1195     {
1196         //init
1197         piIndex[i] = 0;
1198
1199         if (iLeadDims > piDims[i])
1200         {
1201             iLeadDims = piDims[i];
1202         }
1203     }
1204
1205     for (int i = 0 ; i < iLeadDims ; ++i)
1206     {
1207         for (int j = 0 ; j < iDims ; ++j)
1208         {
1209             piIndex[j] = i;
1210         }
1211
1212         int index = _pR->getIndex(piIndex);
1213         sub(&dblLeftR, &dblLeftI, (size_t)1, pdblRight[index], pdblOutR + index, pdblOutI + index);
1214     }
1215
1216     delete[] piIndex;
1217     return pOut;
1218 }
1219
1220 template<class T, class U, class O>
1221 InternalType* sub_IC_MC(T *_pL, U *_pR)
1222 {
1223     int iDims = _pR->getDims();
1224     int* piDims = _pR->getDimsArray();
1225     O* pOut = (O*)opposite_MC<U, O>(_pR);
1226     double dblLeftR = _pL->get(0);
1227     double dblLeftI = _pL->getImg(0);
1228     double* pdblOutR = pOut->get();
1229     double* pdblOutI = pOut->getImg();
1230     double* pdblRightR = _pR->get();
1231     double* pdblRightI = _pR->getImg();
1232     int iLeadDims = piDims[0];
1233     int* piIndex = new int[iDims];
1234     piIndex[0] = 0;
1235     //find smaller dims
1236     for (int i = 1 ; i < iDims ; ++i)
1237     {
1238         //init
1239         piIndex[i] = 0;
1240
1241         if (iLeadDims > piDims[i])
1242         {
1243             iLeadDims = piDims[i];
1244         }
1245     }
1246
1247     for (int i = 0 ; i < iLeadDims ; ++i)
1248     {
1249         for (int j = 0 ; j < iDims ; ++j)
1250         {
1251             piIndex[j] = i;
1252         }
1253
1254         int index = _pR->getIndex(piIndex);
1255         sub(dblLeftR, dblLeftI, pdblRightR[index], pdblRightI[index], pdblOutR + index, pdblOutI + index);
1256     }
1257
1258     delete[] piIndex;
1259     return pOut;
1260 }
1261
1262 template<class T, class U, class O>
1263 InternalType* sub_I_S(T *_pL, U *_pR)
1264 {
1265     return sub_S_S<T, U, O>(_pL, _pR);
1266 }
1267
1268 template<class T, class U, class O>
1269 InternalType* sub_IC_S(T *_pL, U *_pR)
1270 {
1271     return sub_SC_S<T, U, O>(_pL, _pR);
1272 }
1273
1274 template<class T, class U, class O>
1275 InternalType* sub_I_SC(T *_pL, U *_pR)
1276 {
1277     return sub_S_SC<T, U, O>(_pL, _pR);
1278 }
1279
1280 template<class T, class U, class O>
1281 InternalType* sub_IC_SC(T *_pL, U *_pR)
1282 {
1283     return sub_SC_SC<T, U, O>(_pL, _pR);
1284 }
1285
1286 template<class T, class U, class O> InternalType* sub_M_I(T *_pL, U *_pR)
1287 {
1288     int iDims = _pL->getDims();
1289     int* piDims = _pL->getDimsArray();
1290     O* pOut = (O*)_pL->clone();
1291     double* pdblOutR = pOut->get();
1292     double* pdblLeft = _pL->get();
1293     double dblRight = _pR->get(0);
1294     int iLeadDims = piDims[0];
1295     int* piIndex = new int[iDims];
1296     piIndex[0] = 0;
1297
1298     //find smaller dims
1299     for (int i = 1 ; i < iDims ; ++i)
1300     {
1301         //init
1302         piIndex[i] = 0;
1303
1304         if (iLeadDims > piDims[i])
1305         {
1306             iLeadDims = piDims[i];
1307         }
1308     }
1309
1310     for (int i = 0 ; i < iLeadDims ; ++i)
1311     {
1312         for (int j = 0 ; j < iDims ; ++j)
1313         {
1314             piIndex[j] = i;
1315         }
1316
1317         int index = _pL->getIndex(piIndex);
1318         sub(pdblLeft[index], (size_t)1, &dblRight, pdblOutR + index);
1319     }
1320
1321     delete[] piIndex;
1322     return pOut;
1323 }
1324
1325 template<class T, class U, class O> InternalType* sub_MC_I(T *_pL, U *_pR)
1326 {
1327     return sub_M_I<T, U, O>(_pL, _pR);
1328 }
1329
1330 template<class T, class U, class O> InternalType* sub_M_IC(T *_pL, U *_pR)
1331 {
1332     int iDims = _pL->getDims();
1333     int* piDims = _pL->getDimsArray();
1334     O* pOut = (O*)_pL->clone();
1335     pOut->setComplex(true);
1336     double* pdblOutR = pOut->get();
1337     double* pdblOutI = pOut->getImg();
1338     double* pdblLeft = _pL->get();
1339     double dblRightR = _pR->get(0);
1340     double dblRightI = _pR->getImg(0);
1341     int iLeadDims = piDims[0];
1342     int* piIndex = new int[iDims];
1343     piIndex[0] = 0;
1344
1345     //find smaller dims
1346     for (int i = 1 ; i < iDims ; ++i)
1347     {
1348         //init
1349         piIndex[i] = 0;
1350
1351         if (iLeadDims > piDims[i])
1352         {
1353             iLeadDims = piDims[i];
1354         }
1355     }
1356
1357     for (int i = 0 ; i < iLeadDims ; ++i)
1358     {
1359         for (int j = 0 ; j < iDims ; ++j)
1360         {
1361             piIndex[j] = i;
1362         }
1363
1364         int index = _pL->getIndex(piIndex);
1365         sub(pdblLeft[index], (size_t)1, &dblRightR, &dblRightI, pdblOutR + index, pdblOutI + index);
1366     }
1367
1368     delete[] piIndex;
1369     return pOut;
1370 }
1371
1372 template<class T, class U, class O> InternalType* sub_MC_IC(T *_pL, U *_pR)
1373 {
1374     int iDims = _pL->getDims();
1375     int* piDims = _pL->getDimsArray();
1376     O* pOut = (O*)_pL->clone();
1377     pOut->setComplex(true);
1378     double* pdblOutR = pOut->get();
1379     double* pdblOutI = pOut->getImg();
1380     double* pdblLeftR = _pL->get();
1381     double* pdblLeftI = _pL->getImg();
1382     double dblRightR = _pR->get(0);
1383     double dblRightI = _pR->getImg(0);
1384     int iLeadDims = piDims[0];
1385     int* piIndex = new int[iDims];
1386     piIndex[0] = 0;
1387
1388     //find smaller dims
1389     for (int i = 1 ; i < iDims ; ++i)
1390     {
1391         //init
1392         piIndex[i] = 0;
1393
1394         if (iLeadDims > piDims[i])
1395         {
1396             iLeadDims = piDims[i];
1397         }
1398     }
1399
1400     for (int i = 0 ; i < iLeadDims ; ++i)
1401     {
1402         for (int j = 0 ; j < iDims ; ++j)
1403         {
1404             piIndex[j] = i;
1405         }
1406
1407         int index = _pL->getIndex(piIndex);
1408         sub(pdblLeftR[index], pdblLeftI[index], (size_t)1, &dblRightR, &dblRightI, pdblOutR + index, pdblOutI + index);
1409     }
1410
1411     delete[] piIndex;
1412     return pOut;
1413 }
1414
1415 template<class T, class U, class O> InternalType* sub_S_I(T *_pL, U *_pR)
1416 {
1417     return sub_S_S<T, U, O>(_pL, _pR);
1418 }
1419
1420 template<class T, class U, class O> InternalType* sub_SC_I(T *_pL, U *_pR)
1421 {
1422     return sub_SC_SC<T, U, O>(_pL, _pR);
1423 }
1424
1425 template<class T, class U, class O> InternalType* sub_S_IC(T *_pL, U *_pR)
1426 {
1427     return sub_S_SC<T, U, O>(_pL, _pR);
1428 }
1429
1430 template<class T, class U, class O> InternalType* sub_SC_IC(T *_pL, U *_pR)
1431 {
1432     return sub_SC_SC<T, U, O>(_pL, _pR);
1433 }
1434
1435 template<class T, class U, class O> InternalType* sub_I_I(T *_pL, U *_pR)
1436 {
1437     O* pOut = types::Double::Identity(-1, -1);
1438     sub(_pL->get(0), _pR->get(0), pOut->get());
1439     return pOut;
1440 }
1441
1442 template<class T, class U, class O> InternalType* sub_I_IC(T *_pL, U *_pR)
1443 {
1444     O* pOut = types::Double::Identity(-1, -1);
1445     pOut->setComplex(true);
1446     sub(_pL->get(), (size_t)1, _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1447     return pOut;
1448 }
1449
1450 template<class T, class U, class O> InternalType* sub_IC_I(T *_pL, U *_pR)
1451 {
1452     O* pOut = types::Double::Identity(-1, -1);
1453     pOut->setComplex(true);
1454     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1455     return pOut;
1456 }
1457
1458 template<class T, class U, class O> InternalType* sub_IC_IC(T *_pL, U *_pR)
1459 {
1460     O* pOut = types::Double::Identity(-1, -1);
1461     pOut->setComplex(true);
1462     sub(_pL->get(0), _pL->getImg(0), _pR->get(0), _pR->getImg(0), pOut->get(), pOut->getImg());
1463     return pOut;
1464 }
1465
1466 template<class T, class U, class O> types::InternalType* sub_I_E(T *_pL, U * /*_pR*/)
1467 {
1468     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1469     Double* pOut = Double::Empty();
1470     sub();
1471     return pOut;
1472 }
1473
1474 template<class T, class U, class O> types::InternalType* sub_IC_E(T *_pL, U * /*_pR*/)
1475 {
1476     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1477     Double* pOut = Double::Empty();
1478     sub();
1479     return pOut;
1480 }
1481
1482 template<class T, class U, class O> types::InternalType* sub_E_I(T * /*_pL*/, U *_pR)
1483 {
1484     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1485     Double* pOut = Double::Empty();
1486     sub();
1487     return pOut;
1488 }
1489
1490 template<class T, class U, class O> types::InternalType* sub_E_IC(T * /*_pL*/, U *_pR)
1491 {
1492     Sciwarning(_("operation -: Warning adding a matrix with the empty matrix will give an empty matrix result.\n")); 
1493     Double* pOut = Double::Empty();
1494     sub();
1495     return pOut;
1496 }
1497
1498 template<> InternalType* sub_M_M<Polynom, Polynom, Polynom>(Polynom* _pL, Polynom* _pR)
1499 {
1500     Polynom* pOut = NULL;
1501     if (_pL->isScalar())
1502     {
1503         SinglePoly* p1Coef  = _pL->get(0);
1504         int iRank1          = p1Coef->getRank();
1505         int* pRank2         = new int[_pR->getSize()];
1506         int* pRankOut       = new int[_pR->getSize()];
1507
1508         _pR->getRank(pRank2);
1509         for (int i = 0 ; i < _pR->getSize() ; i++)
1510         {
1511             pRankOut[i] = std::max(iRank1, pRank2[i]);
1512         }
1513
1514         pOut = new Polynom(_pR->getVariableName(), _pR->getDims(), _pR->getDimsArray(), pRankOut);
1515         bool isOutComplex = _pL->isComplex() || _pR->isComplex();
1516
1517         //Result P1(0) - P2(i)
1518         double* p1R = p1Coef->get();
1519         if (isOutComplex)
1520         {
1521             double* p1I = p1Coef->getImg();
1522             for (int i = 0 ; i < _pR->getSize() ; i++)
1523             {
1524                 SinglePoly* p2Coef   = _pR->get(i);
1525                 double* p2R          = p2Coef->get();
1526                 double* p2I          = p2Coef->getImg();
1527
1528                 SinglePoly* pOutCoef = pOut->get(i);
1529                 pOutCoef->setComplex(isOutComplex);
1530                 double* pOutR        = pOutCoef->get();
1531                 double* pOutI        = pOutCoef->getImg();
1532
1533                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1534                 {
1535                     if (j > iRank1)
1536                     {
1537                         pOutR[j] = - p2R[j];
1538                         pOutI[j] = - (p2I ? p2I[j] : 0);
1539                     }
1540                     else if (j > pRank2[i])
1541                     {
1542                         pOutR[j] = p1R[j];
1543                         pOutI[j] = (p1I ? p1I[j] : 0);
1544                     }
1545                     else
1546                     {
1547                         pOutR[j] = p1R[j] - p2R[j];
1548                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
1549                     }
1550                 }
1551             }
1552         }
1553         else
1554         {
1555             for (int i = 0 ; i < _pR->getSize() ; i++)
1556             {
1557                 SinglePoly* p2Coef   = _pR->get(i);
1558                 double* p2R          = p2Coef->get();
1559
1560                 SinglePoly* pOutCoef = pOut->get(i);
1561                 double* pOutR        = pOutCoef->get();
1562
1563                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1564                 {
1565                     if (j > iRank1)
1566                     {
1567                         pOutR[j] = - p2R[j];
1568                     }
1569                     else if (j > pRank2[i])
1570                     {
1571                         pOutR[j] = p1R[j];
1572                     }
1573                     else
1574                     {
1575                         pOutR[j] = p1R[j] - p2R[j];
1576                     }
1577                 }
1578             }
1579         }
1580
1581         delete[] pRankOut;
1582         delete[] pRank2;
1583         pOut->updateRank();
1584         return pOut;
1585     }
1586
1587     if (_pR->isScalar())
1588     {
1589         //size(p2) == 1
1590         int *pRank1     = new int[_pL->getSize()];
1591         int iRank2      = _pR->get(0)->getRank();
1592         int *pRankOut   = new int[_pL->getSize()];
1593
1594         _pL->getRank(pRank1);
1595         for (int i = 0 ; i < _pL->getSize() ; i++)
1596         {
1597             pRankOut[i] = std::max(pRank1[i], iRank2);
1598         }
1599
1600         pOut = new Polynom(_pL->getVariableName(), _pL->getDims(), _pL->getDimsArray(), pRankOut);
1601         bool isOutComplex = _pL->isComplex() || _pR->isComplex();
1602
1603         //Result P1(i) - P2(0)
1604         SinglePoly *p2Coef          = _pR->get(0);
1605         double *p2R                 = p2Coef->get();
1606
1607         if (isOutComplex)
1608         {
1609             double *p2I             = p2Coef->getImg();
1610             for (int i = 0 ; i < _pL->getSize() ; i++)
1611             {
1612                 SinglePoly *p1Coef      = _pL->get(i);
1613                 double *p1R             = p1Coef->get();
1614                 double *p1I             = p1Coef->getImg();
1615
1616                 SinglePoly *pOutCoef    = pOut->get(i);
1617                 pOutCoef->setComplex(isOutComplex);
1618                 double *pOutR           = pOutCoef->get();
1619                 double *pOutI           = pOutCoef->getImg();
1620
1621                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1622                 {
1623                     if (j > pRank1[j])
1624                     {
1625                         pOutR[j] = - p2R[j];
1626                         pOutI[j] = - (p2I ? p2I[j] : 0);
1627                     }
1628                     else if (j > iRank2)
1629                     {
1630                         pOutR[j] = p1R[j];
1631                         pOutI[j] = (p1I ? p1I[j] : 0);
1632                     }
1633                     else
1634                     {
1635                         pOutR[j] = p1R[j] - p2R[j];
1636                         pOutI[j] = (p1I ? p1I[j] : 0) - (p2I ? p2I[j] : 0);
1637                     }
1638                 }
1639             }
1640         }
1641         else
1642         {
1643             for (int i = 0 ; i < _pL->getSize() ; i++)
1644             {
1645                 SinglePoly *p1Coef      = _pL->get(i);
1646                 double *p1R             = p1Coef->get();
1647
1648                 SinglePoly *pOutCoef    = pOut->get(i);
1649                 double *pOutR           = pOutCoef->get();
1650
1651                 for (int j = 0 ; j < pRankOut[i] + 1 ; j++)
1652                 {
1653                     if (j > pRank1[j])
1654                     {
1655                         pOutR[j] = - p2R[j];
1656                     }
1657                     else if (j > iRank2)
1658                     {
1659                         pOutR[j] = p1R[j];
1660                     }
1661                     else
1662                     {
1663                         pOutR[j] = p1R[j] - p2R[j];
1664                     }
1665                 }
1666             }
1667         }
1668
1669         delete[] pRankOut;
1670         delete[] pRank1;
1671         pOut->updateRank();
1672         return pOut;
1673     }
1674
1675     //check dimension compatibilities
1676     int iDims1  = _pL->getDims();
1677     int iDims2  = _pR->getDims();
1678
1679     if (iDims1 != iDims2)
1680     {
1681         return nullptr;
1682     }
1683
1684     int* piDims1    = _pL->getDimsArray();
1685     int* piDims2    = _pR->getDimsArray();
1686
1687     for (int i = 0 ; i < iDims1 ; i++)
1688     {
1689         if (piDims1[i] != piDims2[i])
1690         {
1691             throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
1692         }
1693     }
1694
1695     int *pRankOut   = new int[_pL->getSize()];
1696     int *pRank1     = new int[_pL->getSize()];
1697     int *pRank2     = new int[_pR->getSize()];
1698
1699     _pL->getRank(pRank1);
1700     _pR->getRank(pRank2);
1701     for (int i = 0 ; i < _pL->getSize() ; i++)
1702     {
1703         pRankOut[i] = std::max(pRank1[i], pRank2[i]);
1704     }
1705
1706     pOut = new Polynom(_pR->getVariableName(), iDims1, piDims1, pRankOut);
1707     bool isOutComplex =  _pL->isComplex() || _pR->isComplex();
1708
1709     //Result P1(i) - P2(i)
1710     for (int i = 0 ; i < _pL->getSize() ; i++)
1711     {
1712         SinglePoly* pOutCoef = pOut->get(i);
1713         pOutCoef->setComplex(isOutComplex);
1714
1715         double *p1R     = _pL->get(i)->get();
1716         double *p2R     = _pR->get(i)->get();
1717         double *pOutR   = pOutCoef->get();
1718         int iMin        = std::min(pRank1[i], pRank2[i]);
1719         int iMax        = std::max(pRank1[i], pRank2[i]);
1720
1721         for (int j = 0 ; j < iMin + 1 ; j++)
1722         {
1723             pOutR[j]    = p1R[j] - p2R[j];
1724         }
1725
1726         double *pTemp   = NULL;
1727         int iCoef       = 1;
1728         if (pRank1[i] > pRank2[i])
1729         {
1730             pTemp       = p1R;
1731             iCoef       = 1;
1732         }
1733         else
1734         {
1735             pTemp       = p2R;
1736             iCoef       = -1;
1737         }
1738
1739         for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1740         {
1741             pOutR[j]    = pTemp[j] * iCoef;
1742         }
1743
1744         if (isOutComplex)
1745         {
1746             double *p1I     = _pL->get(i)->getImg();
1747             double *p2I     = _pR->get(i)->getImg();
1748             double *pOutI   = pOutCoef->getImg();
1749
1750             for (int j = 0 ; j < iMin + 1 ; j++)
1751             {
1752                 pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
1753             }
1754
1755             for (int j = iMin + 1 ; j < iMax + 1 ; j++)
1756             {
1757                 pOutI[j]  = pTemp[j] * iCoef;
1758             }
1759         }
1760     }
1761
1762     delete[] pRankOut;
1763     delete[] pRank1;
1764     delete[] pRank2;
1765
1766     pOut->updateRank();
1767     return pOut;
1768 }
1769
1770 // P - D
1771 template<> InternalType* sub_M_M<Polynom, Double, Polynom>(Polynom* _pL, Double* _pR)
1772 {
1773     Polynom* pOut = NULL;
1774
1775     bool bComplex1 = _pL->isComplex();
1776     bool bComplex2 = _pR->isComplex();
1777     bool isComplexOut = bComplex2 || bComplex1;
1778
1779     // P - []
1780     if (_pR->isEmpty())
1781     {
1782         return _pL;
1783     }
1784
1785     if (_pR->isIdentity())
1786     {
1787         double dblRightR = _pR->get(0);
1788         double dblRightI = 0;
1789         if (_pR->isComplex())
1790         {
1791             dblRightI = _pR->getImg(0);
1792         }
1793
1794         int iDims = _pL->getDims();
1795         int* piDims = _pL->getDimsArray();
1796         pOut = (Polynom*)_pL->clone();
1797         pOut->setComplex(isComplexOut);
1798         SinglePoly** pSPOut = pOut->get();
1799         SinglePoly** pSPLeft = _pL->get();
1800         int iLeadDims = piDims[0];
1801         int* piIndex = new int[iDims];
1802         piIndex[0] = 0;
1803
1804         //find smaller dims
1805         for (int i = 1 ; i < iDims ; ++i)
1806         {
1807             //init
1808             piIndex[i] = 0;
1809
1810             if (iLeadDims > piDims[i])
1811             {
1812                 iLeadDims = piDims[i];
1813             }
1814         }
1815
1816         if (isComplexOut)
1817         {
1818             for (int i = 0 ; i < iLeadDims ; ++i)
1819             {
1820                 for (int j = 0 ; j < iDims ; ++j)
1821                 {
1822                     piIndex[j] = i;
1823                 }
1824
1825                 int index = _pL->getIndex(piIndex);
1826                 sub(pSPLeft[index]->get(0), pSPLeft[index]->getImg(0), (size_t)1, &dblRightR, &dblRightI, pSPOut[index]->get(), pSPOut[index]->getImg());
1827             }
1828         }
1829         else
1830         {
1831             for (int i = 0 ; i < iLeadDims ; ++i)
1832             {
1833                 for (int j = 0 ; j < iDims ; ++j)
1834                 {
1835                     piIndex[j] = i;
1836                 }
1837
1838                 int index = _pL->getIndex(piIndex);
1839                 sub(pSPLeft[index]->get(0), (size_t)1, &dblRightR, pSPOut[index]->get());
1840             }
1841         }
1842
1843         delete[] piIndex;
1844         return pOut;
1845
1846     }
1847
1848     if (_pR->isScalar())
1849     {
1850         //subtract same value to all polynoms
1851         pOut = (Polynom*)_pL->clone();
1852         pOut->setComplex(isComplexOut);
1853
1854         SinglePoly** pSP = pOut->get();
1855         int iSize = pOut->getSize();
1856
1857         double dblR = _pR->get(0);
1858         if (isComplexOut)
1859         {
1860             double dblI = _pR->getImg(0);
1861             for (int i = 0 ; i < iSize ; ++i)
1862             {
1863                 pSP[i]->get()[0] -= dblR;
1864                 pSP[i]->getImg()[0] -= dblI;
1865             }
1866         }
1867         else
1868         {
1869             for (int i = 0 ; i < iSize ; ++i)
1870             {
1871                 pSP[i]->get()[0] -= dblR;
1872             }
1873         }
1874
1875         return pOut;
1876     }
1877
1878     if (_pL->isScalar())
1879     {
1880         //and _pR is not !
1881         int iDims = _pR->getDims();
1882         int* piDims = _pR->getDimsArray();
1883         int iSize = _pR->getSize();
1884
1885         pOut = new Polynom(_pL->getVariableName(), iDims, piDims);
1886
1887         SinglePoly* pSPL = _pL->get(0);
1888         if (_pR->isComplex())
1889         {
1890             double* pdblR = _pR->get();
1891             double* pdblI = _pR->getImg();
1892
1893             for (int i = 0 ; i < iSize ; ++i)
1894             {
1895                 SinglePoly* pSPOut = pSPL->clone();
1896                 //in case of original is not complex
1897                 pSPOut->setComplex(isComplexOut);
1898                 pSPOut->get()[0] -= pdblR[i];
1899                 pSPOut->getImg()[0] -= pdblI[i];
1900                 pOut->set(i, pSPOut);
1901                 delete pSPOut;
1902             }
1903         }
1904         else
1905         {
1906             double* pdblR = _pR->get();
1907
1908             for (int i = 0 ; i < iSize ; ++i)
1909             {
1910                 SinglePoly* pSPOut = pSPL->clone();
1911                 //update 0th rank value
1912                 pSPOut->get()[0] -= pdblR[i];
1913                 pOut->set(i, pSPOut);
1914                 delete pSPOut;
1915             }
1916         }
1917
1918         return pOut;
1919     }
1920
1921     //P - D
1922
1923     //check dimensions
1924     int iDims1 = _pR->getDims();
1925     int iDims2 = _pL->getDims();
1926
1927     if (iDims1 != iDims2)
1928     {
1929         return nullptr;
1930     }
1931
1932     int* piDims1 = _pR->getDimsArray();
1933     int* piDims2 = _pL->getDimsArray();
1934
1935     for (int i = 0 ; i < iDims1 ; i++)
1936     {
1937         if (piDims1[i] != piDims2[i])
1938         {
1939             wchar_t pMsg[bsiz];
1940             os_swprintf(pMsg, bsiz, _W("Error: operator %ls: Matrix dimensions must agree (op1 is %ls, op2 is %ls).\n").c_str(),  L"-", _pL->DimToString().c_str(), _pR->DimToString().c_str());
1941             throw ast::InternalError(pMsg);
1942         }
1943     }
1944
1945     double* pInDblR = _pR->get();
1946     pOut = (Polynom*)_pL->clone();
1947     if (bComplex1 && bComplex2)
1948     {
1949         double* pInDblI = _pR->getImg();
1950         for (int i = 0 ; i < pOut->getSize() ; i++)
1951         {
1952             SinglePoly *pSPOut   = pOut->get(i);
1953             double *pOutPolyR    = pSPOut->get();
1954             double *pOutPolyI    = pSPOut->getImg();
1955
1956             pOutPolyR[0] -= pInDblR[i];
1957             pOutPolyI[0] -= pInDblI[i];
1958         }
1959     }
1960     else if (bComplex2)
1961     {
1962         double* pInDblI = _pR->getImg();
1963         pOut->setComplex(true);
1964         for (int i = 0 ; i < pOut->getSize() ; i++)
1965         {
1966             SinglePoly *pSPOut   = pOut->get(i);
1967             double *pOutPolyR    = pSPOut->get();
1968             double *pOutPolyI    = pSPOut->getImg();
1969
1970             pOutPolyR[0] -= pInDblR[i];
1971             pOutPolyI[0] = -pInDblI[i];
1972         }
1973     }
1974     else
1975     {
1976         for (int i = 0 ; i < pOut->getSize() ; i++)
1977         {
1978             SinglePoly *pSPOut = pOut->get(i);
1979             double *pOutPolyR  = pSPOut->get();
1980
1981             pOutPolyR[0] -= pInDblR[i];
1982         }
1983     }
1984
1985     return pOut;
1986 }
1987
1988 template<> InternalType* sub_I_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
1989 {
1990     Polynom* pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
1991     double dblLeft = _pL->get(0);
1992
1993     int iDims = _pR->getDims();
1994     int* piDims = _pR->getDimsArray();
1995     int iLeadDims = piDims[0];
1996     int* piIndex = new int[iDims];
1997     SinglePoly** pSP = _pR->get();
1998     SinglePoly** pSPOut = pOut->get();
1999     piIndex[0] = 0;
2000
2001     //find smaller dims
2002     for (int i = 1 ; i < iDims ; ++i)
2003     {
2004         //init
2005         piIndex[i] = 0;
2006
2007         if (iLeadDims > piDims[i])
2008         {
2009             iLeadDims = piDims[i];
2010         }
2011     }
2012
2013     for (int i = 0 ; i < iLeadDims ; ++i)
2014     {
2015         for (int j = 0 ; j < iDims ; ++j)
2016         {
2017             piIndex[j] = i;
2018         }
2019
2020         int index = _pR->getIndex(piIndex);
2021         sub(dblLeft, pSP[index]->get(0), pSPOut[index]->get());
2022     }
2023
2024     delete[] piIndex;
2025     return pOut;
2026 }
2027
2028 template<> InternalType* sub_I_MC<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2029 {
2030     Polynom* pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2031     double dblLeft = _pL->get(0);
2032
2033     int iDims = _pR->getDims();
2034     int* piDims = _pR->getDimsArray();
2035     int iLeadDims = piDims[0];
2036     int* piIndex = new int[iDims];
2037     SinglePoly** pSP = _pR->get();
2038     SinglePoly** pSPOut = pOut->get();
2039     piIndex[0] = 0;
2040
2041     for (int i = 0 ; i < iDims ; ++i)
2042     {
2043         //init
2044         piIndex[i] = 0;
2045
2046         if (iLeadDims > piDims[i])
2047         {
2048             iLeadDims = piDims[i];
2049         }
2050     }
2051
2052     for (int i = 0 ; i < iLeadDims ; ++i)
2053     {
2054         for (int j = 0 ; j < iDims ; ++j)
2055         {
2056             piIndex[j] = i;
2057         }
2058
2059         int index = _pR->getIndex(piIndex);
2060         sub(dblLeft, pSP[index]->get(0), pSPOut[index]->get());
2061     }
2062
2063     delete[] piIndex;
2064     return pOut;
2065 }
2066
2067 template<> InternalType* sub_IC_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2068 {
2069     Polynom* pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2070     pOut->setComplex(true);
2071     double dblLeftR = _pL->get(0);
2072     double dblLeftI = _pL->getImg(0);
2073
2074     int iDims = _pR->getDims();
2075     int* piDims = _pR->getDimsArray();
2076     int iLeadDims = piDims[0];
2077     int* piIndex = new int[iDims];
2078     SinglePoly** pSP = _pR->get();
2079     SinglePoly** pSPOut = pOut->get();
2080     piIndex[0] = 0;
2081
2082     for (int i = 0 ; i < iDims ; ++i)
2083     {
2084         //init
2085         piIndex[i] = 0;
2086
2087         if (iLeadDims > piDims[i])
2088         {
2089             iLeadDims = piDims[i];
2090         }
2091     }
2092
2093     for (int i = 0 ; i < iLeadDims ; ++i)
2094     {
2095         for (int j = 0 ; j < iDims ; ++j)
2096         {
2097             piIndex[j] = i;
2098         }
2099
2100         int index = _pR->getIndex(piIndex);
2101         sub(&dblLeftR, &dblLeftI, (size_t)1, pSP[index]->get(0), pSPOut[index]->get(), pSPOut[index]->getImg());
2102     }
2103
2104     delete[] piIndex;
2105     return pOut;
2106 }
2107
2108 template<> InternalType* sub_IC_MC<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2109 {
2110     Polynom* pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2111     double dblLeftR = _pL->get(0);
2112     double dblLeftI = _pL->getImg(0);
2113
2114     int iDims = _pR->getDims();
2115     int* piDims = _pR->getDimsArray();
2116     int iLeadDims = piDims[0];
2117     int* piIndex = new int[iDims];
2118     SinglePoly** pSP = _pR->get();
2119     SinglePoly** pSPOut = pOut->get();
2120     piIndex[0] = 0;
2121
2122     for (int i = 0 ; i < iDims ; ++i)
2123     {
2124         //init
2125         piIndex[i] = 0;
2126
2127         if (iLeadDims > piDims[i])
2128         {
2129             iLeadDims = piDims[i];
2130         }
2131     }
2132
2133     for (int i = 0 ; i < iLeadDims ; ++i)
2134     {
2135         for (int j = 0 ; j < iDims ; ++j)
2136         {
2137             piIndex[j] = i;
2138         }
2139
2140         int index = _pR->getIndex(piIndex);
2141         sub(dblLeftR, dblLeftI, pSP[index]->get(0), pSP[index]->getImg(0), pSPOut[index]->get(), pSPOut[index]->getImg());
2142     }
2143
2144     delete[] piIndex;
2145     return pOut;
2146 }
2147
2148 template<> InternalType* sub_M_M<Double, Polynom, Polynom>(Double* _pL, Polynom* _pR)
2149 {
2150     Polynom* pOut = NULL;
2151     bool bComplex1 = _pL->isComplex();
2152     bool bComplex2 = _pR->isComplex();
2153
2154     double *pInDblR = _pL->getReal();
2155     double *pInDblI = _pL->getImg();
2156
2157     if (_pL->isEmpty())
2158     {
2159         return _pR;
2160     }
2161
2162     if (_pR->isScalar())
2163     {
2164         int *piRank = new int[_pL->getSize()];
2165         for (int i = 0 ; i < _pL->getSize() ; i++)
2166         {
2167             piRank[i] = _pR->get(0)->getRank();
2168         }
2169
2170         pOut = new Polynom(_pR->getVariableName(), _pL->getDims(), _pL->getDimsArray(), piRank);
2171         delete[] piRank;
2172         if (bComplex1 || bComplex2)
2173         {
2174             pOut->setComplex(true);
2175         }
2176
2177         for (int i = 0 ; i < pOut->getSize() ; i++)
2178         {
2179             SinglePoly *pInPoly  = _pR->get(0);
2180             SinglePoly *pOutPoly = pOut->get(i);
2181             double *pInPolyR     = pInPoly->get();
2182             double *pOutPolyR    = pOutPoly->get();
2183
2184             pOutPolyR[0] = pInDblR[i] - pInPolyR[0];
2185
2186             for (int j = 1 ; j < pInPoly->getSize() ; j++)
2187             {
2188                 pOutPolyR[j] = -pInPolyR[j];
2189             }
2190         }
2191
2192         if (pOut->isComplex())
2193         {
2194             for (int i = 0 ; i < pOut->getSize() ; i++)
2195             {
2196                 SinglePoly *pInPoly  = _pR->get(0);
2197                 SinglePoly *pOutPoly = pOut->get(i);
2198                 double *pInPolyI     = pInPoly->getImg();
2199                 double *pOutPolyI    = pOutPoly->getImg();
2200
2201                 pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pInPolyI != NULL ? pInPolyI[0] : 0);
2202
2203                 for (int j = 1 ; j < pInPoly->getSize() ; j++)
2204                 {
2205                     pOutPolyI[j] = (pInPolyI != NULL ? -pInPolyI[j] : 0);
2206                 }
2207             }
2208         }
2209
2210         return pOut;
2211     }
2212
2213     if (_pL->isScalar())
2214     {
2215         if (bComplex2)
2216         {
2217             pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2218         }
2219         else
2220         {
2221             pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2222         }
2223
2224         int iSize = pOut->getSize();
2225
2226         if (bComplex1 && bComplex2)
2227         {
2228             for (int i = 0 ; i < iSize ; i++)
2229             {
2230                 SinglePoly *pSPOut   = pOut->get(i);
2231                 double *pOutPolyR    = pSPOut->get();
2232                 double *pOutPolyI    = pSPOut->getImg();
2233
2234                 pOutPolyR[0] += pInDblR[0];
2235                 pOutPolyI[0] += pInDblI[0];
2236             }
2237         }
2238         else if (bComplex1)
2239         {
2240             pOut->setComplex(true);
2241             for (int i = 0 ; i < iSize ; i++)
2242             {
2243                 SinglePoly *pSPOut   = pOut->get(i);
2244                 double *pOutPolyR    = pSPOut->get();
2245                 double *pOutPolyI    = pSPOut->getImg();
2246
2247                 pOutPolyR[0] += pInDblR[0];
2248                 pOutPolyI[0] = pInDblI[0];
2249             }
2250         }
2251         else
2252         {
2253             for (int i = 0 ; i < iSize ; i++)
2254             {
2255                 SinglePoly *pSPOut = pOut->get(i);
2256                 double *pOutPolyR  = pSPOut->get();
2257
2258                 pOutPolyR[0] += pInDblR[0];
2259             }
2260         }
2261
2262         return pOut;
2263     }
2264
2265     int iDims1 = _pR->getDims();
2266     int iDims2 = _pL->getDims();
2267
2268     if (iDims1 != iDims2)
2269     {
2270         return nullptr;
2271     }
2272
2273     int* piDims1 = _pR->getDimsArray();
2274     int* piDims2 = _pL->getDimsArray();
2275
2276     for (int i = 0 ; i < iDims1 ; i++)
2277     {
2278         if (piDims1[i] != piDims2[i])
2279         {
2280             wchar_t pMsg[bsiz];
2281             os_swprintf(pMsg, bsiz, _W("Error: operator %ls: Matrix dimensions must agree (op1 is %ls, op2 is %ls).\n").c_str(),  L"+", _pL->DimToString().c_str(), _pR->DimToString().c_str());
2282             throw ast::InternalError(pMsg);
2283         }
2284     }
2285
2286     if (bComplex2)
2287     {
2288         pOut = (Polynom*)opposite_MC<Polynom, Polynom>(_pR);
2289     }
2290     else
2291     {
2292         pOut = (Polynom*)opposite_M<Polynom, Polynom>(_pR);
2293     }
2294
2295     if (bComplex1 && bComplex2)
2296     {
2297         for (int i = 0 ; i < pOut->getSize() ; i++)
2298         {
2299             SinglePoly *pSPOut   = pOut->get(i);
2300             double *pOutPolyR    = pSPOut->get();
2301             double *pOutPolyI    = pSPOut->getImg();
2302
2303             pOutPolyR[0] += pInDblR[i];
2304             pOutPolyI[0] += pInDblI[i];
2305         }
2306     }
2307     else if (bComplex1)
2308     {
2309         pOut->setComplex(true);
2310         for (int i = 0 ; i < pOut->getSize() ; i++)
2311         {
2312             SinglePoly *pSPOut   = pOut->get(i);
2313             double *pOutPolyR    = pSPOut->get();
2314             double *pOutPolyI    = pSPOut->getImg();
2315
2316             pOutPolyR[0] += pInDblR[i];
2317             pOutPolyI[0] = pInDblI[i];
2318         }
2319     }
2320     else
2321     {
2322         for (int i = 0 ; i < pOut->getSize() ; i++)
2323         {
2324             SinglePoly *pSPOut = pOut->get(i);
2325             double *pOutPolyR  = pSPOut->get();
2326
2327             pOutPolyR[0] += pInDblR[i];
2328         }
2329     }
2330
2331     return pOut;
2332 }
2333
2334 //sp - sp
2335 template<> InternalType* sub_M_M<Sparse, Sparse, Sparse>(Sparse* _pL, Sparse* _pR)
2336 {
2337     //check scalar hidden in a sparse ;)
2338     if (_pL->isScalar() || _pR->isScalar())
2339     {
2340         // scalar - sp  or  sp - scalar
2341         // call Overload
2342         return NULL;
2343     }
2344
2345     if (_pL->getRows() != _pR->getRows() || _pL->getCols() != _pR->getCols())
2346     {
2347         //dimensions not match
2348         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2349     }
2350
2351     types::Sparse* pSPOut = _pL->substract(*_pR);
2352     pSPOut->finalize();
2353     return pSPOut;
2354 }
2355
2356 //d - sp
2357 template<> InternalType* sub_M_M<Double, Sparse, Double>(Double* _pL, Sparse* _pR)
2358 {
2359     Double* pOut = NULL;
2360     int iOne = 1; //fortran
2361     bool bComplex1 = _pL->isComplex();
2362     bool bComplex2 = _pR->isComplex();
2363
2364     if (_pL->isScalar() && _pR->isScalar())
2365     {
2366         //d - sp
2367         pOut = (Double*)_pL->clone();
2368         pOut->setComplex(bComplex1 || bComplex2);
2369         if (bComplex2)
2370         {
2371             std::complex<double> dbl = _pR->getImg(0, 0);
2372             pOut->set(0, dbl.real() - pOut->get(0));
2373             pOut->setImg(0, pOut->getImg(0) - dbl.imag());
2374         }
2375         else
2376         {
2377             pOut->set(0, pOut->get(0) - _pL->get(0, 0));
2378         }
2379
2380         return pOut;
2381     }
2382
2383     if (_pL->isScalar())
2384     {
2385         //d - SP
2386         pOut = new Double(_pR->getRows(), _pR->getCols(), bComplex1 || bComplex2);
2387         int iSize = _pR->getSize();
2388         double dblVal = _pL->get(0);
2389         double dblValI = 0;
2390         C2F(dset)(&iSize, &dblVal, pOut->get(), &iOne);
2391         if (bComplex1)
2392         {
2393             //initialize imag part at 0
2394             dblValI = _pL->getImg(0);
2395             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2396         }
2397         else if (bComplex2)
2398         {
2399             dblValI = 0;
2400             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2401         }
2402
2403         int nonZeros = static_cast<int>(_pR->nonZeros());
2404         int* pRows = new int[nonZeros * 2];
2405         _pR->outputRowCol(pRows);
2406         int* pCols = pRows + nonZeros;
2407
2408         if (pOut->isComplex())
2409         {
2410             for (int i = 0 ; i < nonZeros ; i++)
2411             {
2412                 int iRow = static_cast<int>(pRows[i]) - 1;
2413                 int iCol = static_cast<int>(pCols[i]) - 1;
2414                 std::complex<double> dbl = _pR->getImg(iRow, iCol);
2415                 pOut->set(iRow, iCol, dblVal - dbl.real());
2416                 pOut->setImg(iRow, iCol, dblValI - dbl.imag());
2417             }
2418         }
2419         else
2420         {
2421             for (int i = 0 ; i < nonZeros ; i++)
2422             {
2423                 int iRow = static_cast<int>(pRows[i]) - 1;
2424                 int iCol = static_cast<int>(pCols[i]) - 1;
2425                 pOut->set(iRow, iCol, dblVal - _pR->get(iRow, iCol));
2426             }
2427         }
2428
2429         //clear
2430         delete[] pRows;
2431
2432         return pOut;
2433     }
2434
2435     if (_pR->isScalar())
2436     {
2437         //D - sp
2438         pOut = (Double*)_pR->clone();
2439         pOut->setComplex(bComplex1 || bComplex2);
2440
2441         if (pOut->isComplex())
2442         {
2443             double* pReal = pOut->get();
2444             double* pImg = pOut->getImg();
2445             int size = pOut->getSize();
2446             for (int i = 0 ; i < size ; i++)
2447             {
2448                 std::complex<double> dbl = _pR->getImg(0, 0);
2449                 pReal[i] -= dbl.real();
2450                 pImg[i] -= dbl.imag();
2451             }
2452         }
2453         else
2454         {
2455             double* pReal = pOut->get();
2456             int size = pOut->getSize();
2457             for (int i = 0 ; i < size ; i++)
2458             {
2459                 pReal[i] -= pReal[i];
2460             }
2461         }
2462
2463         return pOut;
2464     }
2465
2466     if (_pL->getDims() > 2)
2467     {
2468         return nullptr;
2469     }
2470
2471     if (_pL->getRows() == _pR->getRows() && _pL->getCols() == _pR->getCols())
2472     {
2473         //D - SP
2474         pOut = (Double*)_pL->clone();
2475         pOut->setComplex(bComplex1 || bComplex2);
2476
2477         int nonZeros = static_cast<int>(_pR->nonZeros());
2478         int* pRows = new int[nonZeros * 2];
2479         _pR->outputRowCol(pRows);
2480         int* pCols = pRows + nonZeros;
2481
2482         if (bComplex1 || bComplex2)
2483         {
2484             for (int i = 0 ; i < nonZeros ; i++)
2485             {
2486                 int iRow = static_cast<int>(pRows[i]) - 1;
2487                 int iCol = static_cast<int>(pCols[i]) - 1;
2488                 std::complex<double> dbl = _pR->getImg(iRow, iCol);
2489                 pOut->set(iRow, iCol, pOut->get(iRow, iCol) - dbl.real());
2490                 pOut->setImg(iRow, iCol, pOut->getImg(iRow, iCol) - dbl.imag());
2491             }
2492         }
2493         else
2494         {
2495             for (int i = 0 ; i < nonZeros ; i++)
2496             {
2497                 int iRow = static_cast<int>(pRows[i]) - 1;
2498                 int iCol = static_cast<int>(pCols[i]) - 1;
2499                 pOut->set(iRow, iCol, pOut->get(iRow, iCol) - _pR->get(iRow, iCol));
2500             }
2501         }
2502
2503         //clear
2504         delete[] pRows;
2505         return pOut;
2506     }
2507     else
2508     {
2509         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2510     }
2511 }
2512
2513 //sp - d
2514 template<> InternalType* sub_M_M<Sparse, Double, Double>(Sparse* _pL, Double* _pR)
2515 {
2516     Double* pOut = NULL;
2517     int iOne = 1; //fortran
2518     bool bComplex1 = _pL->isComplex();
2519     bool bComplex2 = _pR->isComplex();
2520     bool bComplexOut = bComplex1 || bComplex2;
2521
2522     if (_pL->isScalar() && _pR->isScalar())
2523     {
2524         //sp - d
2525         pOut = (Double*)_pR->clone();
2526         pOut->setComplex(bComplexOut);
2527         if (bComplex1)
2528         {
2529             std::complex<double> dbl = _pL->getImg(0, 0);
2530             pOut->set(0, pOut->get(0) - dbl.real());
2531             pOut->setImg(0, dbl.imag() - pOut->getImg(0));
2532         }
2533         else
2534         {
2535             pOut->set(0, _pL->get(0, 0) - pOut->get(0));
2536         }
2537
2538         return pOut;
2539     }
2540
2541     if (_pR->isScalar())
2542     {
2543         //SP - d
2544         pOut = new Double(_pL->getRows(), _pL->getCols(), bComplexOut);
2545         int iSize = _pL->getSize();
2546         double dblVal = -_pR->get(0);
2547         double dblValI = 0;
2548         C2F(dset)(&iSize, &dblVal, pOut->get(), &iOne);
2549         if (bComplex2)
2550         {
2551             dblValI = -_pR->getImg(0);
2552             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2553         }
2554         else if (bComplex1)
2555         {
2556             //initialize imag part at 0
2557             dblValI = 0;
2558             C2F(dset)(&iSize, &dblValI, pOut->getImg(), &iOne);
2559         }
2560
2561         int nonZeros = static_cast<int>(_pL->nonZeros());
2562         int* pRows = new int[nonZeros * 2];
2563         _pL->outputRowCol(pRows);
2564         int* pCols = pRows + nonZeros;
2565
2566         if (bComplex1)
2567         {
2568             for (int i = 0 ; i < nonZeros ; i++)
2569             {
2570                 int iRow = static_cast<int>(pRows[i]) - 1;
2571                 int iCol = static_cast<int>(pCols[i]) - 1;
2572                 std::complex<double> dbl = _pL->getImg(iRow, iCol);
2573                 pOut->set(iRow, iCol, dbl.real() - (-dblVal));
2574                 pOut->setImg(iRow, iCol, dbl.imag() - (-dblValI));
2575             }
2576         }
2577         else
2578         {
2579             for (int i = 0 ; i < nonZeros ; i++)
2580             {
2581                 int iRow = static_cast<int>(pRows[i]) - 1;
2582                 int iCol = static_cast<int>(pCols[i]) - 1;
2583                 pOut->set(iRow, iCol, _pL->get(iRow, iCol) - (-dblVal));
2584             }
2585         }
2586
2587         //clear
2588         delete[] pRows;
2589
2590         return pOut;
2591     }
2592
2593     if (_pL->isScalar())
2594     {
2595         //sp - D
2596         pOut = (Double*)_pR->clone();
2597         pOut->setComplex(bComplexOut);
2598
2599         if (bComplex1)
2600         {
2601             double* pReal = pOut->get();
2602             double* pImg = pOut->getImg();
2603             int size = pOut->getSize();
2604             for (int i = 0 ; i < size ; i++)
2605             {
2606                 std::complex<double> dbl = _pL->getImg(0, 0);
2607                 pReal[i] = dbl.real() - pReal[i];
2608                 pImg[i] = dbl.imag() - pImg[i];
2609             }
2610         }
2611         else
2612         {
2613             double* pReal = pOut->get();
2614             int size = pOut->getSize();
2615             for (int i = 0 ; i < size ; i++)
2616             {
2617                 pReal[i] = _pL->get(0, 0) - pReal[i];
2618             }
2619         }
2620
2621         return pOut;
2622     }
2623
2624
2625     if (_pR->getDims() > 2)
2626     {
2627         return nullptr;
2628     }
2629
2630     if (_pL->getRows() != _pR->getRows() || _pL->getCols() != _pR->getCols())
2631     {
2632         throw ast::InternalError(_W("Inconsistent row/column dimensions.\n"));
2633     }
2634
2635     //SP - D
2636     pOut = new types::Double(_pL->getRows(), _pL->getCols(), _pL->isComplex());
2637     _pL->fill(*pOut);
2638     int iSize = pOut->getSize();
2639     pOut->setComplex(bComplexOut);
2640     double* pdblOutR = pOut->get();
2641     double* pdblOutI = pOut->getImg(); //can be null for non complex output
2642     double* pdblRR = _pR->get();
2643     double* pdblRI = _pR->getImg(); //can be null for non complex output
2644
2645     if (bComplex1)
2646     {
2647         if (bComplex2)
2648         {
2649             sub(pdblOutR, pdblOutI, iSize, pdblRR, pdblRI, pdblOutR, pdblOutI);
2650         }
2651         else
2652         {
2653             sub(pdblOutR, pdblOutI, iSize, pdblRR, pdblOutR, pdblOutI);
2654         }
2655     }
2656     else
2657     {
2658         if (bComplex2)
2659         {
2660             sub(pdblOutR, iSize, pdblRR, pdblRI, pdblOutR, pdblOutI);
2661         }
2662         else
2663         {
2664             sub(pdblOutR, iSize, pdblRR, pdblOutR);
2665         }
2666     }
2667
2668     return pOut;
2669 }
2670
2671 //[] - sp
2672 template<> InternalType* sub_M_M<Double, Sparse, Sparse>(Double* _pL, Sparse* _pR)
2673 {
2674     return sub_M_M<Sparse, Double, Sparse>(_pR, _pL);
2675 }
2676
2677 //sp - []
2678 template<> InternalType* sub_M_M<Sparse, Double, Sparse>(Sparse* _pL, Double* _pR)
2679 {
2680     Sparse* pOut = NULL;
2681     if (_pR->isIdentity())
2682     {
2683         //convert to _pL
2684         Sparse* pS = new Sparse(_pL->getRows(), _pL->getCols(), _pR->isComplex());
2685         if (pS->isComplex())
2686         {
2687             int size = std::min(_pL->getRows(), _pL->getCols());
2688             for (int i = 0 ; i < size ; i++)
2689             {
2690                 pS->set(i, i, std::complex<double>(_pR->get(0), _pR->getImg(0)));
2691             }
2692         }
2693         else
2694         {
2695             int size = std::min(_pL->getRows(), _pL->getCols());
2696             for (int i = 0 ; i < size ; i++)
2697             {
2698                 pS->set(i, i, _pR->get(0));
2699             }
2700         }
2701
2702         //AddSparseToSparse(_pL, pS, (Sparse**)pOut);
2703         delete pS;
2704         return pOut;
2705     }
2706     else
2707     {
2708         //is []
2709         return _pL;
2710     }
2711 }