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