adding sparse management via Eigen library
[scilab.git] / scilab / modules / operations / src / cpp / types_multiplication.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 "types_multiplication.hxx"
15 #include "types_addition.hxx"
16 #include "arrayof.hxx"
17
18 #include "scilabexception.hxx"
19
20 extern "C"
21 {
22 #include "matrix_multiplication.h"
23 #include "matrix_addition.h"
24 #include "operation_f.h"
25 #include "localization.h"
26 #include "charEncoding.h"
27 }
28
29 using namespace types;
30
31 InternalType *GenericDotTimes(InternalType *_pLeftOperand, InternalType *_pRightOperand)
32 {
33     InternalType *pResult = NULL;
34
35     /*
36     ** DOUBLE .* DOUBLE
37     */
38     if(_pLeftOperand->isDouble() && _pRightOperand->isDouble())
39     {
40         Double *pL                      = _pLeftOperand->getAs<Double>();
41         Double *pR                      = _pRightOperand->getAs<Double>();
42
43         int iResult = DotMultiplyDoubleByDouble(pL, pR, (Double**)&pResult);
44         if(iResult)
45         {
46             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
47         }
48
49         return pResult;
50     }
51
52     /*
53     ** SPARSE .* SPARSE
54     */
55     if(_pLeftOperand->isSparse() && _pRightOperand->isSparse())
56     {
57         Sparse *pL                      = _pLeftOperand->getAs<Sparse>();
58         Sparse *pR                      = _pRightOperand->getAs<Sparse>();
59
60         int iResult = DotMultiplySparseBySparse(pL, pR, (Sparse**)&pResult);
61         if(iResult)
62         {
63             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
64         }
65
66         return pResult;
67     }
68
69     /*
70     ** SPARSE .* DOUBLE
71     */
72     if(_pLeftOperand->isSparse() && _pRightOperand->isDouble())
73     {
74         Sparse *pL                      = _pLeftOperand->getAs<Sparse>();
75         Double *pR                      = _pRightOperand->getAs<Double>();
76
77         int iResult = DotMultiplySparseByDouble(pL, pR, (GenericType**)&pResult);
78         if(iResult)
79         {
80             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
81         }
82
83         return pResult;
84     }
85
86     /*
87     ** DOUBLE .* SPARSE
88     */
89     if(_pLeftOperand->isDouble() && _pRightOperand->isSparse())
90     {
91         Double *pL                      = _pLeftOperand->getAs<Double>();
92         Sparse *pR                      = _pRightOperand->getAs<Sparse>();
93
94         int iResult = DotMultiplyDoubleBySparse(pL, pR, (GenericType**)&pResult);
95         if(iResult)
96         {
97             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
98         }
99
100         return pResult;
101     }
102
103     /*
104     ** Default case : Return NULL will Call Overloading.
105     */
106     return NULL;
107
108 }
109
110 InternalType *GenericTimes(InternalType *_pLeftOperand, InternalType *_pRightOperand)
111 {
112     InternalType *pResult = NULL;
113     GenericType::RealType TypeL = _pLeftOperand->getType();
114     GenericType::RealType TypeR = _pRightOperand->getType();
115
116     /*
117     ** DOUBLE * DOUBLE
118     */
119     if(TypeL == GenericType::RealDouble && TypeR == GenericType::RealDouble)
120     {
121         Double *pL                      = _pLeftOperand->getAs<Double>();
122         Double *pR                      = _pRightOperand->getAs<Double>();
123
124         int iResult = MultiplyDoubleByDouble(pL, pR, (Double**)&pResult);
125         if(iResult)
126         {
127             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
128         }
129
130         return pResult;
131     }
132
133     /*
134     ** DOUBLE * POLY
135     */
136     else if(TypeL == InternalType::RealDouble && TypeR == InternalType::RealPoly)
137     {
138         Double *pL                      = _pLeftOperand->getAs<Double>();
139         Polynom *pR         = _pRightOperand->getAs<types::Polynom>();
140
141         int iResult = MultiplyDoubleByPoly(pL, pR, (Polynom**)&pResult);
142         if(iResult)
143         {
144             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
145         }
146
147         return pResult;
148     }
149
150     /*
151     ** POLY * DOUBLE
152     */
153     else if(TypeL == InternalType::RealPoly && TypeR == InternalType::RealDouble)
154     {
155         Polynom *pL          = _pLeftOperand->getAs<types::Polynom>();
156         Double *pR              = _pRightOperand->getAs<Double>();
157
158         int iResult = MultiplyPolyByDouble(pL, pR, (Polynom**)&pResult);
159         if(iResult)
160         {
161             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
162         }
163
164         return pResult;
165     }
166
167     /*
168     ** POLY * POLY
169     */
170     else if(TypeL == InternalType::RealPoly && TypeR == InternalType::RealPoly)
171     {
172         Polynom *pL          = _pLeftOperand->getAs<types::Polynom>();
173         Polynom *pR          = _pRightOperand->getAs<types::Polynom>();
174
175         int iResult = MultiplyPolyByPoly(pL, pR, (Polynom**)&pResult);
176         if(iResult)
177         {
178             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
179         }
180
181         return pResult;
182     }
183
184     /*
185     ** SPARSE * SPARSE
186     */
187     if(TypeL == GenericType::RealSparse && TypeR == GenericType::RealSparse)
188     {
189         Sparse *pL                      = _pLeftOperand->getAs<Sparse>();
190         Sparse *pR                      = _pRightOperand->getAs<Sparse>();
191
192         int iResult = MultiplySparseBySparse(pL, pR, (Sparse**)&pResult);
193         if(iResult)
194         {
195             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
196         }
197
198         return pResult;
199     }
200
201     /*
202     ** DOUBLE * SPARSE
203     */
204     if(TypeL == GenericType::RealDouble && TypeR == GenericType::RealSparse)
205     {
206         Double *pL                      = _pLeftOperand->getAs<Double>();
207         Sparse *pR                      = _pRightOperand->getAs<Sparse>();
208
209         int iResult = MultiplyDoubleBySparse(pL, pR, (GenericType**)&pResult);
210         if(iResult)
211         {
212             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
213         }
214
215         return pResult;
216     }
217
218     /*
219     ** SPARSE * DOUBLE
220     */
221     if(TypeL == GenericType::RealSparse && TypeR == GenericType::RealDouble)
222     {
223         Sparse *pL                      = _pLeftOperand->getAs<Sparse>();
224         Double *pR                      = _pRightOperand->getAs<Double>();
225
226         int iResult = MultiplySparseByDouble(pL, pR, (GenericType**)&pResult);
227         if(iResult)
228         {
229             throw ast::ScilabError(_W("Inconsistent row/column dimensions.\n"));
230         }
231
232         return pResult;
233     }
234
235     /*
236     ** Default case : Return NULL will Call Overloading.
237     */
238     return NULL;
239
240 }
241
242 int MultiplyDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double** _pDoubleOut)
243 {
244     bool bComplex1  = _pDouble1->isComplex();
245     bool bComplex2  = _pDouble2->isComplex();
246     bool bScalar1   = _pDouble1->getRows() == 1 && _pDouble1->getCols() == 1;
247     bool bScalar2   = _pDouble2->getRows() == 1 && _pDouble2->getCols() == 1;
248
249     int iRowResult  = 0;
250     int iColResult  = 0;
251
252     if(bScalar1)
253     {
254         iRowResult = _pDouble2->getRows();
255         iColResult = _pDouble2->getCols();
256     }
257     else if(bScalar2)
258     {
259         iRowResult = _pDouble1->getRows();
260         iColResult = _pDouble1->getCols();
261     }
262     else if(_pDouble1->getCols() == _pDouble2->getRows())
263     {
264         iRowResult = _pDouble1->getRows();
265         iColResult = _pDouble2->getCols();
266     }
267     else
268     {
269         return 1;
270     }
271
272     //Output variables
273     bool bComplexOut    = bComplex1 || bComplex2;
274     (*_pDoubleOut)      = new Double(iRowResult, iColResult, bComplexOut);
275
276     double *pReal       = (*_pDoubleOut)->getReal();
277     double *pImg        = (*_pDoubleOut)->getImg();
278
279     if(bScalar1)
280     {//cst*b
281         if(bComplex1 == false && bComplex2 == false)
282         {
283             iMultiRealScalarByRealMatrix(_pDouble1->getReal(0,0), _pDouble2->getReal(), iRowResult, iColResult, pReal);
284         }
285         else if(bComplex1 == false && bComplex2 == true)
286         {
287             iMultiRealScalarByComplexMatrix(_pDouble1->getReal(0,0), _pDouble2->getReal(), _pDouble2->getImg(), iRowResult, iColResult, pReal, pImg);
288         }
289         else if(bComplex1 == true && bComplex2 == false)
290         {
291             iMultiComplexScalarByRealMatrix(_pDouble1->getReal(0,0), _pDouble1->getImg(0,0), _pDouble2->getReal(), iRowResult, iColResult, pReal, pImg);
292         }
293         else //if(bComplex1 == true && bComplex2 == true)
294         {
295             iMultiComplexScalarByComplexMatrix(_pDouble1->getReal(0,0), _pDouble1->getImg(0,0), _pDouble2->getReal(), _pDouble2->getImg(), iRowResult, iColResult, pReal, pImg);
296         }
297     }
298     else if(bScalar2)
299     {//a * cst
300         if(bComplex1 == false && bComplex2 == false)
301         {//Real Matrix by Real Scalar
302             iMultiRealScalarByRealMatrix(_pDouble2->getReal(0,0), _pDouble1->getReal(), iRowResult, iColResult, pReal);
303         }
304         else if(bComplex1 == false && bComplex2 == true)
305         {//Real Matrix by Scalar Complex
306             iMultiComplexScalarByRealMatrix(_pDouble2->getReal(0,0), _pDouble2->getImg(0,0), _pDouble1->getReal(), iRowResult, iColResult, pReal, pImg);
307         }
308         else if(bComplex1 == true && bComplex2 == false)
309         {
310             iMultiRealScalarByComplexMatrix(_pDouble2->getReal(0,0), _pDouble1->getReal(), _pDouble1->getImg(), iRowResult, iColResult, pReal, pImg);
311         }
312         else //if(bComplex1 == true && bComplex2 == true)
313         {
314             iMultiComplexScalarByComplexMatrix(_pDouble2->getReal(0,0), _pDouble2->getImg(0,0), _pDouble1->getReal(), _pDouble1->getImg(), iRowResult, iColResult, pReal, pImg);
315         }
316     }
317     else if(_pDouble1->getCols() == _pDouble2->getRows())
318     {//a * b
319         if(bComplex1 == false && bComplex2 == false)
320         {//Real Matrix by Real Matrix
321             iMultiRealMatrixByRealMatrix(
322                 _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
323                 _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
324                 pReal);
325         }
326         else if(bComplex1 == false && bComplex2 == true)
327         {//Real Matrix by Matrix Complex
328             iMultiRealMatrixByComplexMatrix(
329                 _pDouble1->getReal(), _pDouble1->getRows(), _pDouble1->getCols(),
330                 _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
331                 pReal, pImg);
332         }
333         else if(bComplex1 == true && bComplex2 == false)
334         {//Complex Matrix by Real Matrix
335             iMultiComplexMatrixByRealMatrix(
336                 _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
337                 _pDouble2->getReal(), _pDouble2->getRows(), _pDouble2->getCols(),
338                 pReal, pImg);
339         }
340         else //if(bComplex1 == true && bComplex2 == true)
341         {//Complex Matrix by Complex Matrix
342             iMultiComplexMatrixByComplexMatrix(
343                 _pDouble1->getReal(), _pDouble1->getImg(), _pDouble1->getRows(), _pDouble1->getCols(),
344                 _pDouble2->getReal(), _pDouble2->getImg(), _pDouble2->getRows(), _pDouble2->getCols(),
345                 pReal, pImg);
346         }
347     }
348     else
349     {
350         return 1;
351     }
352     return 0;
353 }
354
355 int DotMultiplyDoubleByDouble(Double* _pDouble1, Double* _pDouble2, Double**     _pDoubleOut)
356 {
357     bool bComplex1  = _pDouble1->isComplex();
358     bool bComplex2  = _pDouble2->isComplex();
359     bool bScalar1   = _pDouble1->getRows() == 1 && _pDouble1->getCols() == 1;
360     bool bScalar2   = _pDouble2->getRows() == 1 && _pDouble2->getCols() == 1;
361
362     int iRowResult  = 0;
363     int iColResult  = 0;
364
365     if(bScalar1)
366     {
367         iRowResult = _pDouble2->getRows();
368         iColResult = _pDouble2->getCols();
369     }
370     else if(bScalar2)
371     {
372         iRowResult = _pDouble1->getRows();
373         iColResult = _pDouble1->getCols();
374     }
375     else if(_pDouble1->getRows() == _pDouble2->getRows() && _pDouble1->getCols() == _pDouble2->getCols())
376     {
377         iRowResult = _pDouble2->getRows();
378         iColResult = _pDouble2->getCols();
379     }
380     else
381     {//Error
382         return 1;
383     }
384
385     bool bComplexOut    = bComplex1 || bComplex2;
386     (*_pDoubleOut)      = new Double(iRowResult, iColResult, bComplexOut);
387
388     double *pReal       = (*_pDoubleOut)->getReal();
389     double *pImg        = (*_pDoubleOut)->getImg();
390
391     if(bScalar1)
392     {
393         if(bComplex1 == false && bComplex2 == false)
394         {
395             iMultiRealScalarByRealMatrix(_pDouble1->getReal(0,0), _pDouble2->getReal(), iRowResult, iColResult, pReal);
396         }
397         else if(bComplex1 == false && bComplex2 == true)
398         {
399             iMultiRealScalarByComplexMatrix(_pDouble1->getReal(0,0), _pDouble2->getReal(), _pDouble2->getImg(), iRowResult, iColResult, pReal, pImg);
400         }
401         else if(bComplex1 == true && bComplex2 == false)
402         {
403             iMultiComplexScalarByRealMatrix(_pDouble1->getReal(0,0), _pDouble1->getImg(0,0), _pDouble2->getReal(), iRowResult, iColResult, pReal, pImg);
404         }
405         else //if(bComplex1 == true && bComplex2 == true)
406         {
407             iMultiComplexScalarByComplexMatrix(_pDouble1->getReal(0,0), _pDouble1->getImg(0,0), _pDouble2->getReal(), _pDouble2->getImg(), iRowResult, iColResult, pReal, pImg);
408         }
409     }
410     else if(bScalar2)
411     {
412         if(bComplex1 == false && bComplex2 == false)
413         {//Real Matrix by Real Scalar
414             iMultiRealScalarByRealMatrix(_pDouble2->getReal(0,0), _pDouble1->getReal(), iRowResult, iColResult, pReal);
415         }
416         else if(bComplex1 == false && bComplex2 == true)
417         {//Real Matrix by Scalar Complex
418             iMultiComplexScalarByRealMatrix(_pDouble2->getReal(0,0), _pDouble2->getImg(0,0), _pDouble1->getReal(), iRowResult, iColResult, pReal, pImg);
419         }
420         else if(bComplex1 == true && bComplex2 == false)
421         {
422             iMultiRealScalarByComplexMatrix(_pDouble2->getReal(0,0), _pDouble1->getReal(), _pDouble1->getImg(), iRowResult, iColResult, pReal, pImg);
423         }
424         else //if(bComplex1 == true && bComplex2 == true)
425         {
426             iMultiComplexScalarByComplexMatrix(_pDouble2->getReal(0,0), _pDouble2->getImg(0,0), _pDouble1->getReal(), _pDouble1->getImg(), iRowResult, iColResult, pReal, pImg);
427         }
428     }
429     else
430     {//multiply element wise element
431         if(bComplex1 == false && bComplex2 == false)
432         {
433             iDotMultiplyRealMatrixByRealMatrix(_pDouble1->getReal(), _pDouble2->getReal(), pReal, iRowResult, iColResult);
434         }
435         else if(bComplex1 == false && bComplex2 == true)
436         {
437             iDotMultiplyRealMatrixByComplexMatrix(_pDouble1->getReal(), _pDouble2->getReal(), _pDouble2->getImg(), pReal, pImg, iRowResult, iColResult);
438         }
439         else if(bComplex1 == true && bComplex2 == false)
440         {
441             iDotMultiplyComplexMatrixByRealMatrix(_pDouble1->getReal(), _pDouble1->getImg(), _pDouble2->getReal(), pReal, pImg, iRowResult, iColResult);
442         }
443         else //if(bComplex1 == true && bComplex2 == true)
444         {
445             iDotMultiplyComplexMatrixByComplexMatrix(_pDouble1->getReal(), _pDouble1->getImg(), _pDouble2->getReal(), _pDouble2->getImg(), pReal, pImg, iRowResult, iColResult);
446         }
447     }
448
449     return 0;
450 }
451 int MultiplyDoubleByPoly(Double* _pDouble, Polynom* _pPoly, Polynom** _pPolyOut)
452 {
453     bool bComplex1  = _pDouble->isComplex();
454     bool bComplex2  = _pPoly->isComplex();
455     bool bScalar1   = _pDouble->getRows() == 1 && _pDouble->getCols() == 1;
456     bool bScalar2   = _pPoly->getRows() == 1 && _pPoly->getCols() == 1;
457
458     int iRowResult  = 0;
459     int iColResult  = 0;
460     int *piRank     = NULL;
461
462     if(bScalar1)
463     {
464         iRowResult = _pPoly->getRows();
465         iColResult = _pPoly->getCols();
466
467         piRank = new int[iRowResult * iColResult];
468         for(int i = 0 ; i < iRowResult * iColResult ; i++)
469         {
470             piRank[i] = _pPoly->get(i)->getRank();
471         }
472     }
473     else if (bScalar2)
474     {
475         iRowResult = _pDouble->getRows();
476         iColResult = _pDouble->getCols();
477
478         piRank = new int[iRowResult * iColResult];
479         for(int i = 0 ; i < iRowResult * iColResult ; i++)
480         {
481             piRank[i] = _pPoly->get(0)->getRank();
482         }
483     }
484     else if(_pDouble->getCols() == _pPoly->getRows())
485     {
486         iRowResult = _pDouble->getRows();
487         iColResult = _pPoly->getCols();
488         piRank = new int[iRowResult * iColResult];
489         for(int i = 0 ; i < iRowResult * iColResult ; i++)
490         {
491             piRank[i] = _pPoly->getMaxRank();
492         }
493     }
494     else
495     {
496         return 1;
497     }
498
499     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
500     delete[] piRank;
501     if(bComplex1 || bComplex2)
502     {
503         (*_pPolyOut)->setComplex(true);
504     }
505
506     if(bScalar1)
507     {//cst * p
508         for(int i = 0 ; i < _pPoly->getSize() ; i++)
509         {
510             SinglePoly *pPolyIn     = _pPoly->get(i);
511             double* pRealIn         = pPolyIn->getCoef()->getReal();
512             double* pImgIn          = pPolyIn->getCoef()->getImg();
513
514             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
515             double* pRealOut        = pPolyOut->getCoef()->getReal();
516             double* pImgOut         = pPolyOut->getCoef()->getImg();
517
518             if(bComplex1 == false && bComplex2 == false)
519             {
520                 iMultiRealScalarByRealMatrix(
521                     _pDouble->getReal(0,0),
522                     pRealIn, 1, pPolyIn->getRank(),
523                     pRealOut);
524             }
525             else if(bComplex1 == false && bComplex2 == true)
526             {
527                 iMultiRealScalarByComplexMatrix(
528                     _pDouble->getReal(0,0),
529                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
530                     pRealOut, pImgOut);
531             }
532             else if(bComplex1 == true && bComplex2 == false)
533             {
534                 iMultiComplexScalarByRealMatrix(
535                     _pDouble->getReal(0,0), _pDouble->getImg(0,0),
536                     pRealIn, 1, pPolyIn->getRank(),
537                     pRealOut, pImgOut);
538             }
539             else if(bComplex1 == true && bComplex2 == true)
540             {
541                 iMultiComplexScalarByComplexMatrix(
542                     _pDouble->getReal(0,0), _pDouble->getImg(0,0),
543                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
544                     pRealOut, pImgOut);
545             }
546         }
547     }
548     else if(bScalar2)
549     {
550         double *pDoubleR    = _pDouble->getReal();
551         double *pDoubleI    = _pDouble->getImg();
552
553         SinglePoly *pPolyIn = _pPoly->get(0);
554         double* pRealIn     = pPolyIn->getCoef()->getReal();
555         double* pImgIn      = pPolyIn->getCoef()->getImg();
556
557         for(int i = 0 ; i < _pDouble->getSize() ; i++)
558         {
559             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
560             double* pRealOut        = pPolyOut->getCoef()->getReal();
561             double* pImgOut         = pPolyOut->getCoef()->getImg();
562
563             if(bComplex1 == false && bComplex2 == false)
564             {
565                 iMultiRealScalarByRealMatrix(
566                     pDoubleR[i],
567                     pRealIn, 1, pPolyIn->getRank(),
568                     pRealOut);
569             }
570             else if(bComplex1 == false && bComplex2 == true)
571             {
572                 iMultiRealScalarByComplexMatrix(
573                     pDoubleR[i],
574                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
575                     pRealOut, pImgOut);
576             }
577             else if(bComplex1 == true && bComplex2 == false)
578             {
579                 iMultiComplexScalarByRealMatrix(
580                     pDoubleR[i], pDoubleI[i],
581                     pRealIn, 1, pPolyIn->getRank(),
582                     pRealOut, pImgOut);
583             }
584             else if(bComplex1 == true && bComplex2 == true)
585             {
586                 iMultiComplexScalarByComplexMatrix(
587                     pDoubleR[i], pDoubleI[i],
588                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
589                     pRealOut, pImgOut);
590             }
591         }
592     }
593     else if(_pDouble->getCols() == _pPoly->getRows())
594     {
595         Double *pCoef = _pPoly->getCoef();
596         Double *pTemp = new Double(_pDouble->getRows(), pCoef->getCols(), bComplex1 || bComplex2);
597
598         if(bComplex1 == false && bComplex2 == false)
599         {
600             iMultiRealMatrixByRealMatrix(
601                 _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
602                 pCoef->getReal(), pCoef->getRows(), pCoef->getCols(),
603                 pTemp->getReal());
604         }
605         else if(bComplex1 == false && bComplex2 == true)
606         {
607             iMultiRealMatrixByComplexMatrix(
608                 _pDouble->getReal(), _pDouble->getRows(), _pDouble->getCols(),
609                 pCoef->getReal(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
610                 pTemp->getReal(), pTemp->getImg());
611
612         }
613         else if(bComplex1 == true && bComplex2 == false)
614         {
615             iMultiComplexMatrixByRealMatrix(
616                 _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
617                 pCoef->getReal(), pCoef->getRows(), pCoef->getCols(),
618                 pTemp->getReal(), pTemp->getImg());
619         }
620         else //if(bComplex1 == true && bComplex2 == true)
621         {
622             iMultiComplexMatrixByComplexMatrix(
623                 _pDouble->getReal(), _pDouble->getImg(), _pDouble->getRows(), _pDouble->getCols(),
624                 pCoef->getReal(), pCoef->getImg(), pCoef->getRows(), pCoef->getCols(),
625                 pTemp->getReal(), pTemp->getImg());
626         }
627
628         (*_pPolyOut)->setCoef(pTemp);
629         delete pTemp;
630     }
631     else
632     {
633         return 1;
634     }
635
636     return 0;
637 }
638
639 int MultiplyPolyByDouble(Polynom* _pPoly, Double* _pDouble, Polynom **_pPolyOut)
640 {
641     bool bComplex1  = _pPoly->isComplex();
642     bool bComplex2  = _pDouble->isComplex();
643     bool bScalar1   = _pPoly->getRows() == 1 && _pPoly->getCols() == 1;
644     bool bScalar2   = _pDouble->getRows() == 1 && _pDouble->getCols() == 1;
645
646     int iRowResult  = 0;
647     int iColResult  = 0;
648     int *piRank     = NULL;
649
650     if(bScalar1)
651     {
652         iRowResult = _pDouble->getRows();
653         iColResult = _pDouble->getCols();
654
655         piRank = new int[iRowResult * iColResult];
656         for(int i = 0 ; i < iRowResult * iColResult ; i++)
657         {
658             piRank[i] = _pPoly->get(0)->getRank();
659         }
660     }
661     else if (bScalar2)
662     {
663         iRowResult = _pPoly->getRows();
664         iColResult = _pPoly->getCols();
665
666         piRank = new int[iRowResult * iColResult];
667         for(int i = 0 ; i < iRowResult * iColResult ; i++)
668         {
669             piRank[i] = _pPoly->get(i)->getRank();
670         }
671     }
672     else if(_pPoly->getCols() == _pDouble->getRows())
673     {
674         iRowResult = _pPoly->getRows();
675         iColResult = _pDouble->getCols();
676         piRank = new int[iRowResult * iColResult];
677         for(int i = 0 ; i < iRowResult * iColResult ; i++)
678         {
679             piRank[i] = _pPoly->getMaxRank();
680         }
681     }
682     else
683     {
684         return 1;
685     }
686
687     (*_pPolyOut) = new Polynom(_pPoly->getVariableName(), iRowResult, iColResult, piRank);
688     delete[] piRank;
689     if(bComplex1 || bComplex2)
690     {
691         (*_pPolyOut)->setComplex(true);
692     }
693     //Output variables
694     //Polynom *pResult  = NULL;
695     bool bComplexOut = bComplex1 || bComplex2;
696
697     if(bScalar1)
698     {//p[1] * A
699         double *pDoubleR    = _pDouble->getReal();
700         double *pDoubleI    = _pDouble->getImg();
701
702         SinglePoly *pPolyIn = _pPoly->get(0);
703         double* pRealIn     = pPolyIn->getCoef()->getReal();
704         double* pImgIn      = pPolyIn->getCoef()->getImg();
705
706         for(int i = 0 ; i < _pDouble->getSize() ; i++)
707         {
708             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
709             double* pRealOut        = pPolyOut->getCoef()->getReal();
710             double* pImgOut         = pPolyOut->getCoef()->getImg();
711
712             if(bComplex1 == false && bComplex2 == false)
713             {
714                 iMultiRealScalarByRealMatrix(
715                     pDoubleR[i],
716                     pRealIn, 1, pPolyIn->getRank(),
717                     pRealOut);
718             }
719             else if(bComplex1 == false && bComplex2 == true)
720             {
721                 iMultiComplexScalarByRealMatrix(
722                     pDoubleR[i], pDoubleI[i],
723                     pRealIn, 1, pPolyIn->getRank(),
724                     pRealOut, pImgOut);
725             }
726             else if(bComplex1 == true && bComplex2 == false)
727             {
728                 iMultiRealScalarByComplexMatrix(
729                     pDoubleR[i],
730                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
731                     pRealOut, pImgOut);
732             }
733             else if(bComplex1 == true && bComplex2 == true)
734             {
735                 iMultiComplexScalarByComplexMatrix(
736                     pDoubleR[i], pDoubleI[i],
737                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
738                     pRealOut, pImgOut);
739             }
740         }
741     }
742     else if(bScalar2)
743     {//p * cst
744         for(int i = 0 ; i < _pPoly->getSize() ; i++)
745         {
746             SinglePoly *pPolyIn = _pPoly->get(i);
747             double* pRealIn     = pPolyIn->getCoef()->getReal();
748             double* pImgIn      = pPolyIn->getCoef()->getImg();
749
750             SinglePoly *pPolyOut    = (*_pPolyOut)->get(i);
751             double* pRealOut        = pPolyOut->getCoef()->getReal();
752             double* pImgOut         = pPolyOut->getCoef()->getImg();
753
754             if(bComplex1 == false && bComplex2 == false)
755             {
756                 iMultiRealScalarByRealMatrix(
757                     _pDouble->getReal(0,0),
758                     pRealIn, 1, pPolyIn->getRank(),
759                     pRealOut);
760             }
761             else if(bComplex1 == false && bComplex2 == true)
762             {
763                 iMultiComplexScalarByRealMatrix(
764                     _pDouble->getReal(0,0), _pDouble->getImg(0,0),
765                     pRealIn, 1, pPolyIn->getRank(),
766                     pRealOut, pImgOut);
767             }
768             else if(bComplex1 == true && bComplex2 == false)
769             {
770                 iMultiRealScalarByComplexMatrix(
771                     _pDouble->getReal(0,0),
772                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
773                     pRealOut, pImgOut);
774             }
775             else if(bComplex1 == true && bComplex2 == true)
776             {
777                 iMultiComplexScalarByComplexMatrix(
778                     _pDouble->getReal(0,0), _pDouble->getImg(0,0),
779                     pRealIn, pImgIn, 1, pPolyIn->getRank(),
780                     pRealOut, pImgOut);
781             }
782         }
783     }
784     else if(_pPoly->getCols() == _pDouble->getRows())
785     {
786         //Distribution a la mano par appels a des sous-fonctions ( iMulti...ScalarBy...Scalar ) plus iAdd...To... )
787
788         //for each line of _pPoly
789         for(int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
790         {//for each col of _pDouble
791             for(int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
792             {//for each rows of _pDouble / cols of _pPoly
793                 for(int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
794                 {
795                     Double *pPolyCoef = _pPoly->get(iRow1, iRow2)->getCoef();
796
797                     Double *pDouble = NULL;
798                     if(bComplex2 == false)
799                     {
800                         pDouble = new Double(_pDouble->getReal(iRow2, iCol2));
801                     }
802                     else
803                     {
804                         pDouble = new Double(_pDouble->getReal(iRow2, iCol2), _pDouble->getImg(iRow2, iCol2));
805                     }
806
807                     Double *TimeDouble = NULL; //(pPolyCoef->getRows(), pPolyCoef->getCols(), bComplexOut);
808                     MultiplyDoubleByDouble(pPolyCoef, pDouble, &TimeDouble);
809                     Double *pAddDouble = NULL;
810
811                     //Adjust size to allow vector multiplication
812                     Double* pCoef = (*_pPolyOut)->get(iRow1, iCol2)->getCoef();
813
814                     if(TimeDouble->getRows() > pCoef->getRows())
815                     {
816                         pCoef->resize(TimeDouble->getRows(), pCoef->getCols());
817                     }
818                     else if(TimeDouble->getRows() < pCoef->getRows())
819                     {
820                         TimeDouble->resize(pCoef->getRows(), TimeDouble->getCols());
821                     }
822
823                     if(TimeDouble->getCols() > pCoef->getCols())
824                     {
825                         pCoef->resize(pCoef->getRows(), TimeDouble->getCols());
826                     }
827                     else if(TimeDouble->getCols() < pCoef->getCols())
828                     {
829                         TimeDouble->resize(TimeDouble->getRows(), pCoef->getCols());
830                     }
831
832                     AddDoubleToDouble(TimeDouble, pCoef, &pAddDouble);
833                     (*_pPolyOut)->setCoef(iRow1, iCol2, pAddDouble);
834
835                     delete pAddDouble;
836                     delete pDouble;
837                 }//for(int iRow2 = 0 ; iRow2 < _pDouble->getRows() ; iRow2++)
838             }//for(int iCol2 = 0 ; iCol2 < _pDouble->getCols() ; iCol2++)
839         }//for(int iRow1 = 0 ; iRow1 < _pPoly->getRows() ; iRow1++)
840     }
841     else //Wrong dimensions.
842     {
843         return 1;
844     }
845     return 0;
846 }
847
848 int MultiplyPolyByPoly(Polynom* _pPoly1, Polynom* _pPoly2, Polynom** _pPolyOut)
849 {
850     bool bComplex1      = _pPoly1->isComplex();
851     bool bComplex2      = _pPoly2->isComplex();
852     bool bScalar1               = _pPoly1->getRows() == 1 && _pPoly1->getCols() == 1;
853     bool bScalar2               = _pPoly2->getRows() == 1 && _pPoly2->getCols() == 1;
854
855     int iRowResult      = 0;
856     int iColResult      = 0;
857     int *piRank                 = NULL;
858
859     if(_pPoly1->getSize() == 1 && _pPoly2->getSize() == 1)
860     {
861         iRowResult = 1;
862         iColResult = 1;
863
864         piRank = new int[1];
865         piRank[0] = _pPoly1->get(0)->getRank() + _pPoly2->get(0)->getRank() - 1;
866     }
867     else if(_pPoly1->getSize() == 1)
868     {
869         iRowResult = _pPoly2->getRows();
870         iColResult = _pPoly2->getCols();
871
872         piRank = new int[iRowResult * iColResult];
873         for(int i = 0 ; i < iRowResult * iColResult ; i++)
874         {
875             piRank[i] = _pPoly1->get(0)->getRank() + _pPoly2->get(i)->getRank() - 1;
876         }
877     }
878     else if (_pPoly2->getSize() == 1)
879     {
880         iRowResult = _pPoly1->getRows();
881         iColResult = _pPoly1->getCols();
882
883         piRank = new int[iRowResult * iColResult];
884         for(int i = 0 ; i < iRowResult * iColResult ; i++)
885         {
886             piRank[i] = _pPoly2->get(0)->getRank() * _pPoly1->get(i)->getRank();
887         }
888     }
889     else if(_pPoly1->getCols() == _pPoly2->getRows())
890     {
891         iRowResult = _pPoly1->getRows();
892         iColResult = _pPoly2->getCols();
893         piRank = new int[iRowResult * iColResult];
894         for(int i = 0 ; i < iRowResult * iColResult ; i++)
895         {
896             piRank[i] = _pPoly1->getMaxRank() * _pPoly2->getMaxRank();
897         }
898     }
899     else
900     {
901         return 1;
902     }
903
904     (*_pPolyOut) = new Polynom(_pPoly1->getVariableName(), iRowResult, iColResult, piRank);
905
906     delete[] piRank;
907
908     if(bComplex1 || bComplex2)
909     {
910         (*_pPolyOut)->setComplex(true);
911     }
912
913     if(bScalar1 && bScalar2)
914     {//poly1(0) * poly2(0)
915         if(bComplex1 == false && bComplex2 == false)
916         {
917             SinglePoly *pPoly1          = _pPoly1->get(0);
918             SinglePoly *pPoly2          = _pPoly2->get(0);
919             SinglePoly *pPolyOut        = (*_pPolyOut)->get(0);
920
921             pPolyOut->getCoef()->setZeros();
922
923             iMultiRealPolyByRealPoly(
924                 pPoly1->getCoef()->getReal(), pPoly1->getRank(),
925                 pPoly2->getCoef()->getReal(), pPoly2->getRank(),
926                 pPolyOut->getCoef()->getReal(), pPolyOut->getRank());
927         }
928         else if(bComplex1 == false && bComplex2 == true)
929         {
930             SinglePoly *pPoly1          = _pPoly1->get(0);
931             SinglePoly *pPoly2          = _pPoly2->get(0);
932             SinglePoly *pPolyOut        = (*_pPolyOut)->get(0);
933
934             pPolyOut->getCoef()->setZeros();
935
936             iMultiRealPolyByComplexPoly(
937                 pPoly1->getCoef()->getReal(), pPoly1->getRank(),
938                 pPoly2->getCoef()->getReal(), pPoly2->getCoef()->getImg(), pPoly2->getRank(),
939                 pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
940         }
941         else if(bComplex1 == true && bComplex2 == false)
942         {
943             SinglePoly *pPoly1          = _pPoly1->get(0);
944             SinglePoly *pPoly2          = _pPoly2->get(0);
945             SinglePoly *pPolyOut        = (*_pPolyOut)->get(0);
946
947             pPolyOut->getCoef()->setZeros();
948
949             iMultiComplexPolyByRealPoly(
950                 pPoly1->getCoef()->getReal(), pPoly1->getCoef()->getImg(), pPoly1->getRank(),
951                 pPoly2->getCoef()->getReal(), pPoly2->getRank(),
952                 pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
953         }
954         else if(bComplex1 == true && bComplex2 == true)
955         {
956             SinglePoly *pPoly1                  = _pPoly1->get(0);
957             SinglePoly *pPoly2                  = _pPoly2->get(0);
958             SinglePoly *pPolyOut                = (*_pPolyOut)->get(0);
959             Double *pCoef1              = pPoly1->getCoef();
960             Double *pCoef2              = pPoly2->getCoef();
961             Double *pCoefOut    = pPolyOut->getCoef();
962
963             pCoefOut->setZeros();
964             iMultiComplexPolyByComplexPoly(
965                 pCoef1->getReal(), pCoef1->getImg(), pPoly1->getRank(),
966                 pCoef2->getReal(), pCoef2->getImg(), pPoly2->getRank(),
967                 pCoefOut->getReal(), pCoefOut->getImg(), pPolyOut->getRank());
968         }
969     }
970     else if(bScalar1)
971     {//poly1(0) * poly2(n)
972         SinglePoly *pPoly1              = _pPoly1->get(0);
973         if(bComplex1 == false && bComplex2 == false)
974         {
975             for(int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
976             {
977                 SinglePoly *pPoly2              = _pPoly2->get(iPoly);
978                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
979
980                 pPolyOut->getCoef()->setZeros();
981
982                 iMultiRealPolyByRealPoly(
983                     pPoly1->getCoef()->getReal(), pPoly1->getRank(),
984                     pPoly2->getCoef()->getReal(), pPoly2->getRank(),
985                     pPolyOut->getCoef()->getReal(), pPolyOut->getRank());
986             }
987         }
988         else if(bComplex1 == false && bComplex2 == true)
989         {
990             for(int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
991             {
992                 SinglePoly *pPoly2              = _pPoly2->get(iPoly);
993                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
994
995                 pPolyOut->getCoef()->setZeros();
996
997                 iMultiRealPolyByComplexPoly(
998                     pPoly1->getCoef()->getReal(), pPoly1->getRank(),
999                     pPoly2->getCoef()->getReal(), pPoly2->getCoef()->getImg(), pPoly2->getRank(),
1000                     pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
1001             }
1002         }
1003         else if(bComplex1 == true && bComplex2 == false)
1004         {
1005             for(int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
1006             {
1007                 SinglePoly *pPoly2              = _pPoly2->get(iPoly);
1008                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
1009
1010                 pPolyOut->getCoef()->setZeros();
1011
1012                 iMultiComplexPolyByRealPoly(
1013                     pPoly1->getCoef()->getReal(), pPoly1->getCoef()->getImg(), pPoly1->getRank(),
1014                     pPoly2->getCoef()->getReal(), pPoly2->getRank(),
1015                     pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
1016             }
1017         }
1018         else if(bComplex1 == true && bComplex2 == true)
1019         {
1020             Double *pCoef1              = pPoly1->getCoef();
1021             for(int iPoly = 0 ; iPoly < _pPoly2->getSize() ; iPoly++)
1022             {
1023                 SinglePoly *pPoly2                      = _pPoly2->get(iPoly);
1024                 SinglePoly *pPolyOut            = (*_pPolyOut)->get(iPoly);
1025                 Double *pCoef2          = pPoly2->getCoef();
1026                 Double *pCoefOut        = pPolyOut->getCoef();
1027
1028
1029                 pCoefOut->setZeros();
1030
1031                 iMultiComplexPolyByComplexPoly(
1032                     pCoef1->getReal(), pCoef1->getImg(), pPoly1->getRank(),
1033                     pCoef2->getReal(), pCoef2->getImg(), pPoly2->getRank(),
1034                     pCoefOut->getReal(), pCoefOut->getImg(), pPolyOut->getRank());
1035             }
1036         }
1037     }
1038     else if(bScalar2)
1039     {//poly1(n) * poly2(0)
1040         SinglePoly *pPoly2              = _pPoly2->get(0);
1041         if(bComplex1 == false && bComplex2 == false)
1042         {
1043             for(int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
1044             {
1045                 SinglePoly *pPoly1              = _pPoly1->get(iPoly);
1046                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
1047
1048                 pPolyOut->getCoef()->setZeros();
1049
1050                 iMultiRealPolyByRealPoly(
1051                     pPoly1->getCoef()->getReal(), pPoly1->getRank(),
1052                     pPoly2->getCoef()->getReal(), pPoly2->getRank(),
1053                     pPolyOut->getCoef()->getReal(), pPolyOut->getRank());
1054             }
1055         }
1056         else if(bComplex1 == false && bComplex2 == true)
1057         {
1058             for(int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
1059             {
1060                 SinglePoly *pPoly1              = _pPoly1->get(iPoly);
1061                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
1062
1063                 pPolyOut->getCoef()->setZeros();
1064
1065                 iMultiRealPolyByComplexPoly(
1066                     pPoly1->getCoef()->getReal(), pPoly1->getRank(),
1067                     pPoly2->getCoef()->getReal(), pPoly2->getCoef()->getImg(), pPoly2->getRank(),
1068                     pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
1069             }
1070         }
1071         else if(bComplex1 == true && bComplex2 == false)
1072         {
1073             for(int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
1074             {
1075                 SinglePoly *pPoly1              = _pPoly1->get(iPoly);
1076                 SinglePoly *pPolyOut    = (*_pPolyOut)->get(iPoly);
1077
1078                 pPolyOut->getCoef()->setZeros();
1079
1080                 iMultiComplexPolyByRealPoly(
1081                     pPoly1->getCoef()->getReal(), pPoly1->getCoef()->getImg(), pPoly1->getRank(),
1082                     pPoly2->getCoef()->getReal(), pPoly2->getRank(),
1083                     pPolyOut->getCoef()->getReal(), pPolyOut->getCoef()->getImg(), pPolyOut->getRank());
1084             }
1085         }
1086         else if(bComplex1 == true && bComplex2 == true)
1087         {
1088             Double *pCoef2              = pPoly2->getCoef();
1089             for(int iPoly = 0 ; iPoly < _pPoly1->getSize() ; iPoly++)
1090             {
1091                 SinglePoly *pPoly1                      = _pPoly1->get(iPoly);
1092                 SinglePoly *pPolyOut            = (*_pPolyOut)->get(iPoly);
1093                 Double *pCoef1          = pPoly1->getCoef();
1094                 Double *pCoefOut        = pPolyOut->getCoef();
1095
1096
1097                 pCoefOut->setZeros();
1098
1099                 iMultiComplexPolyByComplexPoly(
1100                     pCoef1->getReal(), pCoef1->getImg(), pPoly1->getRank(),
1101                     pCoef2->getReal(), pCoef2->getImg(), pPoly2->getRank(),
1102                     pCoefOut->getReal(), pCoefOut->getImg(), pPolyOut->getRank());
1103             }
1104         }
1105     }
1106     else
1107     {// matrix by matrix
1108         if(bComplex1 == false && bComplex2 == false)
1109         {
1110             double *pReal       = NULL;
1111             SinglePoly *pTemp           = new SinglePoly(&pReal, (*_pPolyOut)->getMaxRank());
1112
1113             for(int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1114             {
1115                 for(int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1116                 {
1117                     SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1118                     pResult->getCoef()->setZeros();
1119
1120                     for(int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1121                     {
1122                         SinglePoly *pL                  = _pPoly1->get(iRow, iCommon);
1123                         SinglePoly *pR                  = _pPoly2->get(iCommon, iCol);
1124
1125                         pTemp->getCoef()->setZeros();
1126
1127                         iMultiRealPolyByRealPoly(
1128                             pL->getCoef()->getReal(), pL->getRank(),
1129                             pR->getCoef()->getReal(), pR->getRank(),
1130                             pTemp->getCoef()->getReal(), pL->getRank() + pR->getRank() - 1);
1131
1132                         iAddRealPolyToRealPoly(
1133                             pResult->getCoef()->getReal(), pResult->getRank(),
1134                             pTemp->getCoef()->getReal(), pResult->getRank(),
1135                             pResult->getCoef()->getReal(), pResult->getRank());
1136                     }
1137                 }
1138             }
1139         }
1140         else if(bComplex1 == false && bComplex2 == true)
1141         {
1142             double *pReal       = NULL;
1143             double *pImg        = NULL;
1144             SinglePoly *pTemp           = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1145
1146             for(int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1147             {
1148                 for(int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1149                 {
1150                     SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1151                     pResult->getCoef()->setZeros();
1152
1153                     for(int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1154                     {
1155                         SinglePoly *pL                  = _pPoly1->get(iRow, iCommon);
1156                         SinglePoly *pR                  = _pPoly2->get(iCommon, iCol);
1157
1158                         pTemp->getCoef()->setZeros();
1159
1160                         iMultiRealPolyByComplexPoly(
1161                             pL->getCoef()->getReal(), pL->getRank(),
1162                             pR->getCoef()->getReal(), pR->getCoef()->getImg(), pR->getRank(),
1163                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pL->getRank() + pR->getRank() - 1);
1164
1165                         iAddComplexPolyToComplexPoly(
1166                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank(),
1167                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pResult->getRank(),
1168                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank());
1169                     }
1170                 }
1171             }
1172         }
1173         else if(bComplex1 == true && bComplex2 == false)
1174         {
1175             double *pReal       = NULL;
1176             double *pImg        = NULL;
1177             SinglePoly *pTemp           = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1178
1179             for(int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1180             {
1181                 for(int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1182                 {
1183                     SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1184                     pResult->getCoef()->setZeros();
1185
1186                     for(int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1187                     {
1188                         SinglePoly *pL                  = _pPoly1->get(iRow, iCommon);
1189                         SinglePoly *pR                  = _pPoly2->get(iCommon, iCol);
1190
1191                         pTemp->getCoef()->setZeros();
1192
1193                         iMultiRealPolyByComplexPoly(
1194                             pR->getCoef()->getReal(), pR->getRank(),
1195                             pL->getCoef()->getReal(), pL->getCoef()->getImg(), pL->getRank(),
1196                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pL->getRank() + pR->getRank() - 1);
1197
1198                         iAddComplexPolyToComplexPoly(
1199                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank(),
1200                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pResult->getRank(),
1201                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank());
1202                     }
1203                 }
1204             }
1205         }
1206         else if(bComplex1 == true && bComplex2 == true)
1207         {
1208             double *pReal       = NULL;
1209             double *pImg        = NULL;
1210             SinglePoly *pTemp           = new SinglePoly(&pReal, &pImg, (*_pPolyOut)->getMaxRank());
1211
1212             for(int iRow = 0 ; iRow < _pPoly1->getRows() ; iRow++)
1213             {
1214                 for(int iCol = 0 ; iCol < _pPoly2->getCols() ; iCol++)
1215                 {
1216                     SinglePoly *pResult = (*_pPolyOut)->get(iRow, iCol);
1217                     pResult->getCoef()->setZeros();
1218
1219                     for(int iCommon = 0 ; iCommon < _pPoly1->getCols() ; iCommon++)
1220                     {
1221                         SinglePoly *pL                  = _pPoly1->get(iRow, iCommon);
1222                         SinglePoly *pR                  = _pPoly2->get(iCommon, iCol);
1223
1224                         pTemp->getCoef()->setZeros();
1225
1226                         iMultiComplexPolyByComplexPoly(
1227                             pL->getCoef()->getReal(), pL->getCoef()->getImg(), pL->getRank(),
1228                             pR->getCoef()->getReal(), pR->getCoef()->getImg(), pR->getRank(),
1229                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pL->getRank() + pR->getRank() - 1);
1230
1231                         iAddComplexPolyToComplexPoly(
1232                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank(),
1233                             pTemp->getCoef()->getReal(), pTemp->getCoef()->getImg(), pResult->getRank(),
1234                             pResult->getCoef()->getReal(), pResult->getCoef()->getImg(), pResult->getRank());
1235                     }
1236                 }
1237             }
1238         }
1239     }
1240     return 0; //No Error;
1241 }
1242
1243 int MultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pSparseOut)
1244 {
1245     if(_pSparse1->getRows() == 1 && _pSparse1->getCols() == 1)
1246     {//scalar * sp
1247         Double* pDbl = NULL;
1248         if(_pSparse1->isComplex())
1249         {
1250             std::complex<double> dbl = _pSparse1->getImg(0,0);
1251             pDbl = new Double(dbl.real(), dbl.imag());
1252         }
1253         else
1254         {
1255             pDbl = new Double(_pSparse1->get(0,0));
1256         }
1257         
1258         MultiplyDoubleBySparse(pDbl, _pSparse2, (GenericType**)_pSparseOut);
1259         delete pDbl;
1260         return 0;
1261     }
1262
1263     if(_pSparse2->getRows() == 1 && _pSparse2->getCols() == 1)
1264     {//sp * scalar
1265         Double* pDbl = NULL;
1266         if(_pSparse2->isComplex())
1267         {
1268             std::complex<double> dbl = _pSparse2->getImg(0,0);
1269             pDbl = new Double(dbl.real(), dbl.imag());
1270         }
1271         else
1272         {
1273             pDbl = new Double(_pSparse2->get(0,0));
1274         }
1275         
1276         MultiplySparseByDouble(_pSparse1, pDbl, (GenericType**)_pSparseOut);
1277         delete pDbl;
1278         return 0;
1279     }
1280
1281     if(_pSparse1->getCols() != _pSparse2->getRows())
1282     {
1283         return 1;
1284     }
1285
1286     *_pSparseOut = _pSparse1->multiply(*_pSparse2);
1287     return 0;
1288 }
1289
1290 int MultiplyDoubleBySparse(Double* _pDouble, Sparse *_pSparse, GenericType** _pOut)
1291 {//D * SP
1292     if(_pDouble->isEmpty())
1293     {
1294         *_pOut = Double::Empty();
1295         return 0;
1296     }
1297
1298     if(_pDouble->isScalar())
1299     {//d * SP -> SP
1300         Sparse* pOut = NULL;
1301         if(_pDouble->isComplex())
1302         {
1303             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1304             pOut = _pSparse->multiply(dbl);
1305         }
1306         else
1307         {
1308             pOut = _pSparse->multiply(_pDouble->get(0));
1309         }
1310         *_pOut = pOut;
1311         return 0;
1312     }
1313
1314     if(_pSparse->isScalar())
1315     {//D * sp -> D .* d
1316         Double* pD = NULL;
1317         
1318         if(_pSparse->isComplex())
1319         {
1320             std::complex<double> dbl(_pSparse->getImg(0, 0));
1321             pD = new Double(dbl.real(), dbl.imag());
1322         }
1323         else
1324         {
1325             pD = new Double(_pSparse->get(0, 0));
1326         }
1327
1328         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1329         *_pOut = pIT->getAs<GenericType>();
1330         delete pD;
1331         return 0;
1332     }
1333
1334     if(_pDouble->getCols() != _pSparse->getRows())
1335     {
1336         return 1;
1337     }
1338
1339     //try to be smart and only compute for non zero values
1340
1341     //get some information
1342     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1343     double* pRows = new double[iNonZeros * 2];
1344     _pSparse->outputRowCol(pRows);
1345     double* pCols = pRows + iNonZeros;
1346     double* pValR = new double[iNonZeros];
1347     double* pValI = new double[iNonZeros];
1348     _pSparse->outputValues(pValR, pValI);
1349
1350     Double* pOut = new Double(_pDouble->getRows(), _pSparse->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1351     pOut->setZeros();
1352
1353     if(_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1354     {
1355         for(int i = 0 ; i < iNonZeros ; i++)
1356         {
1357             int iRow = static_cast<int>(pRows[i]) - 1;
1358             int iCol = static_cast<int>(pCols[i]) - 1;
1359             double dbl = pValR[i];
1360
1361             for(int j = 0 ; j < _pDouble->getRows() ; j++)
1362             {
1363                 double dblVal = _pDouble->get(j, iRow) * dbl;
1364                 pOut->set(j, iCol, pOut->get(j, iCol) + dblVal);
1365             }
1366         }
1367     }
1368     else if(_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1369     {//a * (b ci) -> ab ac
1370         for(int i = 0 ; i < iNonZeros ; i++)
1371         {
1372             int iRow = static_cast<int>(pRows[i]) - 1;
1373             int iCol = static_cast<int>(pCols[i]) - 1;
1374             double dblR = pValR[i];
1375             double dblI = pValI[i];
1376
1377             for(int j = 0 ; j < _pDouble->getRows() ; j++)
1378             {
1379                 double dblValR = _pDouble->get(j, iRow) * dblR;
1380                 double dblValI = _pDouble->get(j, iRow) * dblI;
1381                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1382                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1383             }
1384         }
1385     }
1386     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1387     {//(a bi) * c -> ac + bc
1388         for(int i = 0 ; i < iNonZeros ; i++)
1389         {
1390             int iRow = static_cast<int>(pRows[i]) - 1;
1391             int iCol = static_cast<int>(pCols[i]) - 1;
1392             double dblR = pValR[i];
1393
1394             for(int j = 0 ; j < _pDouble->getRows() ; j++)
1395             {
1396                 double dblValR = _pDouble->get(j, iRow) * dblR;
1397                 double dblValI = _pDouble->getImg(j, iRow) * dblR;
1398                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1399                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1400             }
1401         }
1402     }
1403     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1404     {
1405         for(int i = 0 ; i < iNonZeros ; i++)
1406         {
1407             int iRow = static_cast<int>(pRows[i]) - 1;
1408             int iCol = static_cast<int>(pCols[i]) - 1;
1409             double dblR = pValR[i];
1410             double dblI = pValI[i];
1411
1412             for(int j = 0 ; j < _pDouble->getRows() ; j++)
1413             {
1414                 double dblValR = _pDouble->get(j, iRow) * dblR - _pDouble->getImg(j, iRow) * dblI;
1415                 double dblValI = _pDouble->get(j, iRow) * dblI + _pDouble->getImg(j, iRow) * dblR;
1416                 pOut->set(j, iCol, pOut->get(j, iCol) + dblValR);
1417                 pOut->setImg(j, iCol, pOut->getImg(j, iCol) + dblValI);
1418             }
1419         }
1420     }
1421     
1422     *_pOut = pOut;
1423     delete[] pRows;
1424     delete[] pValR;
1425     delete[] pValI;
1426
1427     return 0;
1428 }
1429
1430 int MultiplySparseByDouble(Sparse *_pSparse, Double*_pDouble, GenericType** _pOut)
1431 {
1432     if(_pDouble->isEmpty())
1433     {//SP * [] -> []
1434         *_pOut = Double::Empty();
1435         return 0;
1436     }
1437
1438     if(_pDouble->isScalar())
1439     {//SP * d -> SP
1440         Sparse* pOut = NULL;
1441         if(_pDouble->isComplex())
1442         {
1443             std::complex<double> dbl(_pDouble->get(0), _pDouble->getImg(0));
1444             pOut = _pSparse->multiply(dbl);
1445         }
1446         else
1447         {
1448             pOut = _pSparse->multiply(_pDouble->get(0));
1449         }
1450         *_pOut = pOut;
1451         return 0;
1452     }
1453
1454     if(_pSparse->isScalar())
1455     {//D * sp -> D .* d
1456         Double* pD = NULL;
1457         
1458         if(_pSparse->isComplex())
1459         {
1460             std::complex<double> dbl(_pSparse->getImg(0, 0));
1461             pD = new Double(dbl.real(), dbl.imag());
1462         }
1463         else
1464         {
1465             pD = new Double(_pSparse->get(0, 0));
1466         }
1467
1468         InternalType* pIT = GenericDotTimes(_pDouble, pD);
1469         *_pOut = pIT->getAs<GenericType>();
1470         delete pD;
1471         return 0;
1472     }
1473
1474     if(_pSparse->getCols() != _pDouble->getRows())
1475     {
1476         return 1;
1477     }
1478
1479     //try to be smart and only compute for non zero values
1480
1481     //get some information
1482     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1483     double* pRows = new double[iNonZeros * 2];
1484     _pSparse->outputRowCol(pRows);
1485     double* pCols = pRows + iNonZeros;
1486     double* pValR = new double[iNonZeros];
1487     double* pValI = new double[iNonZeros];
1488     _pSparse->outputValues(pValR, pValI);
1489
1490     Double* pOut = new Double(_pSparse->getRows(), _pDouble->getCols(), _pDouble->isComplex() | _pSparse->isComplex());
1491     pOut->setZeros();
1492
1493     if(_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1494     {
1495         for(int i = 0 ; i < iNonZeros ; i++)
1496         {
1497             int iRow    = static_cast<int>(pRows[i]) - 1;
1498             int iCol    = static_cast<int>(pCols[i]) - 1;
1499             double dbl  = pValR[i];
1500
1501             for(int j = 0 ; j < _pDouble->getCols() ; j++)
1502             {
1503                 double dblVal = _pDouble->get(iCol, j) * dbl;
1504                 pOut->set(iRow, j, pOut->get(iRow, j) + dblVal);
1505             }
1506         }
1507     }
1508     else if(_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1509     {//a * (b ci) -> ab ac
1510         for(int i = 0 ; i < iNonZeros ; i++)
1511         {
1512             int iRow = static_cast<int>(pRows[i]) - 1;
1513             int iCol = static_cast<int>(pCols[i]) - 1;
1514             double dblR = pValR[i];
1515             double dblI = pValI[i];
1516
1517             for(int j = 0 ; j < _pDouble->getCols() ; j++)
1518             {
1519                 double dblValR = _pDouble->get(iCol, j) * dblR;
1520                 double dblValI = _pDouble->get(iCol, j) * dblI;
1521                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1522                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1523             }
1524         }
1525     }
1526     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1527     {//(a bi) * c -> ac + bc
1528         for(int i = 0 ; i < iNonZeros ; i++)
1529         {
1530             int iRow = static_cast<int>(pRows[i]) - 1;
1531             int iCol = static_cast<int>(pCols[i]) - 1;
1532             double dblR = pValR[i];
1533
1534             for(int j = 0 ; j < _pDouble->getCols() ; j++)
1535             {
1536                 double dblValR = _pDouble->get(iCol, j) * dblR;
1537                 double dblValI = _pDouble->getImg(iCol, j) * dblR;
1538                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1539                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1540             }
1541         }
1542     }
1543     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1544     {
1545         for(int i = 0 ; i < iNonZeros ; i++)
1546         {
1547             int iRow = static_cast<int>(pRows[i]) - 1;
1548             int iCol = static_cast<int>(pCols[i]) - 1;
1549             double dblR = pValR[i];
1550             double dblI = pValI[i];
1551
1552             for(int j = 0 ; j < _pDouble->getCols() ; j++)
1553             {
1554                 double dblValR = _pDouble->get(iCol, j) * dblR - _pDouble->getImg(iCol, j) * dblI;
1555                 double dblValI = _pDouble->get(iCol, j) * dblI + _pDouble->getImg(iCol, j) * dblR;
1556                 pOut->set(iRow, j, pOut->get(iRow, j) + dblValR);
1557                 pOut->setImg(iRow, j, pOut->getImg(iRow, j) + dblValI);
1558             }
1559         }
1560     }
1561     
1562     *_pOut = pOut;
1563     delete[] pRows;
1564     delete[] pValR;
1565     delete[] pValI;
1566
1567     return 0;
1568 }
1569
1570 int DotMultiplySparseBySparse(Sparse* _pSparse1, Sparse* _pSparse2, Sparse** _pOut)
1571 {
1572     if(_pSparse1->isScalar() || _pSparse2->isScalar())
1573     {//SP .* sp or sp .* SP
1574         return MultiplySparseBySparse(_pSparse1, _pSparse2, _pOut);
1575     }
1576
1577     if(_pSparse1->getRows() != _pSparse2->getRows() || _pSparse1->getCols() != _pSparse2->getCols())
1578     {
1579         return 1;
1580     }
1581
1582     *_pOut = _pSparse1->dotMultiply(*_pSparse2);
1583
1584     return 0;
1585 }
1586
1587 int DotMultiplyDoubleBySparse(Double* _pDouble, Sparse* _pSparse, GenericType**  _pOut)
1588 {
1589     if(_pDouble->isEmpty())
1590     {
1591         *_pOut = Double::Empty();
1592         return 0;
1593     }
1594
1595     if(_pDouble->isScalar())
1596     {
1597         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1598     }
1599
1600     if(_pSparse->isScalar())
1601     {
1602         return MultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1603     }
1604
1605     if(_pSparse->getRows() != _pDouble->getRows() || _pSparse->getCols() != _pDouble->getCols())
1606     {
1607         return 1;
1608     }
1609
1610     Sparse* pOut = new Sparse(_pDouble->getRows(), _pDouble->getCols(), _pSparse->isComplex() || _pDouble->isComplex());
1611     //get some information
1612     int iNonZeros = static_cast<int>(_pSparse->nonZeros());
1613     double* pRows = new double[iNonZeros * 2];
1614     _pSparse->outputRowCol(pRows);
1615     double* pCols = pRows + iNonZeros;
1616
1617     if(_pDouble->isComplex() == false && _pSparse->isComplex() == false)
1618     {
1619         for(int i = 0 ; i < iNonZeros ; i++)
1620         {
1621             int iRow = static_cast<int>(pRows[i]) - 1;
1622             int iCol = static_cast<int>(pCols[i]) - 1;
1623             pOut->set(iRow, iCol, _pSparse->get(iRow, iCol) * _pDouble->get(iRow, iCol));
1624         }
1625     }
1626     else if(_pDouble->isComplex() == false && _pSparse->isComplex() == true)
1627     {
1628         for(int i = 0 ; i < iNonZeros ; i++)
1629         {
1630             int iRow = static_cast<int>(pRows[i]) - 1;
1631             int iCol = static_cast<int>(pCols[i]) - 1;
1632             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1633             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.imag() * _pDouble->get(iRow, iCol));
1634             pOut->set(iRow, iCol, newVal);
1635         }
1636     }
1637     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == false)
1638     {
1639         for(int i = 0 ; i < iNonZeros ; i++)
1640         {
1641             int iRow = static_cast<int>(pRows[i]) - 1;
1642             int iCol = static_cast<int>(pCols[i]) - 1;
1643             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1644             std::complex<double> newVal(dbl.real() * _pDouble->get(iRow, iCol), dbl.real() * _pDouble->getImg(iRow, iCol));
1645             pOut->set(iRow, iCol, newVal);
1646         }
1647     }
1648     else if(_pDouble->isComplex() == true && _pSparse->isComplex() == true)
1649     {
1650         for(int i = 0 ; i < iNonZeros ; i++)
1651         {
1652             int iRow = static_cast<int>(pRows[i]) - 1;
1653             int iCol = static_cast<int>(pCols[i]) - 1;
1654             std::complex<double> dbl = _pSparse->getImg(iRow, iCol);
1655             double dblR = _pDouble->get(iRow, iCol) * dbl.real() - _pDouble->getImg(iRow, iCol) * dbl.imag();
1656             double dblI = _pDouble->getImg(iRow, iCol) * dbl.real() + _pDouble->get(iRow, iCol) * dbl.imag();
1657
1658             std::complex<double> newVal(dblR, dblI);
1659             pOut->set(iRow, iCol, newVal);
1660         }
1661     }
1662     
1663     *_pOut = pOut;
1664     delete[] pRows;
1665
1666     return 0;
1667 }
1668
1669 int DotMultiplySparseByDouble(Sparse* _pSparse, Double* _pDouble, GenericType** _pOut)
1670 {
1671     return DotMultiplyDoubleBySparse(_pDouble, _pSparse, _pOut);
1672 }