d3ba531f4b17a6029f819fc82e8fd0ae3634f534
[scilab.git] / scilab / modules / ast / src / cpp / operations / types_substraction.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 <stdio.h>
15
16 #include "types_substraction.hxx"
17 //#include "scilabexception.hxx"
18 //#include "core_math.h"
19 #include "double.hxx"
20 #include "polynom.hxx"
21 #include "int.hxx"
22 #include "sparse.hxx"
23 #include "generic_operations.hxx"
24
25 //
26 extern "C"
27 {
28 #include "matrix_substraction.h"
29     //#include "localization.h"
30     //#include "charEncoding.h"
31     //#include "os_swprintf.h"
32 #include "elem_common.h" //dset
33 }
34
35 using namespace types;
36
37 InternalType* GenericUnaryMinus(InternalType* _pRightOperand)
38 {
39     // - Double
40     if (_pRightOperand->isDouble())
41     {
42         Double *pR = dynamic_cast<Double*>(_pRightOperand->clone());
43         bool bComplex = pR->isComplex();
44         double* pReal = pR->get();
45
46         if (bComplex)
47         {
48             double* pImg = pR->getImg();
49             for (int i = 0 ; i < pR->getSize() ; i++)
50             {
51                 pReal[i] *= -1;
52                 pImg[i]  *= -1;
53             }
54         }
55         else
56         {
57             for (int i = 0 ; i < pR->getSize() ; i++)
58             {
59                 pReal[i] *= -1;
60             }
61         }
62
63         return pR;
64     }
65
66     // - Polynom
67     if (_pRightOperand->isPoly())
68     {
69         double* pReal = NULL;
70         double* pImg  = NULL;
71
72         Double* pDblCoef = NULL;
73
74         Polynom *pR = dynamic_cast<Polynom*>(_pRightOperand->clone());
75         bool bComplex = pR->isComplex();
76
77         if (bComplex)
78         {
79             for (int i = 0; i < pR->getSize(); i++)
80             {
81                 pDblCoef = pR->get(i)->getCoef();
82                 pReal = pDblCoef->getReal();
83                 pImg = pDblCoef->getImg();
84
85                 for (int j = 0; j < pDblCoef->getSize(); j++)
86                 {
87                     pReal[j] *= -1;
88                     pImg[j] *= -1;
89                 }
90             }
91         }
92         else
93         {
94             for (int i = 0; i < pR->getSize(); i++)
95             {
96                 pDblCoef = pR->get(i)->getCoef();
97                 pReal = pDblCoef->getReal();
98
99                 for (int j = 0; j < pDblCoef->getSize(); j++)
100                 {
101                     pReal[j] *= -1;
102                 }
103             }
104         }
105
106         return pR;
107     }
108
109     // - Int
110     if (_pRightOperand->isInt())
111     {
112         switch (_pRightOperand->getType())
113         {
114             case InternalType::ScilabInt8 :
115             {
116                 Int8 *pR = _pRightOperand->clone()->getAs<Int8>();
117                 char* pReal = pR->get();
118
119                 for (int i = 0 ; i < pR->getSize() ; i++)
120                 {
121                     pReal[i] *= -1;
122                 }
123                 return pR;
124             }
125             case InternalType::ScilabInt16 :
126             {
127                 Int16 *pR = _pRightOperand->clone()->getAs<Int16>();
128                 short* pReal = pR->get();
129
130                 for (int i = 0 ; i < pR->getSize() ; i++)
131                 {
132                     pReal[i] *= -1;
133                 }
134                 return pR;
135             }
136             case InternalType::ScilabInt32 :
137             {
138                 Int32 *pR = _pRightOperand->clone()->getAs<Int32>();
139                 int* pReal = pR->get();
140
141                 for (int i = 0 ; i < pR->getSize() ; i++)
142                 {
143                     pReal[i] *= -1;
144                 }
145                 return pR;
146             }
147             case InternalType::ScilabInt64 :
148             {
149                 Int64 *pR = _pRightOperand->clone()->getAs<Int64>();
150                 long long* pReal = pR->get();
151
152                 for (int i = 0 ; i < pR->getSize() ; i++)
153                 {
154                     pReal[i] *= -1;
155                 }
156                 return pR;
157             }
158             case InternalType::ScilabUInt8 :
159             {
160                 UInt8 *pR = _pRightOperand->clone()->getAs<UInt8>();
161                 unsigned char* pReal = pR->get();
162
163                 for (int i = 0 ; i < pR->getSize() ; i++)
164                 {
165                     pReal[i] *= -1;
166                 }
167                 return pR;
168             }
169             case InternalType::ScilabUInt16 :
170             {
171                 UInt16 *pR = _pRightOperand->clone()->getAs<UInt16>();
172                 unsigned short* pReal = pR->get();
173
174                 for (int i = 0 ; i < pR->getSize() ; i++)
175                 {
176                     pReal[i] *= -1;
177                 }
178                 return pR;
179             }
180             case InternalType::ScilabUInt32 :
181             {
182                 UInt32 *pR = _pRightOperand->clone()->getAs<UInt32>();
183                 unsigned int* pReal = pR->get();
184
185                 for (int i = 0 ; i < pR->getSize() ; i++)
186                 {
187                     pReal[i] *= -1;
188                 }
189                 return pR;
190             }
191             case InternalType::ScilabUInt64 :
192             {
193                 UInt64 *pR = _pRightOperand->clone()->getAs<UInt64>();
194                 unsigned long long* pReal = pR->get();
195
196                 for (int i = 0 ; i < pR->getSize() ; i++)
197                 {
198                     pReal[i] *= -1;
199                 }
200                 return pR;
201             }
202             default :
203                 return NULL;
204         }
205     }
206
207     // - Sparse
208     if (_pRightOperand->isSparse())
209     {
210         Sparse* pSp = dynamic_cast<Sparse*>(_pRightOperand)->multiply(-1);
211         return pSp;
212     }
213
214     /*
215     ** Default case : Return NULL will Call Overloading.
216     */
217     return NULL;
218 }
219
220 InternalType* GenericMinus(InternalType* _pLeftOperand, InternalType* _pRightOperand)
221 {
222     InternalType *pResult = NULL;
223
224     // InternalType - [] and [] - []
225     if (_pRightOperand->isDouble() && _pRightOperand->getAs<Double>()->isEmpty())
226     {
227         return _pLeftOperand->clone();
228     }
229
230     // [] - InternalType
231     if (_pLeftOperand->isDouble() && _pLeftOperand->getAs<Double>()->isEmpty())
232     {
233         if (_pRightOperand->isSparse())
234         {
235             InternalType* pSpOut = NULL;
236             Sparse* pSp = new Sparse(1, 1);
237             pSp->set(0, 0, -1.0);
238             pSpOut = GenericTimes((InternalType*)pSp, _pRightOperand);
239             delete pSp;
240             return pSpOut;
241         }
242
243         Double pDblZero = Double(0);
244         return GenericMinus((InternalType*)&pDblZero, _pRightOperand);
245     }
246
247     /*
248     ** DOUBLE - DOUBLE
249     */
250     if (_pLeftOperand->isDouble() && _pRightOperand->isDouble())
251     {
252         Double *pL = _pLeftOperand->getAs<Double>();
253         Double *pR = _pRightOperand->getAs<Double>();
254
255         int iResult = SubstractDoubleToDouble(pL, pR, (Double**)&pResult);
256         if (iResult != 0)
257         {
258             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
259         }
260         return pResult;
261     }
262
263     /*
264     ** DOUBLE - POLY
265     */
266     else if (_pLeftOperand->isDouble() && _pRightOperand->isPoly())
267     {
268         Double *pL              = _pLeftOperand->getAs<Double>();
269         Polynom *pR          = _pRightOperand->getAs<types::Polynom>();
270
271         int iResult = SubstractPolyToDouble(pL, pR, (Polynom**)&pResult);
272         if (iResult != 0)
273         {
274             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
275         }
276         return pResult;
277     }
278
279     /*
280     ** POLY - DOUBLE
281     */
282     else if (_pLeftOperand->isPoly() && _pRightOperand->isDouble())
283     {
284         Polynom *pL   = _pLeftOperand->getAs<types::Polynom>();
285         Double *pR    = _pRightOperand->getAs<Double>();
286
287         int iResult = SubstractDoubleToPoly(pL, pR, (Polynom**)&pResult);
288         if (iResult != 0)
289         {
290             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
291         }
292         return pResult;
293     }
294
295     /*
296     ** POLY - POLY
297     */
298     else if (_pLeftOperand->isPoly() && _pRightOperand->isPoly())
299     {
300         Polynom *pL   = _pLeftOperand->getAs<types::Polynom>();
301         Polynom *pR   = _pRightOperand->getAs<types::Polynom>();
302
303         int iResult = SubstractPolyToPoly(pL, pR, (Polynom**)&pResult);
304         if (iResult != 0)
305         {
306             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
307         }
308         return pResult;
309     }
310
311     /*
312     ** SPARSE - SPARSE
313     */
314     else if (_pLeftOperand->isSparse() && _pRightOperand->isSparse())
315     {
316         Sparse *pL = _pLeftOperand->getAs<Sparse>();
317         Sparse *pR = _pRightOperand->getAs<Sparse>();
318
319         int iResult = SubstractSparseToSparse(pR, pL, (GenericType**)&pResult);
320         if (iResult != 0)
321         {
322             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
323         }
324
325         return pResult;
326     }
327
328     /*
329     ** SPARSE - DOUBLE
330     */
331     else if (_pLeftOperand->isSparse() && _pRightOperand->isDouble())
332     {
333         Sparse *pL = _pLeftOperand->getAs<Sparse>();
334         Double *pR = _pRightOperand->getAs<Double>();
335
336         int iResult = SubstractDoubleToSparse(pR, pL, (GenericType**)&pResult);
337         if (iResult != 0)
338         {
339             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
340         }
341
342         return pResult;
343     }
344
345     /*
346     ** DOUBLE - SPARSE
347     */
348     else if (_pLeftOperand->isDouble() && _pRightOperand->isSparse())
349     {
350         Double *pL = _pLeftOperand->getAs<Double>();
351         Sparse *pR = _pRightOperand->getAs<Sparse>();
352
353         int iResult = SubstractSparseToDouble(pR, pL, (GenericType**)&pResult);
354         if (iResult != 0)
355         {
356             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
357         }
358
359         return pResult;
360     }
361     /*
362     ** Default case : Return NULL will Call Overloading.
363     */
364     return NULL;
365 }
366
367 int SubstractDoubleToDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
368 {
369     bool bComplex1  = _pDouble1->isComplex();
370     bool bComplex2  = _pDouble2->isComplex();
371
372     if (_pDouble1->isIdentity())
373     {
374         if (_pDouble2->getDims() > 2)
375         {
376             //unable to substract identity matrix and greater than 2 dimensions matrix
377             return 1;
378         }
379
380         (*_pDoubleOut) = new Double(_pDouble2->getRows(), _pDouble2->getCols(), bComplex1 || bComplex2);
381
382         if (bComplex1 == false && bComplex2 == false)
383         {
384             iSubstractRealIdentityToRealMatrix(
385                 _pDouble1->get(0),
386                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
387                 (*_pDoubleOut)->get());
388         }
389         else if (bComplex1 == false && bComplex2 == true)
390         {
391             iSubstractRealIdentityToComplexMatrix(
392                 _pDouble1->get(0),
393                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
394                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
395         }
396         else if (bComplex1 == true && bComplex2 == false)
397         {
398             iSubstractComplexIdentityToRealMatrix(
399                 _pDouble1->get(0), _pDouble1->getImg(0),
400                 _pDouble2->get(), _pDouble2->getRows(), _pDouble2->getCols(),
401                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
402         }
403         else if (bComplex1 == true && bComplex2 == true)
404         {
405             iSubstractComplexIdentityToComplexMatrix(
406                 _pDouble1->get(0), _pDouble1->getImg(0),
407                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
408                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
409         }
410
411         //return opposite
412         double* pDbl = (*_pDoubleOut)->get();
413         if ((*_pDoubleOut)->isComplex())
414         {
415             double* pDblImg = (*_pDoubleOut)->getImg();
416             for (int i = 0 ; i < _pDouble2->getSize() ; i++)
417             {
418                 pDbl[i] = -pDbl[i];
419                 pDblImg[i] = -pDblImg[i];
420             }
421         }
422         else
423         {
424             for (int i = 0 ; i < _pDouble2->getSize() ; i++)
425             {
426                 pDbl[i] = -pDbl[i];
427             }
428         }
429         return 0;
430     }
431
432     if (_pDouble2->isIdentity())
433     {
434         if (_pDouble1->getDims() > 2)
435         {
436             //unable to substract identity matrix and greater than 2 dimensions matrix
437             return 1;
438         }
439
440         (*_pDoubleOut) = new Double(_pDouble1->getRows(), _pDouble1->getCols(), bComplex1 || bComplex2);
441         if (bComplex1 == false && bComplex2 == false)
442         {
443             iSubstractRealIdentityToRealMatrix(
444                 _pDouble2->get(0),
445                 _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
446                 (*_pDoubleOut)->get());
447         }
448         else if (bComplex1 == false && bComplex2 == true)
449         {
450             iSubstractComplexIdentityToRealMatrix(
451                 _pDouble2->get(0), _pDouble2->getImg(0),
452                 _pDouble1->get(), _pDouble1->getRows(), _pDouble1->getCols(),
453                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
454         }
455         else if (bComplex1 == true && bComplex2 == false)
456         {
457             iSubstractRealIdentityToComplexMatrix(
458                 _pDouble2->get(0),
459                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
460                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
461         }
462         else if (bComplex1 == true && bComplex2 == true)
463         {
464             iSubstractComplexIdentityToComplexMatrix(
465                 _pDouble2->get(0), _pDouble2->getImg(0),
466                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
467                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
468         }
469
470         return 0;
471     }
472
473     if (_pDouble1->isScalar())
474     {
475         //add pL with each element of pR
476         (*_pDoubleOut) = new Double(_pDouble2->getDims(), _pDouble2->getDimsArray(), bComplex1 || bComplex2);
477         if (bComplex1 == false && bComplex2 == false)
478         {
479             iSubstractRealMatrixToRealScalar(
480                 _pDouble2->get(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
481                 _pDouble1->get(0),
482                 (*_pDoubleOut)->get());
483         }
484         else if (bComplex1 == false && bComplex2 == true)
485         {
486             iSubstractComplexMatrixToRealScalar(
487                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
488                 _pDouble1->get(0),
489                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
490         }
491         else if (bComplex1 == true && bComplex2 == false)
492         {
493             iSubstractRealMatrixToComplexScalar(
494                 _pDouble2->get(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
495                 _pDouble1->get(0), _pDouble1->getImg(0),
496                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
497         }
498         else if (bComplex1 == true && bComplex2 == true)
499         {
500             iSubstractComplexMatrixToComplexScalar(
501                 _pDouble2->get(), _pDouble2->getImg(), _pDouble2->getDimsArray(), _pDouble2->getDims(),
502                 _pDouble1->get(0), _pDouble1->getImg(0),
503                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
504         }
505
506         return 0;
507     }
508
509     if (_pDouble2->isScalar())
510     {
511         //add pL with each element of pR
512         (*_pDoubleOut) = new Double(_pDouble1->getDims(), _pDouble1->getDimsArray(), bComplex1 || bComplex2);
513         if (bComplex1 == false && bComplex2 == false)
514         {
515             iSubstractRealScalarToRealMatrix(
516                 _pDouble2->get(0),
517                 _pDouble1->get(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
518                 (*_pDoubleOut)->get());
519         }
520         else if (bComplex1 == false && bComplex2 == true)
521         {
522             iSubstractComplexScalarToRealMatrix(
523                 _pDouble2->get(0), _pDouble2->getImg(0),
524                 _pDouble1->get(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
525                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
526         }
527         else if (bComplex1 == true && bComplex2 == false)
528         {
529             iSubstractRealScalarToComplexMatrix(
530                 _pDouble2->get(0),
531                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
532                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
533         }
534         else if (bComplex1 == true && bComplex2 == true)
535         {
536             iSubstractComplexScalarToComplexMatrix(
537                 _pDouble2->get(0), _pDouble2->getImg(0),
538                 _pDouble1->get(), _pDouble1->getImg(), _pDouble1->getDimsArray(), _pDouble1->getDims(),
539                 (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
540         }
541
542         return 0;
543     }
544
545
546     //add pL and pR element wise
547
548     int iDims1 = _pDouble1->getDims();
549     int iDims2 = _pDouble2->getDims();
550
551     if (iDims1 != iDims2)
552     {
553         return 1;
554     }
555
556     int* piDims1    = _pDouble1->getDimsArray();
557     int* piDims2    = _pDouble2->getDimsArray();
558     for (int i = 0 ; i < _pDouble1->getDims() ; i++)
559     {
560         if (piDims1[i] != piDims2[i])
561         {
562             return 1;
563         }
564     }
565
566     (*_pDoubleOut) = new Double(iDims2, piDims2, bComplex1 || bComplex2);
567     if (bComplex1 == false && bComplex2 == false)
568     {
569         iSubstractRealMatrixToRealMatrix(
570             _pDouble2->get(),
571             _pDouble1->get(), piDims1, iDims1,
572             (*_pDoubleOut)->get());
573     }
574     else if (bComplex1 == false && bComplex2 == true)
575     {
576         iSubstractComplexMatrixToRealMatrix(
577             _pDouble2->get(), _pDouble2->getImg(),
578             _pDouble1->get(), piDims1, iDims1,
579             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
580     }
581     else if (bComplex1 == true && bComplex2 == false)
582     {
583         iSubstractRealMatrixToComplexMatrix(
584             _pDouble2->get(),
585             _pDouble1->get(), _pDouble1->getImg(), piDims1, iDims1,
586             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
587     }
588     else if (bComplex1 == true && bComplex2 == true)
589     {
590         iSubstractComplexMatrixToComplexMatrix(
591             _pDouble2->get(), _pDouble2->getImg(),
592             _pDouble1->get(), _pDouble1->getImg(), piDims1, iDims1,
593             (*_pDoubleOut)->get(), (*_pDoubleOut)->getImg());
594     }
595
596     return 0;
597 }
598
599 //1 - %s
600 int SubstractPolyToDouble(Double *_pDouble, Polynom *_pPoly, Polynom** _pPolyOut)
601 {
602     double *pInDblR = _pDouble->get();
603     double *pInDblI = _pDouble->getImg();
604
605     if (_pDouble->isScalar())
606     {
607         //Clone original poly
608         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
609
610         if (_pDouble->isComplex())
611         {
612             (*_pPolyOut)->setComplex(true);
613         }
614
615         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
616         {
617             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
618             double *pOutPolyR       = pOutPoly->getCoef()->get();
619
620             pOutPolyR[0] = pInDblR[0] - pOutPolyR[0];
621
622             for (int j = 1 ; j < pOutPoly->getRank() ; j++)
623             {
624                 pOutPolyR[j] = - pOutPolyR[j];
625             }
626         }
627
628         if ((*_pPolyOut)->isComplex())
629         {
630             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
631             {
632                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
633                 double *pOutPolyI       = pOutPoly->getCoef()->getImg();
634
635                 pOutPolyI[0]            = (pInDblI == NULL ? 0 : pInDblI[0]) - (pOutPolyI == NULL ? 0 : pOutPolyI[0]);
636
637                 for (int j = 1 ; j < pOutPoly->getRank() ; j++)
638                 {
639                     pOutPolyI[j] = - pOutPolyI[j];
640                 }
641             }
642         }
643
644         return 0;
645     }
646
647     if (_pPoly->isScalar())
648     {
649         //[] - %s
650         int *piRank = new int[_pDouble->getSize()];
651         for (int i = 0 ; i < _pDouble->getSize() ; i++)
652         {
653             piRank[i] = _pPoly->get(0)->getRank();
654         }
655
656         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
657
658         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
659         {
660             SinglePoly *pInPoly   = _pPoly->get(0);
661             SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
662
663             double *pInPolyR = pInPoly->getCoef()->get();
664             double *pOutPolyR = pOutPoly->getCoef()->get();
665
666             pOutPolyR[0] = pInDblR[i] - pInPolyR[0];
667
668             for (int j = 1 ; j < pOutPoly->getRank() ; j++)
669             {
670                 pOutPolyR[j] = - pInPolyR[j];
671             }
672         }
673
674         if (_pPoly->isComplex() || _pDouble->isComplex())
675         {
676             (*_pPolyOut)->setComplex(true);
677             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
678             {
679                 SinglePoly *pInPoly   = _pPoly->get(0);
680                 SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
681
682                 double *pInPolyI = pInPoly->getCoef()->getImg();
683                 double *pOutPolyI = pOutPoly->getCoef()->getImg();
684
685                 pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pInPolyI != NULL ? pInPolyI[0] : 0);
686
687                 for (int j = 1 ; j < pOutPoly->getRank() ; j++)
688                 {
689                     pOutPolyI[j] = -pInPolyI[j];
690                 }
691             }
692         }
693
694         return 0;
695     }
696
697     //check dimension compatibilities
698     int iLeftDims = _pDouble->getDims();
699     int iRightDims = _pPoly->getDims();
700
701     if (iLeftDims != iRightDims)
702     {
703         return 1;
704     }
705
706     int* piLeftDims = _pDouble->getDimsArray();
707     int* piRightDims = _pPoly->getDimsArray();
708
709     for (int i = 0 ; i < iLeftDims ; i++)
710     {
711         if (piLeftDims[i] != piRightDims[i])
712         {
713             return 1;
714         }
715     }
716
717
718     //Clone original poly
719     (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
720
721     for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
722     {
723         SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
724         double *pOutPolyR       = pOutPoly->getCoef()->get();
725
726         pOutPolyR[0] = pInDblR[i] - pOutPolyR[0];
727
728         for (int j = 1 ; j < pOutPoly->getRank() ; j++)
729         {
730             pOutPolyR[j] = - pOutPolyR[j];
731         }
732     }
733
734     if (_pPoly->isComplex() || _pDouble->isComplex())
735     {
736         (*_pPolyOut)->setComplex(true);
737         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
738         {
739             SinglePoly *pOutPoly  = (*_pPolyOut)->get(i);
740             double *pOutPolyI = pOutPoly->getCoef()->getImg();
741
742             pOutPolyI[0] = (pInDblI != NULL ? pInDblI[i] : 0) - (pOutPolyI != NULL ? pOutPolyI[0] : 0);
743
744             for (int j = 1 ; j < pOutPoly->getRank() ; j++)
745             {
746                 pOutPolyI[j] = - pOutPolyI[j];
747             }
748         }
749     }
750
751     return 0;
752 }
753
754 //%s - 1
755 int SubstractDoubleToPoly(Polynom *_pPoly, Double *_pDouble, Polynom **_pPolyOut)
756 {
757     double *pInDblR   = _pDouble->get();
758     double *pInDblI   = _pDouble->getImg();
759
760     if (_pDouble->isScalar())
761     {
762         //Clone original poly
763         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
764         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
765         {
766             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
767             double *pOutPolyR       = pOutPoly->getCoef()->get();
768             pOutPolyR[0]   -= pInDblR[0];
769         }
770
771         if ((*_pPolyOut)->isComplex() || _pDouble->isComplex())
772         {
773             (*_pPolyOut)->setComplex(true);
774             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
775             {
776                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
777                 double *pOutPolyI       = pOutPoly->getCoef()->getImg();
778
779                 pOutPolyI[0]            -=  (pInDblI == NULL ? 0 : pInDblI[0]);
780             }
781         }
782     }
783     else if (_pPoly->isScalar())
784     {
785         int *piRank = new int[_pDouble->getSize()];
786         for (int i = 0 ; i < _pDouble->getSize() ; i++)
787         {
788             piRank[i] = _pPoly->get(0)->getRank();
789         }
790
791         (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), _pDouble->getDims(), _pDouble->getDimsArray(), piRank);
792
793         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
794         {
795             SinglePoly *pInPoly     = _pPoly->get(0);
796             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
797             double *pInPolyR        = pInPoly->getCoef()->get();
798             double *pOutPolyR       = pOutPoly->getCoef()->get();
799
800             pOutPolyR[0]            = pInPolyR[0] - pInDblR[i];
801
802             for (int j = 1 ; j < pOutPoly->getRank() ; j++)
803             {
804                 pOutPolyR[j]        = pInPolyR[j];
805             }
806         }
807
808         if (_pPoly->isComplex() || _pDouble->isComplex())
809         {
810             (*_pPolyOut)->setComplex(true);
811             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
812             {
813                 SinglePoly *pInPoly     = _pPoly->get(0);
814                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
815                 double *pInPolyI        = pInPoly->getCoef()->getImg();
816                 double *pOutPolyI       = pOutPoly->getCoef()->getImg();
817
818                 pOutPolyI[0] = (pInPolyI != NULL ? pInPolyI[0] : 0) - (pInDblI != NULL ? pInDblI[i] : 0);
819
820                 for (int j = 1 ; j < pOutPoly->getRank() ; j++)
821                 {
822                     pOutPolyI[j] = pInPolyI[j];
823                 }
824             }
825         }
826     }
827     else
828     {
829         //check dimension compatibilities
830         int iLeftDims = _pDouble->getDims();
831         int iRightDims = _pPoly->getDims();
832
833         if (iLeftDims != iRightDims)
834         {
835             return 1;
836         }
837
838         int* piLeftDims = _pDouble->getDimsArray();
839         int* piRightDims = _pPoly->getDimsArray();
840
841         for (int i = 0 ; i < iLeftDims ; i++)
842         {
843             if (piLeftDims[i] != piRightDims[i])
844             {
845                 return 1;
846             }
847         }
848
849         //Clone original poly
850         (*_pPolyOut) = _pPoly->clone()->getAs<Polynom>();
851         for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
852         {
853             SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
854             double *pOutPolyR       = pOutPoly->getCoef()->get();
855             pOutPolyR[0]            -= pInDblR[i];
856         }
857
858         if (_pPoly->isComplex() || _pDouble->isComplex())
859         {
860             (*_pPolyOut)->setComplex(true);
861             for (int i = 0 ; i < (*_pPolyOut)->getSize() ; i++)
862             {
863                 SinglePoly *pOutPoly    = (*_pPolyOut)->get(i);
864                 double *pOutPolyI       = pOutPoly->getCoef()->getImg();
865
866                 pOutPolyI[0]            -= (pInDblI != NULL ? pInDblI[i] : 0);
867             }
868         }
869     }
870     return 0;
871 }
872
873 int SubstractPolyToPoly(Polynom *_pPoly1, Polynom *_pPoly2, Polynom **_pPolyOut)
874 {
875     if (_pPoly1->isScalar())
876     {
877         int* pRankOut   = new int[_pPoly2->getSize()];
878         int* pRank1     = new int[_pPoly1->getSize()];
879         int* pRank2     = new int[_pPoly2->getSize()];
880         memset(pRank1, 0x00, _pPoly2->getSize() * sizeof(int));
881
882         _pPoly1->getRank(pRank1);
883         _pPoly2->getRank(pRank2);
884         for (int i = 0 ; i < _pPoly2->getSize() ; i++)
885         {
886             pRankOut[i] = std::max(pRank1[0], pRank2[i]);
887         }
888
889         (*_pPolyOut) = new Polynom(_pPoly2->getVariableName(), _pPoly2->getDims(), _pPoly2->getDimsArray(), pRankOut);
890         if (_pPoly1->isComplex() || _pPoly2->isComplex())
891         {
892             (*_pPolyOut)->setComplex(true);
893         }
894
895         //Result P1(0) + P2(i)
896         Double* p1Coef          = _pPoly1->get(0)->getCoef();
897         double* p1R             = p1Coef->get();
898
899         for (int i = 0 ; i < _pPoly2->getSize() ; i++)
900         {
901             Double* p2Coef      = _pPoly2->get(i)->getCoef();
902             double* p2R         = p2Coef->get();
903
904             Double* pOutCoef    = (*_pPolyOut)->get(i)->getCoef();
905             double* pOutR       = pOutCoef->get();
906
907             for (int j = 0 ; j < pRankOut[i] ; j++)
908             {
909                 if (j >= pRank1[0])
910                 {
911                     pOutR[j] = - p2R[j];
912                 }
913                 else if (j >= pRank2[i])
914                 {
915                     pOutR[j] = p1R[j];
916                 }
917                 else
918                 {
919                     pOutR[j] = p1R[j] - p2R[j];
920                 }
921             }
922
923             if ((*_pPolyOut)->isComplex())
924             {
925                 double *p1I     = p1Coef->getImg();
926                 double *p2I     = p2Coef->getImg();
927                 double *pOutI   = pOutCoef->getImg();
928
929                 for (int j = 0 ; j < pRankOut[i] ; j++)
930                 {
931                     pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
932                 }
933             }
934         }
935
936         delete[] pRankOut;
937         delete[] pRank1;
938         delete[] pRank2;
939         return 0;
940     }
941
942     if (_pPoly2->isScalar())
943     {
944         //size(p2) == 1
945         int *pRankOut   = new int[_pPoly1->getSize()];
946         int *pRank1     = new int[_pPoly1->getSize()];
947         int *pRank2     = new int[_pPoly2->getSize()];
948         memset(pRank2, 0x00, _pPoly1->getSize() * sizeof(int));
949
950         _pPoly1->getRank(pRank1);
951         _pPoly2->getRank(pRank2);
952         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
953         {
954             pRankOut[i] = std::max(pRank1[i], pRank2[0]);
955         }
956
957         (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), _pPoly1->getDims(), _pPoly1->getDimsArray(), pRankOut);
958         if (_pPoly1->isComplex() || _pPoly2->isComplex())
959         {
960             (*_pPolyOut)->setComplex(true);
961         }
962
963         //Result P1(i) + P2(0)
964         Double *p2Coef          = _pPoly2->get(0)->getCoef();
965         double *p2R             = p2Coef->get();
966
967         for (int i = 0 ; i < _pPoly1->getSize() ; i++)
968         {
969             Double *p1Coef      = _pPoly1->get(i)->getCoef();
970             double *p1R   = p1Coef->get();
971
972             Double *pOutCoef    = (*_pPolyOut)->get(i)->getCoef();
973             double *pOutR       = pOutCoef->get();
974
975             for (int j = 0 ; j < pRankOut[i] ; j++)
976             {
977                 if (j >= pRank1[j])
978                 {
979                     pOutR[j] = - p2R[j];
980                 }
981                 else if (j >= pRank2[0])
982                 {
983                     pOutR[j] = p1R[j];
984                 }
985                 else
986                 {
987                     pOutR[j] = p1R[j] - p2R[j];
988                 }
989             }
990
991             if ((*_pPolyOut)->isComplex())
992             {
993                 double *p2I  = p2Coef->getImg();
994                 double *p1I     = p1Coef->getImg();
995                 double *pOutI   = pOutCoef->getImg();
996                 for (int j = 0 ; j < pRankOut[i] ; j++)
997                 {
998                     pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
999                 }
1000             }
1001         }
1002
1003         delete[] pRankOut;
1004         delete[] pRank1;
1005         delete[] pRank2;
1006         return 0;
1007     }
1008
1009     //check dimension compatibilities
1010     int iDims1  = _pPoly1->getDims();
1011     int iDims2  = _pPoly2->getDims();
1012
1013     if (iDims1 != iDims2)
1014     {
1015         return 1;
1016     }
1017
1018     int* piDims1    = _pPoly1->getDimsArray();
1019     int* piDims2    = _pPoly2->getDimsArray();
1020
1021     for (int i = 0 ; i < iDims1 ; i++)
1022     {
1023         if (piDims1[i] != piDims2[i])
1024         {
1025             return 1;
1026         }
1027     }
1028
1029     int *pRankOut   = new int[_pPoly1->getSize()];
1030     int *pRank1     = new int[_pPoly1->getSize()];
1031     int *pRank2     = new int[_pPoly2->getSize()];
1032
1033     _pPoly1->getRank(pRank1);
1034     _pPoly2->getRank(pRank2);
1035     for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1036     {
1037         pRankOut[i] = std::max(pRank1[i], pRank2[i]);
1038     }
1039
1040     (*_pPolyOut) = new Polynom(_pPoly2->getVariableName(), iDims1, piDims1, pRankOut);
1041     if (_pPoly1->isComplex() || _pPoly2->isComplex())
1042     {
1043         (*_pPolyOut)->setComplex(true);
1044     }
1045
1046     //Result P1(i) + P2(i)
1047     for (int i = 0 ; i < _pPoly1->getSize() ; i++)
1048     {
1049         double *p1R     = _pPoly1->get(i)->getCoef()->get();
1050         double *p2R     = _pPoly2->get(i)->getCoef()->get();
1051         double *pOutR   = (*_pPolyOut)->get(i)->getCoef()->get();
1052
1053         for (int j = 0 ; j < std::min(pRank1[i], pRank2[i]) ; j++)
1054         {
1055             pOutR[j]    = p1R[j] - p2R[j];
1056         }
1057
1058         double *pTemp   = NULL;
1059         int iCoef       = 1;
1060         if (pRank1[i] > pRank2[i])
1061         {
1062             pTemp       = p1R;
1063             iCoef       = 1;
1064         }
1065         else
1066         {
1067             pTemp       = p2R;
1068             iCoef       = -1;
1069         }
1070
1071         for (int j = std::min(pRank1[i], pRank2[i]) ; j < std::max(pRank1[i], pRank2[i]) ; j++)
1072         {
1073             pOutR[j]    = pTemp[j] * iCoef;
1074         }
1075
1076         if ((*_pPolyOut)->isComplex())
1077         {
1078             double *p1I     = _pPoly1->get(i)->getCoef()->getImg();
1079             double *p2I     = _pPoly2->get(i)->getCoef()->getImg();
1080             double *pOutI   = (*_pPolyOut)->get(i)->getCoef()->getImg();
1081
1082             for (int j = 0 ; j < std::min(pRank1[i], pRank2[i]) ; j++)
1083             {
1084                 pOutI[j]    = (p1I == NULL ? 0 : p1I[j]) - (p2I == NULL ? 0 : p2I[j]);
1085             }
1086
1087             for (int j = std::min(pRank1[i], pRank2[i]) ; j < std::max(pRank1[i], pRank2[i]) ; j++)
1088             {
1089                 pOutI[j]  = pTemp[j] * iCoef;
1090             }
1091         }
1092     }
1093
1094     delete[] pRankOut;
1095     delete[] pRank1;
1096     delete[] pRank2;
1097
1098     if ((*_pPolyOut) != NULL)
1099     {
1100         (*_pPolyOut)->updateRank();
1101     }
1102
1103     return 0;
1104 }
1105
1106 //_pSparse2 - _pSparse1
1107 int SubstractSparseToSparse(Sparse* _pSparse1, Sparse* _pSparse2, GenericType **_pSparseOut)
1108 {
1109     //check scalar hidden in a sparse ;)
1110     if (_pSparse1->getRows() == 1 && _pSparse1->getCols() == 1)
1111     {
1112         //do scalar - sp
1113         Double* pDbl = NULL;
1114         if (_pSparse1->isComplex())
1115         {
1116             std::complex<double> dbl = _pSparse1->getImg(0, 0);
1117             pDbl = new Double(dbl.real(), dbl.imag());
1118         }
1119         else
1120         {
1121             pDbl = new Double(_pSparse1->get(0, 0));
1122         }
1123
1124         SubstractSparseToDouble(_pSparse2, pDbl, (GenericType**)_pSparseOut);
1125         delete pDbl;
1126         return 0;
1127     }
1128
1129     if (_pSparse2->getRows() == 1 && _pSparse2->getCols() == 1)
1130     {
1131         //do sp - scalar
1132         Double* pDbl = NULL;
1133         if (_pSparse2->isComplex())
1134         {
1135             std::complex<double> dbl = _pSparse2->getImg(0, 0);
1136             pDbl = new Double(dbl.real(), dbl.imag());
1137         }
1138         else
1139         {
1140             pDbl = new Double(_pSparse2->get(0, 0));
1141         }
1142
1143         SubstractDoubleToSparse(pDbl, _pSparse1, (GenericType**)_pSparseOut);
1144         delete pDbl;
1145         return 0;
1146     }
1147
1148     if (_pSparse1->getRows() != _pSparse2->getRows() || _pSparse1->getCols() != _pSparse2->getCols())
1149     {
1150         //dimensions not match
1151         return 1;
1152     }
1153
1154     if (_pSparse1->nonZeros() == 0)
1155     {
1156         //sp - sp([])
1157         *_pSparseOut = new Sparse(*_pSparse2);
1158         return 0;
1159     }
1160
1161     if (_pSparse2->nonZeros() == 0)
1162     {
1163         //sp([]) - sp
1164         Sparse* pOut = new Sparse(*_pSparse1);
1165         pOut->opposite();
1166         *_pSparseOut = pOut;
1167         return 0;
1168     }
1169
1170     //copy _pSparse1 in _pSparseOut
1171     *_pSparseOut = _pSparse2->substract(*_pSparse1);
1172     return 0;
1173 }
1174
1175 int SubstractSparseToDouble(Sparse* _pSparse, Double* _pDouble, GenericType **_pOut)
1176 {
1177     //D - SP
1178     int iOne = 1; //fortran
1179     bool bComplex1 = _pSparse->isComplex();
1180     bool bComplex2 = _pDouble->isComplex();
1181
1182     if (_pDouble->isIdentity())
1183     {
1184         //convert to sp
1185         Sparse* pS = new Sparse(_pSparse->getRows(), _pSparse->getCols(), _pDouble->isComplex());
1186         if (pS->isComplex())
1187         {
1188             for (int i = 0 ; i < std::min(_pSparse->getRows() , _pSparse->getCols()) ; i++)
1189             {
1190                 pS->set(i, i, std::complex<double>(_pDouble->get(0), _pDouble->getImg(0)));
1191             }
1192         }
1193         else
1194         {
1195             for (int i = 0 ; i < std::min(_pSparse->getRows() , _pSparse->getCols()) ; i++)
1196             {
1197                 pS->set(i, i, _pDouble->get(0));
1198             }
1199         }
1200
1201         SubstractSparseToSparse(_pSparse, pS, _pOut);
1202         delete pS;
1203         return 0;
1204     }
1205
1206     if (_pSparse->isScalar() && _pDouble->isScalar())
1207     {
1208         //d - sp
1209         Double* pRes = _pDouble->clone()->getAs<Double>();
1210         pRes->setComplex(bComplex1 | bComplex2);
1211         if (bComplex1)
1212         {
1213             std::complex<double> dbl = _pSparse->getImg(0, 0);
1214             pRes->set(0, pRes->get(0) - dbl.real());
1215             pRes->setImg(0, pRes->getImg(0) - dbl.imag());
1216         }
1217         else
1218         {
1219             pRes->set(0, pRes->get(0) - _pSparse->get(0, 0));
1220         }
1221
1222         *_pOut = pRes;
1223         return 0;
1224     }
1225
1226     if (_pDouble->isScalar())
1227     {
1228         //d - SP
1229         Double* pRes = new Double(_pSparse->getRows(), _pSparse->getCols(), bComplex1 | bComplex2);
1230         int iSize = _pSparse->getSize();
1231         double dblVal = _pDouble->get(0);
1232         C2F(dset)(&iSize, &dblVal, pRes->get(), &iOne);
1233         if (bComplex2)
1234         {
1235             double dblValI = _pDouble->getImg(0);
1236             C2F(dset)(&iSize, &dblValI, pRes->getImg(), &iOne);
1237         }
1238         else if (bComplex1)
1239         {
1240             //initialize imag part at 0
1241             double dblValI = 0;
1242             C2F(dset)(&iSize, &dblValI, pRes->getImg(), &iOne);
1243         }
1244
1245         int nonZeros = static_cast<int>(_pSparse->nonZeros());
1246         int* pRows = new int[nonZeros * 2];
1247         _pSparse->outputRowCol(pRows);
1248         int* pCols = pRows + nonZeros;
1249
1250         if (bComplex1)
1251         {
1252             for (int i = 0 ; i < nonZeros ; i++)
1253             {
1254                 int iRow = static_cast<int>(pRows[i]) - 1;
1255                 int iCol = static_cast<int>(pCols[i]) - 1;
1256                 std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1257                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - dbl.real());
1258                 pRes->setImg(iRow, iCol, pRes->getImg(iRow, iCol) - dbl.imag());
1259             }
1260         }
1261         else
1262         {
1263             for (int i = 0 ; i < nonZeros ; i++)
1264             {
1265                 int iRow = static_cast<int>(pRows[i]) - 1;
1266                 int iCol = static_cast<int>(pCols[i]) - 1;
1267                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - _pSparse->get(iRow, iCol));
1268             }
1269         }
1270         *_pOut = pRes;
1271
1272         //clear
1273         delete[] pRows;
1274         return 0;
1275     }
1276
1277     if (_pSparse->isScalar())
1278     {
1279         //D - sp
1280         Double* pRes = _pDouble->clone()->getAs<Double>();
1281         pRes->setComplex(bComplex1 | bComplex2);
1282
1283         if (bComplex1)
1284         {
1285             double* pReal = pRes->get();
1286             double* pImg = pRes->getImg();
1287             for (int i = 0 ; i < pRes->getSize() ; i++)
1288             {
1289                 std::complex<double> dbl = _pSparse->getImg(0, 0);
1290                 pReal[i] -= dbl.real();
1291                 pImg[i] -= dbl.imag();
1292             }
1293         }
1294         else
1295         {
1296             double* pReal = pRes->get();
1297             for (int i = 0 ; i < pRes->getSize() ; i++)
1298             {
1299                 pReal[i] -= _pSparse->get(0, 0);
1300             }
1301         }
1302         *_pOut = pRes;
1303         return 0;
1304     }
1305
1306
1307     if (_pSparse->getRows() == _pDouble->getRows() && _pSparse->getCols() == _pDouble->getCols())
1308     {
1309         //D - SP
1310         Double* pRes = _pDouble->clone()->getAs<Double>();
1311         pRes->setComplex(bComplex1 | bComplex2);
1312
1313         int nonZeros = static_cast<int>(_pSparse->nonZeros());
1314         int* pRows = new int[nonZeros * 2];
1315         _pSparse->outputRowCol(pRows);
1316         int* pCols = pRows + nonZeros;
1317
1318         if (bComplex1)
1319         {
1320             for (int i = 0 ; i < nonZeros ; i++)
1321             {
1322                 int iRow = static_cast<int>(pRows[i]) - 1;
1323                 int iCol = static_cast<int>(pCols[i]) - 1;
1324                 std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1325                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - dbl.real());
1326                 pRes->setImg(iRow, iCol, pRes->getImg(iRow, iCol) - dbl.imag());
1327             }
1328         }
1329         else
1330         {
1331             for (int i = 0 ; i < nonZeros ; i++)
1332             {
1333                 int iRow = static_cast<int>(pRows[i]) - 1;
1334                 int iCol = static_cast<int>(pCols[i]) - 1;
1335                 pRes->set(iRow, iCol, pRes->get(iRow, iCol) - _pSparse->get(iRow, iCol));
1336             }
1337         }
1338
1339         //clear
1340         delete[] pRows;
1341         *_pOut = pRes;
1342         return 0;
1343     }
1344     return 1;
1345 }
1346
1347 int SubstractDoubleToSparse(Double* _pDouble, Sparse* _pSparse, GenericType **_pOut)
1348 {
1349     //SP - D => -(D - SP)
1350     int iRet = SubstractSparseToDouble(_pSparse, _pDouble, _pOut);
1351     if (iRet)
1352     {
1353         return iRet;
1354     }
1355
1356     //compute opposite of result
1357     if ((*_pOut)->isDouble())
1358     {
1359         Double* pD          = (*_pOut)->getAs<Double>();
1360         int iSize           = pD->getSize();
1361         double dblMinusOne  = -1;
1362         int iOne            = 1;
1363
1364         C2F(dscal)(&iSize, &dblMinusOne, pD->get(), &iOne);
1365         if (pD->isComplex())
1366         {
1367             C2F(dscal)(&iSize, &dblMinusOne, pD->getImg(), &iOne);
1368         }
1369     }
1370     else if ((*_pOut)->isSparse())
1371     {
1372         Sparse* pS = (*_pOut)->getAs<Sparse>();
1373         pS->opposite();
1374     }
1375     else
1376     {
1377         return 1;
1378     }
1379
1380     return 0;
1381 }
1382